(* 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)`,
  SIMP_TAC[FORALL_INTERFACE_THM;is_valid_interface;LET_DEFs;
    normal_of_interface;plane_of_interface]);;

let NORMAL_OF_INTERFACE_NON_NULL = prove
  (`!i. is_valid_interface i ==> ~(normal_of_interface i = vec 0)`,
  SIMP_TAC[FORALL_INTERFACE_THM;is_valid_interface;is_normal_to_plane;
    normal_of_interface;LET_DEFs]);;
let IS_VALID_INTERFACE_PLANE = prove
  (`!i. is_valid_interface i ==> plane (plane_of_interface i)`,
  SIMP_TAC[FORALL_INTERFACE_THM;is_valid_interface;plane_of_interface;
    LET_DEFs]);;

(* [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`,
  SIMP_TAC[FORALL_INTERFACE_THM;is_plane_wave_at_interface;plane_of_interface;
    LET_DEFs]);;

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)`,
  SIMP_TAC[is_plane_wave_at_interface;IS_VALID_INTERFACE_IS_NORMAL_TO_PLANE]);;

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)`,
  SIMP_TAC[is_plane_wave_at_interface;IS_VALID_INTERFACE_PLANE]);;

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)`,
  SIMP_TAC[FORALL_INTERFACE_THM;is_plane_wave_at_interface;is_valid_interface;
    normal_of_interface;is_normal_to_plane;LET_DEFs]);;


(***********************************)
(* 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)`,
  REWRITE_TAC[TE_mode;TE_mode_axis;TRANSVERSE_MODE_THM]);;

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)`,
  REWRITE_TAC[TM_mode;TM_mode_axis;TRANSVERSE_MODE_THM]);;

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)`,
  SIMP_TAC[TRANSVERSE_MODE_THM;is_mode_axis]);;

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)`,
  REWRITE_TAC[TE_mode;TE_mode_axis;MODE_AXIS_ORTHOGONAL_N]);;

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)`,
  REWRITE_TAC[TM_mode;TM_mode_axis;MODE_AXIS_ORTHOGONAL_N]);;

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`,
  REWRITE_TAC[TE_mode;TE_mode_axis;TRANSVERSE_MODE_PROJ]);;

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`,
  REWRITE_TAC[TM_mode;TM_mode_axis;TRANSVERSE_MODE_PROJ]);;

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]);;