(* Upadted for the latest version of HOL Light (JULY 2014) *)
(* ========================================================================= *)
(*               Formalization of Electromagnetic Optics                     *)
(*                                                                           *)
(*            (c) Copyright, Sanaz Khan Afshar & Vincent Aravantinos 2011-13 *)
(*                           Hardware Verification Group,                    *)
(*                           Concordia University                            *)
(*                                                                           *)
(*                Contact:   <s_khanaf@encs.concordia.ca>                    *)
(*                           <vincent@encs.concordia.ca>                     *)
(*                                                                           *)
(* This file deals with the definition and basic theorems about the          *)
(* electromagnetic model.                                                    *)
(* ========================================================================= *)

new_type_abbrev("time",`:real`);;

(********************************)
(* Electromagnetic fields       *)
(********************************)

new_type_abbrev("single_field",`:point -> time -> complex^3`);;
new_type_abbrev("emf", `:point -> time -> complex^3 # complex^3`);;

let e_of_emf = new_definition
  `e_of_emf (emf:emf) : single_field = \r t. let (e,h) = emf r t in e`;;
let h_of_emf = new_definition
  `h_of_emf (emf:emf) : single_field = \r t. let (e,h) = emf r t in h`;;
let is_valid_emf = new_definition
  `is_valid_emf emf
    <=> !r t. corthogonal (h_of_emf emf r t) (e_of_emf emf r t)`;;
let EMF_EQ = 
prove (`!emf:emf r t e h. emf r t = e,h <=> e_of_emf emf r t = e /\ h_of_emf emf r t = h`,
REPEAT STRIP_TAC THEN EQ_TAC THEN SIMP_TAC[e_of_emf;h_of_emf;LET_DEFs;PAIR_EQ;LAMBDA_PAIR] THEN MESON_TAC[PAIR]);;
let emf_at_point_mul = new_definition
  `emf_at_point_mul (f:complex) (emf:complex^3#complex^3) :complex^3#complex^3
    = (f % FST emf, f % SND emf)`;;
overload_interface("%", `emf_at_point_mul :complex->complex^3#complex^3->complex^3#complex^3`);;
let emf_add = new_definition
  `emf_add (emf1:emf) (emf2:emf) :emf = \r t.
    (e_of_emf emf1 r t + e_of_emf emf2 r t,
      h_of_emf emf1 r t + h_of_emf emf2 r t)`;;
overload_interface("+",`emf_add :emf->emf->emf`);; (********************************) (* Plane waves *) (********************************)
let plane_wave = new_definition
  `plane_wave k w e h : emf = \r t. cexp (--ii * Cx(k dot r - w*t)) % (e,h)`;;
let is_plane_wave = new_definition
  `is_plane_wave emf <=>
    is_valid_emf emf /\ ?k w e h.
      &0 < w /\ ~(k = vec 0) /\ emf = plane_wave k w e h
      /\ corthogonal e (vector_to_cvector k)
      /\ corthogonal h (vector_to_cvector k)`;;
let k_of_wave = new_definition
  `k_of_wave emf = @k. ?w e h . &0 < w /\ ~(k = vec 0)
    /\ emf = plane_wave k w e h
    /\ corthogonal e (vector_to_cvector k)
    /\ corthogonal h (vector_to_cvector k)`;;
let w_of_wave = new_definition
  `w_of_wave emf = @w. ?e h . &0 < w /\ emf = plane_wave (k_of_wave emf) w e h
    /\ corthogonal e (vector_to_cvector (k_of_wave emf))
    /\ corthogonal h (vector_to_cvector (k_of_wave emf))`;;
let e_of_wave = new_definition
  `e_of_wave emf = @e. ?h . emf = plane_wave (k_of_wave emf) (w_of_wave emf) e h
    /\ corthogonal e (vector_to_cvector (k_of_wave emf))
    /\ corthogonal h (vector_to_cvector (k_of_wave emf))`;;
let h_of_wave = new_definition
  `h_of_wave emf = @h.
    emf = plane_wave (k_of_wave emf) (w_of_wave emf) (e_of_wave emf) h
    /\ corthogonal (e_of_wave emf) (vector_to_cvector (k_of_wave emf)) 
    /\ corthogonal h (vector_to_cvector (k_of_wave emf))`;;
let eh_of_wave = new_definition
  `eh_of_wave emf = (e_of_wave emf,h_of_wave emf)`;;
let scalar_of_wave = new_definition
  `scalar_of_wave emf =
    \r t. cexp (--ii * Cx(k_of_wave emf dot r - w_of_wave emf * t))`;;
let IS_PLANE_WAVE = 
prove (`!emf. is_plane_wave emf <=> is_valid_emf emf /\ &0 < w_of_wave emf /\ ~(k_of_wave emf = vec 0) /\ emf = plane_wave (k_of_wave emf) (w_of_wave emf) (e_of_wave emf) (h_of_wave emf) /\ corthogonal (e_of_wave emf) (vector_to_cvector (k_of_wave emf)) /\ corthogonal (h_of_wave emf) (vector_to_cvector (k_of_wave emf))`,
let COMMON_TAC x = GEN_TAC THEN EQ_TAC THEN REWRITE_TAC[is_plane_wave;x] THEN MESON_TAC[] in
let lemma1 = prove
    (`!emf. is_plane_wave emf <=> is_valid_emf emf /\ ?w e h . &0 < w
      /\ ~(k_of_wave emf = vec 0) /\ emf = plane_wave (k_of_wave emf) w e h
      /\ corthogonal e (vector_to_cvector (k_of_wave emf))
      /\ corthogonal h (vector_to_cvector (k_of_wave emf))`,
    COMMON_TAC k_of_wave)
  in
  let lemma2 = prove
    (`!emf. is_plane_wave emf <=>
      is_valid_emf emf /\ ?e h . &0 < w_of_wave emf /\ ~(k_of_wave emf = vec 0)
      /\ emf = plane_wave (k_of_wave emf) (w_of_wave emf) e h
      /\ corthogonal e (vector_to_cvector (k_of_wave emf))
      /\ corthogonal h (vector_to_cvector (k_of_wave emf))`,
    REWRITE_TAC[lemma1] THEN COMMON_TAC w_of_wave)
  in
  let lemma3 = prove
    (`!emf. is_plane_wave emf <=>
      is_valid_emf emf /\ ?h . &0 < w_of_wave emf /\ ~(k_of_wave emf = vec 0) /\
      emf = plane_wave (k_of_wave emf) (w_of_wave emf) (e_of_wave emf) h
      /\ corthogonal (e_of_wave emf) (vector_to_cvector (k_of_wave emf))
      /\ corthogonal h (vector_to_cvector (k_of_wave emf))`,
    REWRITE_TAC[lemma2] THEN COMMON_TAC e_of_wave)
  in
  REWRITE_TAC[lemma3] THEN COMMON_TAC h_of_wave);;
let EH_OF_EMF_PLANE_WAVE = 
prove (`!k w e h. (e_of_emf (plane_wave k w e h) = \r t. cexp (--ii * Cx(k dot r - w*t)) % e ) /\ h_of_emf (plane_wave k w e h) = \r t. cexp (--ii * Cx(k dot r - w*t)) % h`,
REWRITE_TAC[e_of_emf;h_of_emf;plane_wave;emf_at_point_mul;LET_DEFs]);;
let non_null_wave = new_definition
  `non_null_wave emf <=>
    is_plane_wave emf /\ ~(e_of_wave emf = cvector_zero)
    /\ ~(h_of_wave emf = cvector_zero)`;;
(***********************************) (* Plane interface between mediums *) (***********************************) (* A medium is characterized by its refractive index *) new_type_abbrev("medium", `:real`);; (* The real^3 vector is a normal to the plane. It is needed in order to fix an * orientation so that we can clearly say on which side is medium 1 or 2 * respectively. *) new_type_abbrev("interface",`:medium # medium # plane # real^3`);;
let FORALL_INTERFACE_THM = 
prove (`!P. (!(i:interface). P i) <=> (!n1 n2 p n. P (n1,n2,p,n))`,
MESON_TAC[PAIR_SURJECTIVE]);;
let is_valid_interface = new_definition
  `is_valid_interface (i:interface) <=>
    let (n1,n2,p,n) = i in
    &0 < n1 /\ &0 < n2 /\ plane p /\ is_normal_to_plane n p`;;
let plane_of_interface = new_definition
  `plane_of_interface (i:interface) = let (n1,n2,p,n) = i in p`;;
let normal_of_interface = new_definition
  `normal_of_interface (i:interface) = let (n1,n2,p,n) = i in n`;;
let n1_of_interface = new_definition
  `n1_of_interface (i:interface) = let (n1,n2,p,n) = i in n1`;;
let n2_of_interface = new_definition
  `n2_of_interface (i:interface) = let (n1,n2,p,n) = i in n2`;;
(* Helpers *)
let IS_VALID_INTERFACE_IS_NORMAL_TO_PLANE = 
prove (`!i. is_valid_interface i ==> is_normal_to_plane (normal_of_interface i) (plane_of_interface i)`,
let NORMAL_OF_INTERFACE_NON_NULL = 
prove (`!i. is_valid_interface i ==> ~(normal_of_interface i = vec 0)`,
let IS_VALID_INTERFACE_PLANE = 
prove (`!i. is_valid_interface i ==> plane (plane_of_interface i)`,
(* [p] is the point where continuity holds; [n] is the normal to the tangent * plane at that point. *)
let boundary_conditions = new_definition
  `boundary_conditions emf1 emf2 p r t <=>
    !n. is_normal_to_plane n p ==>
    let n_ccross = (ccross) (vector_to_cvector n) in
    n_ccross (e_of_emf emf1 r t) = n_ccross (e_of_emf emf2 r t)
    /\ n_ccross (h_of_emf emf1 r t) = n_ccross (h_of_emf emf2 r t)`;;
let eta0 = new_specification ["eta0"] 
  (EXISTS (`?x. &0 < x`,`&1`) (REAL_ARITH `&0 < &1`));;
let k0 = new_specification ["k0"] 
  (EXISTS (`?x. &0 < x`,`&1`) (REAL_ARITH `&0 < &1`));;
let is_plane_wave_at_interface = new_definition
  `is_plane_wave_at_interface i emf_i emf_r emf_t <=>
    is_valid_interface i /\ non_null_wave emf_i /\ is_plane_wave emf_r
    /\ is_plane_wave emf_t
    /\ let (n1,n2,p,n) = i in
      (!pt. pt IN p ==> !t. boundary_conditions (emf_i + emf_r) emf_t p pt t) /\
      (let (k_i,k_r,k_t) = map_triple k_of_wave (emf_i, emf_r, emf_t) in
      (k_i dot n > &0 /\ k_r dot n <= &0 /\ k_t dot n >= &0) /\
      norm k_i = k0 * n1 /\ norm k_r = k0 * n1 /\ norm k_t = k0 * n2) /\
      let emf_in_medium = \emf n.
        h_of_wave emf = Cx(&1 / (eta0 * k0)) %
          (vector_to_cvector (k_of_wave emf) ccross (e_of_wave emf))
      in
      emf_in_medium emf_i n1 /\ emf_in_medium emf_r n1
      /\ emf_in_medium emf_t n2`;;
let IS_PLANE_WAVE_AT_INTERFACE_MAGNITUDES_RELATION = 
prove( `!i emf_i emf_r emf_t. is_plane_wave_at_interface i emf_i emf_r emf_t ==> let (n1,n2,p,n) = i in let (e_i,e_r,e_t) = map_triple (norm o e_of_wave) (emf_i,emf_r,emf_t) in let (h_i,h_r,h_t) = map_triple (norm o h_of_wave) (emf_i,emf_r,emf_t) in h_i = e_i * n1_of_interface i / eta0 /\ h_r = e_r * n1_of_interface i / eta0 /\ h_t = e_t * n2_of_interface i / eta0`,
SIMP_TAC[FORALL_INTERFACE_THM;LET_DEFs;map_triple;o_DEF;n1_of_interface; n2_of_interface;is_plane_wave_at_interface;GSYM CCROSS_LMUL; GSYM VECTOR_TO_CVECTOR_MUL;CCROSS_LREAL;CNORM_NORM_2; REWRITE_RULE[REAL_ARITH `x+y=z <=> x=z-y:real`] NORM_CROSS_DOT;non_null_wave; REAL_ARITH `(x-y)+(z-t)=(x+z)-(y+t):real`;REAL_POW_MUL; GSYM REAL_ADD_LDISTRIB;IS_PLANE_WAVE;CORTHOGONAL_REAL_CLAUSES;DOT_LMUL] THEN ONCE_REWRITE_TAC[ORTHOGONAL_SYM] THEN SIMP_TAC[orthogonal;REAL_POW_2; REAL_MUL_RZERO;REAL_ADD_RID;REAL_SUB_RZERO] THEN SIMP_TAC[GSYM REAL_POW_2;REAL_LE_POW_2;MESON[REAL_LE_ADD;REAL_LE_POW_2] `!x y. &0 <= x pow 2 + y pow 2`;SQRT_MUL;CX_MUL] THEN REWRITE_TAC[GSYM CNORM_NORM_2;VECTOR_TO_CVECTOR_CVECTOR_RE_IM] (* THEN SIMP_TAC[NORM_POS_LE;POW_2_SQRT;NORM_MUL;real_div;REAL_MUL_LID] *) THEN SIMP_TAC[NORM_POS_LE;SQRT_POW_2;NORM_MUL;real_div;REAL_MUL_LID] THEN REWRITE_TAC[MESON[REAL_LT_IMP_LE;REAL_ABS_REFL;REAL_LE_INV;REAL_LE_MUL; eta0;k0] `abs (inv (eta0 * k0)) = inv (eta0 * k0)`] THEN REWRITE_TAC[REAL_INV_MUL;REAL_ARITH `(x*y)*z*t=x*(y*z)*t:real`;MATCH_MP REAL_MUL_LINV (GSYM (MATCH_MP REAL_LT_IMP_NE k0));REAL_MUL_LID] THEN REWRITE_TAC[COMPLEX_MUL_SYM;REAL_MUL_SYM]);;
(* Helpers *)
let IS_PLANE_WAVE_AT_INTERFACE_BOUNDARY_CONDITIONS = 
prove (`!i emf_i emf_r emf_t. is_plane_wave_at_interface i emf_i emf_r emf_t ==> !pt. pt IN (plane_of_interface i) ==> !t. boundary_conditions (emf_i + emf_r) emf_t (plane_of_interface i) pt t`,
let IS_PLANE_WAVE_AT_INTERFACE_IS_NORMAL_TO_PLANE = 
prove (`!i emf_i emf_r emf_t. is_plane_wave_at_interface i emf_i emf_r emf_t ==> is_normal_to_plane (normal_of_interface i) (plane_of_interface i)`,
let IS_PLANE_WAVE_AT_INTERFACE_PLANE = 
prove (`!i emf_i emf_r emf_t. is_plane_wave_at_interface i emf_i emf_r emf_t ==> plane (plane_of_interface i)`,
let IS_PLANE_WAVE_AT_INTERFACE_NON_NULL_NORMAL = 
prove (`!i emf_i emf_r emf_t. is_plane_wave_at_interface i emf_i emf_r emf_t ==> ~(normal_of_interface i = vec 0)`,
(***********************************) (* TE and TM modes *) (***********************************)
let is_mode_axis = new_definition
  `is_mode_axis field (i:interface) emf_i emf_r emf_t v <=>
    orthogonal v (normal_of_interface i) /\ norm v = &1 /\ !r t.
    collinear_cvectors (field emf_i r t) (vector_to_cvector v)
    /\ collinear_cvectors (field emf_r r t) (vector_to_cvector v)
    /\ collinear_cvectors (field emf_t r t) (vector_to_cvector v)`;;
let transverse_mode = new_definition
  `transverse_mode field (i:interface) emf_i emf_r emf_t <=>
    ?v. is_mode_axis field (i:interface) emf_i emf_r emf_t v`;;
let TE_mode = new_definition `TE_mode = transverse_mode e_of_emf`;;
let TM_mode = new_definition `TM_mode = transverse_mode h_of_emf`;;
let mode_axis = new_definition `mode_axis field i emf_i emf_r emf_t =
  @v. is_mode_axis field i emf_i emf_r emf_t v`;;
let TE_mode_axis = new_definition `TE_mode_axis = mode_axis e_of_emf`;;
let TM_mode_axis = new_definition `TM_mode_axis = mode_axis h_of_emf`;;
let TRANSVERSE_MODE_THM = 
prove (`!field i emf_i emf_r emf_t. transverse_mode field i emf_i emf_r emf_t <=> is_mode_axis field i emf_i emf_r emf_t (mode_axis field i emf_i emf_r emf_t)`,
REWRITE_TAC[transverse_mode;mode_axis] THEN SELECT_ELIM_TAC THEN REWRITE_TAC[]);;
let TE_MODE_THM = 
prove (`!i emf_i emf_r emf_t. TE_mode i emf_i emf_r emf_t <=> is_mode_axis e_of_emf i emf_i emf_r emf_t (TE_mode_axis i emf_i emf_r emf_t)`,
let TM_MODE_THM = 
prove (`!i emf_i emf_r emf_t. TM_mode i emf_i emf_r emf_t <=> is_mode_axis h_of_emf i emf_i emf_r emf_t (TM_mode_axis i emf_i emf_r emf_t)`,
let MODE_AXIS_ORTHOGONAL_N = 
prove (`!field i emf_i emf_r emf_t. transverse_mode field i emf_i emf_r emf_t ==> orthogonal (mode_axis field i emf_i emf_r emf_t) (normal_of_interface i)`,
let TE_MODE_AXIS_ORTHOGONAL_N = 
prove (`!i emf_i emf_r emf_t. TE_mode i emf_i emf_r emf_t ==> orthogonal (TE_mode_axis i emf_i emf_r emf_t) (normal_of_interface i)`,
let TM_MODE_AXIS_ORTHOGONAL_N = 
prove (`!i emf_i emf_r emf_t. TM_mode i emf_i emf_r emf_t ==> orthogonal (TM_mode_axis i emf_i emf_r emf_t) (normal_of_interface i)`,
let TRANSVERSE_MODE_PROJ = 
prove (`!field i emf_i emf_r emf_t. transverse_mode field i emf_i emf_r emf_t ==> !r t. let x = vector_to_cvector (mode_axis field i emf_i emf_r emf_t) in field emf_i r t = (field emf_i r t cdot x) % x /\ field emf_r r t = (field emf_r r t cdot x) % x /\ field emf_t r t = (field emf_t r t cdot x) % x`,
REWRITE_TAC[TRANSVERSE_MODE_THM;is_mode_axis;LET_DEFs] THEN REPEAT STRIP_TAC THEN POP_ASSUM (STRIP_ASSUME_TAC o SPEC_ALL) THEN FIRST_ASSUM (ASSUME_TAC o REWRITE_RULE[NORM_EQ_0;GSYM VECTOR_TO_CVECTOR_ZERO_EQ] o MATCH_MP (REAL_ARITH `!x. x= &1 ==> ~(x= &0)`)) THEN ASM_SIMP_TAC[COLLINEAR_RUNITREAL]);;
let TE_MODE_PROJ = 
prove (`!i emf_i emf_r emf_t. TE_mode i emf_i emf_r emf_t ==> !r t. let x = vector_to_cvector (TE_mode_axis i emf_i emf_r emf_t) in e_of_emf emf_i r t = (e_of_emf emf_i r t cdot x) % x /\ e_of_emf emf_r r t = (e_of_emf emf_r r t cdot x) % x /\ e_of_emf emf_t r t = (e_of_emf emf_t r t cdot x) % x`,
let TM_MODE_PROJ = 
prove (`!i emf_i emf_r emf_t. TM_mode i emf_i emf_r emf_t ==> !r t. let x = vector_to_cvector (TM_mode_axis i emf_i emf_r emf_t) in h_of_emf emf_i r t = (h_of_emf emf_i r t cdot x) % x /\ h_of_emf emf_r r t = (h_of_emf emf_r r t cdot x) % x /\ h_of_emf emf_t r t = (h_of_emf emf_t r t cdot x) % x`,
let TE_MODE_PLANEWAVE_PROJ = 
prove (`!i emf_i emf_r emf_t. TE_mode i emf_i emf_r emf_t /\ is_plane_wave emf_i /\ is_plane_wave emf_r /\ is_plane_wave emf_t ==> let x = vector_to_cvector (TE_mode_axis i emf_i emf_r emf_t) in e_of_wave emf_i = (e_of_wave emf_i cdot x) % x /\ e_of_wave emf_r = (e_of_wave emf_r cdot x) % x /\ e_of_wave emf_t = (e_of_wave emf_t cdot x) % x`,
REWRITE_TAC[IS_PLANE_WAVE] THEN REPEAT STRIP_TAC THEN ASM_CSQ_THEN TE_MODE_PROJ (MP_TAC o REWRITE_RULE[LET_DEFs]) THEN ASSUM_LIST (REWRITE_TAC o mapfilter (REWRITE_RULE[EH_OF_EMF_PLANE_WAVE] o MATCH_MP (MESON[] `x=y ==> e_of_emf x = e_of_emf y`))) THEN REWRITE_TAC[CDOT_LMUL;GSYM CVECTOR_MUL_ASSOC;CVECTOR_MUL_LCANCEL; CEXP_NZ;LET_DEFs]);;
let TM_MODE_PLANEWAVE_PROJ = 
prove (`!i emf_i emf_r emf_t. TM_mode i emf_i emf_r emf_t /\ is_plane_wave emf_i /\ is_plane_wave emf_r /\ is_plane_wave emf_t ==> let x = vector_to_cvector (TM_mode_axis i emf_i emf_r emf_t) in h_of_wave emf_i = (h_of_wave emf_i cdot x) % x /\ h_of_wave emf_r = (h_of_wave emf_r cdot x) % x /\ h_of_wave emf_t = (h_of_wave emf_t cdot x) % x`,
REWRITE_TAC[IS_PLANE_WAVE] THEN REPEAT STRIP_TAC THEN ASM_CSQ_THEN TM_MODE_PROJ (MP_TAC o REWRITE_RULE[LET_DEFs]) THEN ASSUM_LIST (REWRITE_TAC o mapfilter (REWRITE_RULE[EH_OF_EMF_PLANE_WAVE] o MATCH_MP (MESON[] `x=y ==> h_of_emf x = h_of_emf y`))) THEN REWRITE_TAC[CDOT_LMUL;GSYM CVECTOR_MUL_ASSOC;CVECTOR_MUL_LCANCEL;CEXP_NZ; LET_DEFs]);;