(* ========================================================================= *)
(*               Formalization of Electromagnetic Optics                     *)
(*                                                                           *)
(*            (c) Copyright, Sanaz Khan Afshar & Vincent Aravantinos 2011-14 *)
(*                           Hardware Verification Group,                    *)
(*                           Concordia University                            *)
(*                                                                           *)
(*                Contact:   <s_khanaf@encs.concordia.ca>                    *)
(*                           <vincent@encs.concordia.ca>                     *)
(*                                                                           *)
(* This file deals with the non-trivial theorems called "primitive rules".   *)
(* ========================================================================= *)


let WAVE_ORTHOGONALITY = 
prove
  (`!emf. is_plane_wave emf ==> corthogonal (h_of_wave emf) (e_of_wave emf)`,

  REWRITE_TAC[IS_PLANE_WAVE;is_valid_emf;plane_wave;emf_at_point_mul;e_of_emf;
    h_of_emf;LET_DEF;LET_END_DEF]
  THEN REPEAT STRIP_TAC THEN REPLICATE_TAC 2 (POP_ASSUM (K ALL_TAC))
  THEN POP_ASSUM (fun x -> RULE_ASSUM_TAC (ONCE_REWRITE_RULE[x]))
  THEN RULE_ASSUM_TAC (REWRITE_RULE[BETA_THM;LET_DEF;LET_END_DEF])
  THEN REPEAT (POP_ASSUM MP_TAC)
  THEN SIMP_TAC[corthogonal;CDOT3;CVECTOR_MUL_COMPONENT;CNJ_MUL;
    SIMPLE_COMPLEX_ARITH `(x*y)*cnj z*t=(x*cnj z)*(y*t):complex`;
    GSYM COMPLEX_ADD_LDISTRIB;COMPLEX_ENTIRE;CEXP_NZ;CNJ_EQ_0]);;

let NON_NULL_LEMMA = 
prove
  (`!emf v. ~(v = vec 0) /\ non_null_wave emf ==>
    let v_ccross = (ccross) (vector_to_cvector v) in
    ~(v_ccross (e_of_wave emf) = cvector_zero)
    \/ ~(v_ccross (h_of_wave emf) = cvector_zero)`,

  REPEAT GEN_TAC THEN CONV_TAC (DEPTH_CONV let_CONV)
  THEN ASM_CASES_TAC `vector_to_cvector v ccross e_of_wave emf = cvector_zero`
  THEN ASM_REWRITE_TAC[non_null_wave]
  THEN RULE_ASSUM_TAC (REWRITE_RULE[CCROSS_COLLINEAR_CVECTORS])
  THEN ASSUME_CONSEQUENCES (REWRITE_RULE[IMP_CONJ]
    CORTHOGONAL_COLLINEAR_CVECTORS)
  THEN REPEAT STRIP_TAC
  THEN SUBGOAL_THEN `~(vector_to_cvector (v:real^3) = cvector_zero)` ASSUME_TAC
  THENL ON_FIRST_GOAL 
    (REWRITE_TAC[CART_EQ3;CVECTOR_ZERO_COMPONENT;VECTOR_TO_CVECTOR_COMPONENT;
      CX_INJ]
    THEN RULE_ASSUM_TAC (REWRITE_RULE[CART_EQ;VEC_COMPONENT;DIMINDEX_3;FORALL_3])
    THEN ASM_REWRITE_TAC[])
  THEN SUBGOAL_THEN
    `?a . ~(a=Cx(&0)) /\ vector_to_cvector v = a % e_of_wave emf` STRIP_ASSUME_TAC
  THENL ON_FIRST_GOAL (ASM_MESON_TAC[NON_NULL_COLLINEARS;COLLINEAR_CVECTORS_SYM])
  THEN ASSUME_CONSEQUENCES WAVE_ORTHOGONALITY
  THEN POP_ASSUM MP_TAC
  THEN POP_ASSUM (fun x ->
    RULE_ASSUM_TAC (REWRITE_RULE[x;CCROSS_LMUL;CVECTOR_MUL_EQ_0]))
  THEN POP_ASSUM (fun x ->
    RULE_ASSUM_TAC (REWRITE_RULE[x;CCROSS_COLLINEAR_CVECTORS]))
  THEN ASM_MESON_TAC[CORTHOGONAL_COLLINEAR_CVECTORS;CORTHOGONAL_SYM]);;

let NON_NULL_LEMMA_PASTECART = 
prove
  (`!emf v. ~(v = vec 0) /\ non_null_wave emf ==>
    let v_ccross = (ccross) (vector_to_cvector v) in
    ~(pastecart (v_ccross (e_of_wave emf)) (v_ccross (h_of_wave emf)) =
      cvector_zero)`,

  REWRITE_TAC[PASTECART_EQ_CVECTOR_ZERO;DE_MORGAN_THM;NON_NULL_LEMMA]);;

let BOUNDARY_CONDITIONS_FOR_PLANE_WAVES = 
prove
  (`!i emf_i emf_r emf_t.
    is_plane_wave_at_interface i emf_i emf_r emf_t ==>
    let p = plane_of_interface i in
    !n. is_normal_to_plane n p ==>
    let n_ccross = (ccross) (vector_to_cvector n) in
    !pt. pt IN p ==> !t.
      let plane_component = \f_of_wave emf. cexp (--ii * Cx((k_of_wave emf) dot
        pt - w_of_wave emf*t)) % n_ccross (f_of_wave emf) in
      plane_component e_of_wave emf_i + plane_component e_of_wave emf_r
        = plane_component e_of_wave emf_t
      /\ plane_component h_of_wave emf_i + plane_component h_of_wave emf_r
        = plane_component h_of_wave emf_t`,

  REWRITE_TAC[FORALL_INTERFACE_THM;plane_of_interface;LET_DEF;LET_END_DEF;
    is_plane_wave_at_interface;non_null_wave;IS_PLANE_WAVE;plane_wave;
    emf_at_point_mul;e_of_emf;h_of_emf;LET_DEF;LET_END_DEF;map_triple;o_DEF]
  THEN REPEAT STRIP_TAC
  THEN POP_ASSUM (fun x -> FIRST_X_ASSUM (fun y -> MP_TAC (MATCH_MP y x)))
  THEN ASM (GEN_REWRITE_TAC (RATOR_CONV o DEPTH_CONV)) []
  THEN REWRITE_TAC[boundary_conditions;emf_add;e_of_emf;h_of_emf;LET_DEF;
    LET_END_DEF;CCROSS_RADD;CCROSS_RMUL]
  THEN DISCH_THEN (C ASM_CSQ_THEN STRIP_ASSUME_TAC) THEN ASM_REWRITE_TAC[]);;

(* We combine the E and H field as one complex^6 vector. 
 * Convenient for some proofs. *)
let PASTECART_BOUNDARY_CONDITIONS_FOR_PLANE_WAVES = 
prove
  (`!i emf_i emf_r emf_t. is_plane_wave_at_interface i emf_i emf_r emf_t ==>
  let p = plane_of_interface i in
  !n. is_normal_to_plane n p ==>
  let n_ccross = (ccross) (vector_to_cvector n) in
	!pt. pt IN p ==> !t.
    let plane_component = \emf.
      cexp (--ii * Cx((k_of_wave emf) dot pt - w_of_wave emf*t)) % pastecart
        (n_ccross (e_of_wave emf)) (n_ccross (h_of_wave emf)) in
    plane_component emf_i + plane_component emf_r = plane_component emf_t`,

    REWRITE_TAC[LET_DEF;LET_END_DEF] THEN REPEAT STRIP_TAC
    THEN ASM_CSQ_THEN (REWRITE_RULE[LET_DEF;LET_END_DEF]
      BOUNDARY_CONDITIONS_FOR_PLANE_WAVES) (C ASM_CSQ_THEN (C ASM_CSQ_THEN
        (MP_TAC o SPEC_ALL)))
    THEN PASTECART_TAC[GSYM PASTECART_CVECTOR_MUL;PASTECART_CVECTOR_ADD]);;

let EXISTS_NORMAL_OF_PLANE_INTERFACE = 
prove
  (`!i emf_i emf_r emf_t. is_plane_wave_at_interface i emf_i emf_r emf_t ==>
    ?n. is_normal_to_plane n (plane_of_interface i)`,

  REWRITE_TAC[FORALL_INTERFACE_THM;LET_DEF;LET_END_DEF;
    is_plane_wave_at_interface;is_valid_interface;plane_of_interface]
  THEN MESON_TAC[EXISTS_NORMAL_OF_PLANE]);;

let FREQUENCY_CONSERVATION = 
prove
  (`!i emf_i emf_r emf_t. is_plane_wave_at_interface i emf_i emf_r emf_t ==>
  (non_null_wave emf_r ==> w_of_wave emf_r = w_of_wave emf_i)
  /\ (non_null_wave emf_t ==> w_of_wave emf_t = w_of_wave emf_i)`,

  REWRITE_TAC[FORALL_INTERFACE_THM] THEN REPEAT STRIP_TAC
  THEN ASM_CSQ_THEN EXISTS_NORMAL_OF_PLANE_INTERFACE STRIP_ASSUME_TAC
  THEN SUBGOAL_THEN `~(n' = vec 0:real^3) /\ plane (p:plane)`
    STRIP_ASSUME_TAC THEN
  TRY (RULE_ASSUM_TAC (REWRITE_RULE[is_plane_wave_at_interface;LET_DEF;
    LET_END_DEF;is_valid_interface;is_normal_to_plane])
    THEN ASM_REWRITE_TAC[] THEN NO_TAC)
  THEN ASM_CSQ_THEN PLANE_NON_EMPTY (STRIP_ASSUME_TAC o REWRITE_RULE[GSYM
    MEMBER_NOT_EMPTY])
  THEN ASM_CSQ_THEN (REWRITE_RULE[LET_DEF;LET_END_DEF]
    PASTECART_BOUNDARY_CONDITIONS_FOR_PLANE_WAVES) (C ASM_CSQ_THEN ASSUME_TAC)
  THEN RULE_ASSUM_TAC (REWRITE_RULE[plane_of_interface;LET_DEF;LET_END_DEF])
  THEN POP_ASSUM (C ASM_CSQ_THEN MP_TAC)
  THEN REWRITE_TAC[CX_SUB;CX_MUL;SIMPLE_COMPLEX_ARITH
    `--ii*(x-y*z) = (ii*y)*z+ --ii*x`;CEXP_ADD;GSYM CVECTOR_MUL_ASSOC]
  THEN ONCE_REWRITE_TAC[MESON[COMPLEX_EQ_MUL_LCANCEL;II_NZ;CX_INJ]
    `x=y <=> ii * Cx x = ii * Cx y`]
  THENL
    (let APPLY_FREQ_EQ_TAC x =
      DISCH_THEN (MP_TAC o MATCH_MP (REWRITE_RULE[MESON[]
        `(A /\ B /\ C ==> D) <=> (C ==> A ==> B ==> D)`] x)) in
    map APPLY_FREQ_EQ_TAC [ VEC_WAVE_SUM_EQ_WEAK2; VEC_WAVE_SUM_EQ_WEAK1])
  THEN ANTS_TAC
  THEN
  let APPLY_NON_NULL_LEMMA = REWRITE_TAC[CVECTOR_MUL_EQ_0;CEXP_NZ]
    THEN MATCH_MP_TAC (CONV_RULE (DEPTH_CONV let_CONV) NON_NULL_LEMMA_PASTECART) in
  REPEAT (ANTS_TAC ORELSE APPLY_NON_NULL_LEMMA)
  THEN ASM_MESON_TAC[is_plane_wave_at_interface]);;

let IS_PLANE_WAVE_AT_INTERFACE_THMS =
  [is_plane_wave_at_interface;LET_DEF;LET_END_DEF;map_triple;o_DEF;
    is_valid_interface];;

let K_PLANE_PROJECTION_CONSERVATION = 
prove
  (`!i emf_i emf_r emf_t. is_plane_wave_at_interface i emf_i emf_r emf_t ==>
    let n = unit (normal_of_interface i) in
    (non_null_wave emf_r ==>
      projection_on_hyperplane (k_of_wave emf_r) n =
        projection_on_hyperplane (k_of_wave emf_i) n)
    /\ (non_null_wave emf_t ==>
      projection_on_hyperplane (k_of_wave emf_t) n =
        projection_on_hyperplane (k_of_wave emf_i) n)`,

    REWRITE_TAC[FORALL_INTERFACE_THM] THEN REPEAT GEN_TAC
    THEN DISCH_THEN (DISTRIB [MP_TAC; LABEL_TAC "*" o MATCH_MP
      PASTECART_BOUNDARY_CONDITIONS_FOR_PLANE_WAVES])
    THEN REWRITE_TAC([normal_of_interface] @ IS_PLANE_WAVE_AT_INTERFACE_THMS)
    THEN REPEAT STRIP_TAC
    THEN REMOVE_THEN "*" (C ASM_CSQ_THEN ASSUME_TAC o
      REWRITE_RULE[LET_DEF;LET_END_DEF;plane_of_interface])
    THEN ASM_CSQ_THEN FORALL_PLANE_THM_2 STRIP_ASSUME_TAC
    THEN POP_ASSUM (RULE_ASSUM_TAC o SINGLE REWRITE_RULE)
    THEN FIRST_ASSUM (STRIP_ASSUME_TAC o CONV_RULE
      (REWR_CONV is_normal_to_plane))
    THEN ASM_SIMP_TAC[UNIT_THM;PROJECTION_ON_HYPERPLANE_THM]
    THEN ASM_CSQ_THEN PLANE_DECOMP_DOT (C ASM_CSQ_THEN (C
      ASM_CSQ_THEN ASSUME_TAC))
    THEN MAP_EVERY (fun x -> REWRITE_TAC[VECTOR_ARITH x]
      THEN ASSUM_LIST (GEN_REWRITE_TAC (RATOR_CONV o ONCE_DEPTH_CONV)))
      [`a-b=c-d <=> a=c+b-d:real^N`;`a=c+b-d <=> c=a-b+d:real^N`]
    THEN REWRITE_TAC[VECTOR_ARITH `a+b+c=(d+e+f)-f+c <=> a+b=d+e:real^N`]
    THEN ASM_CSQ_THEN BASIS_NON_NULL (fun x -> FIRST_ASSUM
      (STRIP_ASSUME_TAC o REWRITE_RULE[IN_INSERT;IN_SING;
      MESON[] `(!x. x=a \/ x=b ==> p x) <=> p a /\ p b`;UNIT_EQ_ZERO]
      o MATCH_MP x o CONJUNCT1 o CONV_RULE
        (REWR_CONV is_orthogonal_plane_basis)))
    THEN ASM_SIMP_TAC[UNIT_UNIT] THEN MK_COMB_TAC
    THENL [AP_TERM_TAC;ALL_TAC;AP_TERM_TAC;ALL_TAC]
    THEN AP_THM_TAC THEN AP_TERM_TAC 
    THEN ONCE_REWRITE_TAC[MESON[COMPLEX_EQ_MUL_LCANCEL;II_NZ;CX_INJ;
      COMPLEX_NEG_EQ_0] `x=y <=> --ii * Cx x = --ii * Cx y`]
    THEN RULE_ASSUM_TAC (REWRITE_RULE[
      VECTOR_ARITH `x dot (y+a%z+b%t) = x dot y+(x dot z)*a+(x dot t)*b`;
      SIMPLE_COMPLEX_ARITH `--ii * Cx ((a+b*b'+c*c')-d) = (--ii * Cx b)
        * Cx b' + (--ii * Cx c) * Cx c' + --ii * Cx a + ii * Cx d`;
      CEXP_ADD])
    THENL (map (fun f -> FIRST_X_ASSUM (ASSUME_TAC o f)) [
      funpow 2 (SPEC `&0` o ONCE_REWRITE_RULE[SWAP_FORALL_THM]);
      SPEC `&0` o ONCE_REWRITE_RULE[SWAP_FORALL_THM] o SPEC `&0`;
      funpow 2 (SPEC `&0` o ONCE_REWRITE_RULE[SWAP_FORALL_THM]);
      SPEC `&0` o ONCE_REWRITE_RULE[SWAP_FORALL_THM] o SPEC `&0`;
    ])
    THEN RULE_ASSUM_TAC (REWRITE_RULE[COMPLEX_MUL_RZERO;CEXP_0;COMPLEX_MUL_RID;
      REAL_MUL_RZERO;GSYM CVECTOR_MUL_ASSOC;CVECTOR_MUL_ID])
    THENL
    (let lemma = MESON[] `(!x y z t u v. p x y z t v u ==> q x y z /\ r t v u)
      ==> (!x y z t u v. p x y z t v u ==> q x y z)` in
    let APPLY_FREQ_EQ_TAC x =
      POP_ASSUM (DISTRIB [ASSUME_TAC;MP_TAC o MATCH_MP (REWRITE_RULE[MESON[]
        `(A /\ B /\ C ==> D) <=> (C ==> A ==> B ==> D)`](MATCH_MP lemma x))])
    in
    map APPLY_FREQ_EQ_TAC [VEC_WAVE_SUM_EQ_WEAK2;VEC_WAVE_SUM_EQ_WEAK2;
      VEC_WAVE_SUM_EQ_WEAK1;VEC_WAVE_SUM_EQ_WEAK1])
    THEN ANTS_TAC
    THEN
      let APPLY_NON_NULL_LEMMA =
        REWRITE_TAC[CVECTOR_MUL_EQ_0;CEXP_NZ]
        THEN MATCH_MP_TAC (CONV_RULE (DEPTH_CONV let_CONV)
          NON_NULL_LEMMA_PASTECART)
      in
      REPEAT (ANTS_TAC ORELSE APPLY_NON_NULL_LEMMA)
    THEN ASM_REWRITE_TAC[]);;

let LAW_OF_REFLECTION = 
prove
  (`!i emf_i emf_r emf_t. is_plane_wave_at_interface i emf_i emf_r emf_t ==>
    let n = unit (normal_of_interface i) in
    non_null_wave emf_r ==>
      symetric_vectors_wrt (--(k_of_wave emf_r)) (k_of_wave emf_i) n`,

    REWRITE_TAC[FORALL_INTERFACE_THM;LET_DEFs] THEN REPEAT STRIP_TAC
    THEN ASM_CSQ_THEN (REWRITE_RULE[LET_DEFs] K_PLANE_PROJECTION_CONSERVATION)
      (C ASM_CSQ_THEN MP_TAC o CONJUNCT1)
    THEN ASM_CSQ_THEN (MESON[is_plane_wave_at_interface]
      `is_plane_wave_at_interface i emf_i emf_r emf_t ==> is_valid_interface i`)
      (DISTRIB [ASSUME_TAC; STRIP_ASSUME_TAC o
        REWRITE_RULE[LET_DEFs;is_valid_interface]])
    THEN ASM_SIMP_TAC[normal_of_interface;LET_DEFs;
      PROJECTION_ON_HYPERPLANE_THM;symetric_vectors_wrt;
      NORMAL_OF_INTERFACE_NON_NULL;UNIT_THM;UNIT_EQ_ZERO]
    THEN DISCH_TAC
    THEN SUBGOAL_THEN `!x:real^3. orthogonal (x - (x dot unit n) % unit n) ((x
      dot unit n) % unit n)` ASSUME_TAC
    THENL ON_FIRST_GOAL (ASM_MESON_TAC [ORTHOGONAL_RUNIT;ORTHOGONAL_CLAUSES;
        SUB_UNIT_NORMAL_IS_ORTHOGONAL_TO_NORMAL])
    THEN FIRST_ASSUM (REPEAT_TCL STRIP_THM_THEN (fun x ->
      if contains_sub_term_name "k_of_wave" (concl x) then ASSUME_TAC x
      else ALL_TAC)
      o REWRITE_RULE[LET_DEFs;map_triple;o_DEF] o CONV_RULE (REWR_CONV
        is_plane_wave_at_interface))
    THEN REPLICATE_TAC 3 (POP_ASSUM (K ALL_TAC))
    THEN POP_ASSUM (fun _ -> POP_ASSUM (fun x -> POP_ASSUM (fun y -> MP_TAC
      (REWRITE_RULE[NORM_EQ] (TRANS x (GSYM y))))))
    THEN FIRST_ASSUM (ASSUME_TAC o SPEC `k_of_wave emf`)
    THEN ASM_CSQ_THEN ORTHOGONAL_SQR_NORM (SINGLE REWRITE_TAC)
    THEN ASM_REWRITE_TAC[REAL_EQ_ADD_LCANCEL;DOT_LMUL;DOT_RMUL;DOT_SQUARE_NORM]
    THEN SUBGOAL_TAC "" `~(n = vec 0 :real^3)` 
      [ASM_MESON_TAC[is_normal_to_plane]]
    THEN ASM_SIMP_TAC[UNIT_THM;REAL_POW_2;REAL_MUL_RID]
    THEN REWRITE_TAC[GSYM REAL_POW_2;GSYM REAL_EQ_SQUARE_ABS;real_abs]
    THEN REPEAT COND_CASES_TAC THEN RULE_ASSUM_TAC (REWRITE_RULE 
      [real_ge;REWRITE_RULE[real_lt;MESON[] `(~p <=> ~q) <=> (p <=> q)`]
      DOT_RUNIT_LT0])
    THENL [
      FIRST_X_ASSUM (fun x -> FIRST_X_ASSUM (ASSUME_TAC 
        o ONCE_REWRITE_RULE[GSYM DOT_RUNIT_EQ0] o GSYM
        o CONV_RULE (REWR_CONV REAL_LE_ANTISYM) o CONJ x))
      THEN ASM_REWRITE_TAC[DOT_LNEG;REAL_NEG_0;REAL_MUL_RZERO;VECTOR_MUL_LZERO;
        VECTOR_ARITH `--x+y=vec 0 <=> x=y`]
      THEN DISCH_THEN (RULE_ASSUM_TAC o SINGLE REWRITE_RULE o GSYM)
      THEN POP_ASSUM (RULE_ASSUM_TAC o SINGLE REWRITE_RULE)
      THEN RULE_ASSUM_TAC (REWRITE_RULE[VECTOR_MUL_LZERO;VECTOR_SUB_RZERO])
      THEN ASM_REWRITE_TAC[];
      ASM_MESON_TAC[REAL_ARITH `x > &0 ==> &0 <= x`];
      DISCH_THEN (ASSUME_TAC o GSYM)
      THEN FIRST_X_ASSUM (fun x -> ignore (term_match [] `x-y=z-t` (concl x));
        MP_TAC x)
      THEN POP_ASSUM (fun x -> REWRITE_TAC(x::[VECTOR_MUL_LNEG;DOT_LNEG;
        REAL_MUL_RNEG;VECTOR_ARITH `x-y=z- --y <=> --x+z = --(&2%y)`;
        VECTOR_MUL_ASSOC]));
      ASM_MESON_TAC[REAL_ARITH `x > &0 ==> &0 <= x`];
    ]);;

let PLANE_OF_INCIDENCE_LAW = 
prove
  (`!i emf_i emf_r emf_t.
    is_plane_wave_at_interface i emf_i emf_r emf_t /\ non_null_wave emf_r
    /\ non_null_wave emf_t ==>
    coplanar {vec 0, k_of_wave emf_i, k_of_wave emf_r, k_of_wave emf_t,
      normal_of_interface i}`,

  REWRITE_TAC[FORALL_INTERFACE_THM;LET_DEFs;normal_of_interface]
  THEN REPEAT STRIP_TAC
  THEN REWRITE_TAC[coplanar] THEN REPEAT STRIP_TAC
  THEN MAP_EVERY EXISTS_TAC [`vec 0:real^3`;`k_of_wave emf_i`;`unit n:real^3`]
  THEN ASM_CSQ_THEN (MESON[is_plane_wave_at_interface]
    `is_plane_wave_at_interface i emf_i emf_r emf_t ==> is_valid_interface i`)
    (DISTRIB [ASSUME_TAC; STRIP_ASSUME_TAC o
      REWRITE_RULE[LET_DEFs;is_valid_interface]])
  THEN ASM_CSQ_THEN (REWRITE_RULE[LET_DEFs] K_PLANE_PROJECTION_CONSERVATION)
    (DISTRIB (map ((o) (C ASM_CSQ_THEN ASSUME_TAC)) [CONJUNCT1;CONJUNCT2]))
  THEN REWRITE_TAC[INSERT_SUBSET;SING_SUBSET] THEN REPEAT CONJ_TAC
  THENL 
  let IN_SET_TAC = MATCH_MP_TAC HULL_INC THEN SET_TAC[] in
  let COMBINATION_TAC =
    POP_ASSUM MP_TAC THEN POP_ASSUM MP_TAC 
    THEN ASM_SIMP_TAC[AFFINE_HULL_3_ZERO;IN_ELIM_THM;UNIV;normal_of_interface;
      LET_DEFs;PROJECTION_ON_HYPERPLANE_THM;symetric_vectors_wrt;
      NORMAL_OF_INTERFACE_NON_NULL;UNIT_THM;UNIT_EQ_ZERO;
      VECTOR_ARITH `x-a%y=z-b%y <=> x= &1%z+(a-b)%y`]
    THEN MESON_TAC[]
  in
  [ IN_SET_TAC; IN_SET_TAC; COMBINATION_TAC; COMBINATION_TAC;
    REWRITE_TAC[AFFINE_HULL_3_ZERO;IN_ELIM_THM;UNIV]
    THEN MAP_EVERY EXISTS_TAC [`&0`;`norm (n:real^3)`]
    THEN REWRITE_TAC[VECTOR_MUL_LZERO;VECTOR_ADD_LID]
    THEN MATCH_MP_TAC UNIT_INTRO THEN ASM_MESON_TAC[is_normal_to_plane]]);;

let SNELL_LAW = 
prove
  (`!i emf_i emf_r emf_t.
    is_plane_wave_at_interface i emf_i emf_r emf_t /\ non_null_wave emf_t ==>
    let theta = \emf. vector_angle (k_of_wave emf) (normal_of_interface i) in
    n1_of_interface i * sin (theta emf_i) =
      n2_of_interface i * sin (theta emf_t)`,

  REWRITE_TAC[FORALL_INTERFACE_THM;LET_DEFs;n1_of_interface;n2_of_interface;
    normal_of_interface]
  THEN REPEAT STRIP_TAC
  THEN FIRST_ASSUM (STRIP_ASSUME_TAC o REWRITE_RULE[LET_DEFs;map_triple;o_DEF;
    is_valid_interface;IS_PLANE_WAVE;non_null_wave]
    o CONV_RULE (REWR_CONV is_plane_wave_at_interface))
  THEN SUBGOAL_TAC ""
    `~(k_of_wave emf_i = vec 0) /\ ~(k_of_wave emf_t = vec 0)` 
    [ASM_REWRITE_TAC[]]
  THEN SUBGOAL_TAC "" `~(n = vec 0:real^3)` [ASM_MESON_TAC[is_normal_to_plane]]
  THEN ASSUM_LIST (SIMP_TAC o (@) [REWRITE_RULE[DE_MORGAN_THM;
    MESON[] `(p\/q==>r) <=> (p==>r) /\ (q==>r)`] (CONV_RULE (DEPTH_CONV
      COND_ELIM_CONV) vector_angle)]
    o filter (fun x -> not (contains_sub_term_name "norm" (concl x))))
  THEN SUBGOAL_THEN
    `(k_of_wave emf_i dot unit n) / (k0*n1) =
      (k_of_wave emf_i dot n) / (norm (k_of_wave emf_i) * norm n)
    /\ (k_of_wave emf_t dot unit n) / (k0*n2) =
      (k_of_wave emf_t dot n) / (norm (k_of_wave emf_t) * norm n)`
    ASSUME_TAC
  THENL ON_FIRST_GOAL (ASM_REWRITE_TAC[unit;DOT_RMUL;real_div;
    REAL_ARITH `(x*y)*z=y*(z*x):real`;REAL_INV_MUL])
  THEN SIMP_TAC[REWRITE_RULE[REAL_ABS_BOUNDS] NORM_CAUCHY_SCHWARZ_DIV;SIN_ACS;
    GSYM REAL_EQ_SQUARE_ABS]
  THEN ASM_CSQ_THEN (REWRITE_RULE[LET_DEFs] K_PLANE_PROJECTION_CONSERVATION) 
    (C ASM_CSQ_THEN MP_TAC o CONJUNCT2)
  THEN ASM_SIMP_TAC[normal_of_interface;LET_DEFs;PROJECTION_ON_HYPERPLANE_THM;
    NORMAL_OF_INTERFACE_NON_NULL;UNIT_THM;UNIT_EQ_ZERO]
  THEN DISCH_THEN (MP_TAC o MATCH_MP (MESON[] `x=y ==> x dot x = y dot y`))
  THEN SUBGOAL_THEN `!x:real^3. orthogonal (x - (x dot unit n) % unit n) ((x
    dot unit n) % unit n)` ASSUME_TAC
  THENL ON_FIRST_GOAL (ASM_MESON_TAC[ORTHOGONAL_RUNIT;ORTHOGONAL_CLAUSES;
    SUB_UNIT_NORMAL_IS_ORTHOGONAL_TO_NORMAL])
  THEN ASM_SIMP_TAC[REWRITE_RULE[REAL_ARITH `x=y+z <=> y=x-z:real`] 
    ORTHOGONAL_SQR_NORM;DOT_SQUARE_NORM;NORM_MUL;UNIT_THM;REAL_MUL_RID;
    REAL_POW2_ABS]
  THEN SUBGOAL_TAC "" `~(k0*n2 = &0) /\ ~(k0*n1 = &0)`
    [ASM_MESON_TAC[NORM_EQ_0]]
  THEN SUBGOAL_THEN `&0 < inv(k0*n1)` ASSUME_TAC
  THENL ON_FIRST_GOAL (ASM_REWRITE_TAC[REAL_LT_INV_EQ;REAL_LT_LE]
    THEN ASM_MESON_TAC[NORM_POS_LE])
  THEN SUBGOAL_THEN `~((k0*n2) pow 2 = &0) /\ ~(inv ((k0 * n2) pow 2) = &0)`
    STRIP_ASSUME_TAC
  THENL ON_FIRST_GOAL (REWRITE_TAC[REAL_POW_EQ_0;REAL_INV_EQ_0]
    THEN ASM_ARITH_TAC)
  THEN POP_ASSUM (fun x -> DISCH_THEN (MP_TAC o
    ONCE_REWRITE_RULE [REWRITE_RULE[x] (GENL [`x:real`;`y:real`] (SPECL
      [`x:real`;`y:real`;`inv ((k0 * n2:real) pow 2)`] (GSYM
      REAL_EQ_MUL_RCANCEL)))]))
  THEN ASM_SIMP_TAC[REAL_SUB_RDISTRIB;REAL_MUL_RINV;GSYM real_div;
    GSYM REAL_POW_DIV;REAL_DIV_REFL;REAL_POW_ONE]
  THEN DISCH_THEN (K ALL_TAC)
  THEN MATCH_MP_TAC (MESON[SQRT_MUL;POW_2_SQRT;REAL_LE_POW_2]
    `!x y. &0 <= x /\ &0 <= y /\ &0 <= z /\ &0 <= t /\ sqrt(x pow 2 * y) =
      sqrt(z pow 2*t) ==> x * sqrt y = z * sqrt t`)
  THEN REPEAT CONJ_TAC THENL [
    ASM_REAL_ARITH_TAC;
    REWRITE_TAC[REAL_SUB_LE;ABS_SQUARE_LE_1]
    THEN ASM_MESON_TAC[NORM_CAUCHY_SCHWARZ_DIV];
    ASM_REAL_ARITH_TAC;
    REWRITE_TAC[REAL_POW_DIV;REAL_ARITH `x/y-z/y=(x-z)/y:real`]
    THEN MATCH_MP_TAC REAL_LE_DIV THEN CONJ_TAC THENL [
      REWRITE_TAC[REAL_SUB_LE] THEN MATCH_MP_TAC REAL_POW_LE2 THEN CONJ_TAC
      THENL [
        REWRITE_TAC[unit;DOT_RMUL] THEN MATCH_MP_TAC REAL_LE_MUL
        THEN REWRITE_TAC[REAL_LE_INV_EQ;NORM_POS_LE] THEN ASM_ARITH_TAC;
        FIRST_ASSUM (SINGLE ONCE_REWRITE_TAC o REWRITE_RULE[GSYM real_div]
          o MATCH_MP (GSYM REAL_LE_RMUL_EQ))
        THEN ASM_SIMP_TAC[REAL_DIV_REFL]
        THEN ASM_MESON_TAC[REWRITE_RULE[REAL_ABS_BOUNDS] 
          NORM_CAUCHY_SCHWARZ_DIV]
      ];
      MATCH_ACCEPT_TAC REAL_LE_POW_2;
    ];
    AP_TERM_TAC THEN SUBGOAL_TAC "" `~(k0 = &0) /\ ~(n1 = &0) /\ ~(n2 = &0)`
      [ASM_MESON_TAC[REAL_MUL_LZERO;REAL_MUL_RZERO]]
    THEN ASM_SIMP_TAC[real_div;REAL_INV_MUL;
      REAL_ARITH `(x*y)*inv x*inv z=(x*inv x)*y*inv z:real`;REAL_MUL_RINV;
      REAL_MUL_LID;REAL_SUB_LDISTRIB;REAL_MUL_RID]
    THEN ASM_SIMP_TAC[GSYM REAL_POW_MUL;
      REAL_ARITH `x*y*inv x=(x*inv x)*y:real`;
      REAL_ARITH `x*y*(inv z*inv x)*t=(x*inv x)*y*inv z*t:real`;
      REAL_ARITH `x*y*inv z*inv x=(x*inv x)*y*inv z:real`;
      REAL_MUL_RINV;REAL_MUL_LID;REAL_INV_MUL]
    THEN REWRITE_TAC[REAL_ARITH `x*inv y*inv z=(inv z*x)*inv y:real`;
      GSYM DOT_RMUL;unit]
  ]);;

let phase_shift_at_plane = new_definition
  `phase_shift_at_plane p n emf =
    k_of_wave emf dot (@a. a % unit n IN p) % unit n`;;

let PLANE_WAVE_WITH_PHASE_SHIFT_AT_PLANE = 
prove
  (`!emf. is_plane_wave emf ==> !p. plane p ==> !n. is_normal_to_plane n p ==>
  !r. r IN p ==> !t.
  emf r t = cexp (--ii * Cx(projection_on_hyperplane (k_of_wave emf) (unit n)
    dot r - w_of_wave emf * t + phase_shift_at_plane p n emf))
    % eh_of_wave emf`,

  let tmp = DISCH_ALL (GSYM (UNDISCH (ISPECL [`k_of_wave emf`;`unit n:real^3`]
    PROJECTION_ON_HYPERPLANE_DECOMPOS))) in
  let tmp' = UNDISCH_ALL (IMP_TRANS (SPEC `n:real^3` (UNDISCH (SPEC `p:plane`
    NORMAL_OF_PLANE_NON_NULL))) (IMP_TRANS (ISPEC `n:real^3` UNIT_THM) tmp)) in
  REWRITE_TAC[LET_DEF;LET_END_DEF;IS_PLANE_WAVE;plane_wave]
  THEN REPEAT STRIP_TAC
  THEN FIRST_ASSUM (let thm = MESON[] `f=g ==> g x y=z ==> f x y=z` in
    fun x -> MATCH_MP_TAC (MATCH_MP thm x))
  THEN REWRITE_TAC[emf_at_point_mul;PAIR_EQ;FUN_EQ_THM;eh_of_wave;e_of_emf;
    h_of_emf;LET_DEF;LET_END_DEF]
  THEN REPEAT (FIRST[CONJ_TAC;AP_THM_TAC;AP_TERM_TAC;GEN_TAC])
  THEN REWRITE_TAC[REAL_ARITH `x-y=(z-y)+t <=> x=t+z:real`]
  THEN GEN_REWRITE_TAC (RATOR_CONV o ONCE_DEPTH_CONV) 
    [ONCE_REWRITE_RULE[VECTOR_ADD_SYM] tmp']
  THEN REWRITE_TAC[DOT_LADD;REAL_EQ_ADD_RCANCEL;phase_shift_at_plane]
  THEN ABBREV_TAC `a = @a. a % unit n IN (p:plane)`
  THEN GEN_REWRITE_TAC (RATOR_CONV o ONCE_DEPTH_CONV) 
    [VECTOR_ARITH `v dot r = v dot (r - a % unit n + a % unit n:real^3)`]
  THEN REWRITE_TAC[DOT_RADD;DOT_RMUL;DOT_LMUL;REAL_ADD_LDISTRIB]
  THEN ASM_REWRITE_TAC[MATCH_MP (UNIT_DOT_UNIT_SELF)
    (STRIP_RULE NORMAL_OF_PLANE_NON_NULL);REAL_ARITH `x+y*z* &1=z*y <=> x= &0`]
  THEN REWRITE_TAC[REAL_ENTIRE;GSYM orthogonal;ORTHOGONAL_LUNIT] THEN DISJ2_TAC
  THEN ONCE_REWRITE_TAC[ORTHOGONAL_SYM]
  THEN MATCH_MP_TAC (GENL [`pt1:point`;`pt2:point`] 
    (REWRITE_RULE[GSYM IMP_CONJ] (DISCH_ALL (STRIP_RULE 
      NORMAL_OF_PLANE_IS_ORTHOGONAL_TO_SEGMENT))))
  THEN ASM_REWRITE_TAC[]
  THEN EXPAND_TAC "a" THEN SELECT_ELIM_TAC
  THEN ASM_SIMP_TAC[EXISTS_MULTIPLE_OF_NORMAL_IN_PLANE]);;

let PLANE_WAVE_WITH_PHASE_SHIFT_AT_INTERFACE = 
prove
  (`!emf. is_plane_wave emf ==> !i. is_valid_interface i ==> !r.
  r IN plane_of_interface i ==> !t.
  let n = normal_of_interface i in
  emf r t = cexp (--ii * Cx(projection_on_hyperplane (k_of_wave emf) (unit n)
    dot r - w_of_wave emf * t + phase_shift_at_plane (plane_of_interface i) n emf))
    % eh_of_wave emf`,

  REWRITE_TAC[LET_DEFs] THEN REPEAT STRIP_TAC THEN ASM_CSQ_THEN 
    IS_VALID_INTERFACE_IS_NORMAL_TO_PLANE ASSUME_TAC
  THEN ASM_CSQ_THEN IS_VALID_INTERFACE_PLANE ASSUME_TAC
  THEN ASM_CSQS_THEN PLANE_WAVE_WITH_PHASE_SHIFT_AT_PLANE ASSUME_TAC
  THEN ASM_REWRITE_TAC[]);;

let magnitude_shifted_at_plane = new_definition
  `magnitude_shifted_at_plane p n emf =
    cexp (--ii * Cx(phase_shift_at_plane p n emf)) % eh_of_wave emf`;;

let E_PRESERVED_IN_TE_MODE = 
prove
  (`!i emf_i emf_r emf_t.
    is_plane_wave_at_interface i emf_i emf_r emf_t /\ non_null_wave emf_r
    /\ non_null_wave emf_t /\ TE_mode i emf_i emf_r emf_t ==>
    let magnitude =
      \emf. FST (magnitude_shifted_at_plane (plane_of_interface i)
        (normal_of_interface i) emf) cdot (vector_to_cvector (TE_mode_axis
          i emf_i emf_r emf_t))
    in
    magnitude emf_r + magnitude emf_i = magnitude emf_t`,

  REWRITE_TAC[FORALL_INTERFACE_THM;LET_DEFs;plane_of_interface;
    normal_of_interface] THEN REPEAT STRIP_TAC
  THEN ASM_CSQ_THEN (GEQ_IMP is_plane_wave_at_interface) (STRIP_ASSUME_TAC o
    REWRITE_RULE[LET_DEFs])
  THEN FIRST_ASSUM (let tmp = MESON[non_null_wave] `!emf. non_null_wave emf
    ==> is_plane_wave emf` in ASSUME_TAC o MATCH_MP tmp)
  THEN ASM_CSQ_THEN (GEN_ALL (MATCH_MP EQ_IMP (SPEC_ALL is_valid_interface)))
    (STRIP_ASSUME_TAC o REWRITE_RULE[LET_DEFs])
  THEN ASM_CSQ_THEN PLANE_NON_EMPTY (STRIP_ASSUME_TAC o
    REWRITE_RULE[GSYM MEMBER_NOT_EMPTY])
  THEN FIRST_ASSUM (fun x -> FIRST_X_ASSUM (C ASM_CSQ_THEN (MP_TAC o
    REWRITE_RULE[LET_DEFs] o GEN_REWRITE_RULE (RATOR_CONV o ONCE_DEPTH_CONV)
    [e_of_emf] o CONJUNCT1 o SPEC `&0`) o REWRITE_RULE[boundary_conditions;
    emf_add;LET_DEFs] o C MATCH_MP x))
  THEN ASSUM_LIST (MAP_EVERY (C ASM_CSQS_THEN (SINGLE REWRITE_TAC o CONJUNCT1
    o SPEC `&0` o REWRITE_RULE[emf_at_point_mul;eh_of_wave;FST;SND;EMF_EQ]))
    o mapfilter (REWRITE_RULE[FORALL_INTERFACE_THM;LET_DEFs;plane_of_interface;
      normal_of_interface;]
    o MATCH_MP PLANE_WAVE_WITH_PHASE_SHIFT_AT_INTERFACE))
  THEN ASM_CSQ_THEN (REWRITE_RULE[LET_DEFs] K_PLANE_PROJECTION_CONSERVATION)
    (CONJUNCTS_THEN (C ASM_CSQ_THEN (SINGLE REWRITE_TAC o 
      SIMP_RULE[normal_of_interface;LET_DEFs])))
  THEN REWRITE_TAC[REAL_ARITH `x - y * &0 = x`;CX_ADD;COMPLEX_ADD_LDISTRIB;
    CEXP_ADD;GSYM CVECTOR_MUL_ASSOC;GSYM CVECTOR_ADD_LDISTRIB;
    CCROSS_RMUL;CVECTOR_MUL_LCANCEL;CEXP_NZ;magnitude_shifted_at_plane;
    eh_of_wave;emf_at_point_mul;GSYM CDOT_LADD]
  THEN ASM_CSQS_THEN (REWRITE_RULE[IMP_CONJ;LET_DEFs] TE_MODE_PLANEWAVE_PROJ)
    (GEN_REWRITE_TAC (RATOR_CONV o ONCE_DEPTH_CONV) o CONJUNCTS)
  THEN REWRITE_TAC[CCROSS_RADD;CCROSS_RMUL;CVECTOR_MUL_ASSOC;
    GSYM CVECTOR_ADD_RDISTRIB;CVECTOR_MUL_RCANCEL;CCROSS_COLLINEAR_CVECTORS;
    GSYM CDOT_LMUL;GSYM CDOT_LADD;COLLINEAR_CVECTORS_VECTOR_TO_CVECTOR]
  THEN FIRST_X_ASSUM (STRIP_ASSUME_TAC o REWRITE_RULE[is_mode_axis;
    normal_of_interface;LET_DEFs] o MATCH_MP (GEQ_IMP TE_MODE_THM))
  THEN FIRST_ASSUM (ASSUME_TAC o REWRITE_RULE[NORM_EQ_0] o MATCH_MP 
    (REAL_ARITH `!x. x= &1 ==> ~(x= &0)`))
  THEN FIRST_X_ASSUM (STRIP_ASSUME_TAC o MATCH_MP (GEQ_IMP is_normal_to_plane))
  THEN ASM_CSQS_THEN (ONCE_REWRITE_RULE[ORTHOGONAL_SYM] (REWRITE_RULE[IMP_CONJ]
    ORTHOGONAL_NON_COLLINEAR)) (SINGLE REWRITE_TAC)
  THEN REWRITE_TAC[CVECTOR_ADD_SYM]);;

let H_CROSS_Z_WRT_E_IN_TE_MODE = 
prove
  (`!i emf_i emf_r emf_t.
    is_plane_wave_at_interface i emf_i emf_r emf_t
    /\ TE_mode i emf_i emf_r emf_t ==>
    let p = plane_of_interface i in
    let n = normal_of_interface i in
    let h_cross_z_wrt_e =
      \emf. (h_of_wave emf) ccross (vector_to_cvector n) =
        Cx ((n dot k_of_wave emf) / (eta0 * k0)) % e_of_wave emf
    in
    h_cross_z_wrt_e emf_i /\ h_cross_z_wrt_e emf_r /\ h_cross_z_wrt_e emf_t`,

  SIMP_TAC[FORALL_INTERFACE_THM;LET_DEFs;is_plane_wave_at_interface;
    plane_of_interface;non_null_wave;normal_of_interface]
  THEN REPEAT STRIP_TAC THEN ASM_CSQS_THEN (REWRITE_RULE[LET_DEFs;IMP_CONJ]
    TE_MODE_PLANEWAVE_PROJ) (SINGLE ONCE_REWRITE_TAC)
  THEN REWRITE_TAC[CCROSS_LMUL;CCROSS_LREAL;CDOT_RREAL;CVECTOR_ADD_RDISTRIB;
    GSYM CVECTOR_MUL_ASSOC;GSYM VECTOR_TO_CVECTOR_MUL;
    CVECTOR_IM_VECTOR_TO_CVECTOR_RE_IM;CVECTOR_RE_VECTOR_TO_CVECTOR_RE_IM;
    CCROSS_LADD;CVECTOR_ADD_LDISTRIB]
  THEN REWRITE_TAC[REWRITE_RULE[VECTOR_ARITH `--x = y <=> x = --y :real^N`] 
    (ONCE_REWRITE_RULE[CROSS_SKEW] CROSS_LAGRANGE);
      CVECTOR_RE_VECTOR_TO_CVECTOR;
      CVECTOR_IM_VECTOR_TO_CVECTOR;DOT_LZERO;VECTOR_MUL_LZERO;VECTOR_SUB_RZERO;
      VECTOR_NEG_0;VECTOR_TO_CVECTOR_ZERO;CVECTOR_MUL_RZERO;CVECTOR_ADD_ID;
      DOT_RMUL]
  THEN ASM_CSQ_THEN TE_MODE_AXIS_ORTHOGONAL_N
    (SINGLE REWRITE_TAC o REWRITE_RULE[orthogonal;normal_of_interface;LET_DEFs]
      o ONCE_REWRITE_RULE[ORTHOGONAL_SYM])
  THEN REWRITE_TAC[REAL_MUL_RZERO;VECTOR_MUL_LZERO;VECTOR_ARITH
    `--(vec 0 - x) = x`;VECTOR_NEG_0;VECTOR_TO_CVECTOR_ZERO;CVECTOR_MUL_RZERO;
    CVECTOR_ADD_ID]
  THEN REWRITE_TAC[MESON[CVECTOR_MUL_ASSOC;COMPLEX_MUL_SYM]
    `a % ii % v = ii % a % v:complex^N`;GSYM VECTOR_TO_CVECTOR_MUL;CVECTOR_EQ;
    CVECTOR_RE_VECTOR_TO_CVECTOR_RE_IM;CVECTOR_IM_VECTOR_TO_CVECTOR_RE_IM;
    VECTOR_MUL_ASSOC]
  THEN CONJ_TAC THEN REPEAT (AP_THM_TAC ORELSE AP_TERM_TAC)
  THEN REAL_ARITH_TAC);;

let NON_PROJECTED_E_RELATION_IN_TE_MODE = 
prove
  (`!i emf_i emf_r emf_t.
    is_plane_wave_at_interface i emf_i emf_r emf_t /\ non_null_wave emf_r
    /\ non_null_wave emf_t /\ TE_mode i emf_i emf_r emf_t ==>
    let magnitude =
      \emf. FST (magnitude_shifted_at_plane (plane_of_interface i)
        (normal_of_interface i) emf) cdot (vector_to_cvector (TE_mode_axis i
          emf_i emf_r emf_t))
    in
    let n = unit (normal_of_interface i) in
    Cx (n dot k_of_wave emf_i) * (magnitude emf_i - magnitude emf_r) =
      Cx (n dot k_of_wave emf_t) * magnitude emf_t`,

  REWRITE_TAC[FORALL_INTERFACE_THM;LET_DEFs;normal_of_interface;
    plane_of_interface] THEN REPEAT STRIP_TAC 
  THEN ASM_CSQ_THEN (GEQ_IMP is_plane_wave_at_interface) (STRIP_ASSUME_TAC o
    REWRITE_RULE[LET_DEFs])
  THEN FIRST_ASSUM (let tmp = MESON[non_null_wave] `!emf. non_null_wave emf
    ==> is_plane_wave emf` in ASSUME_TAC o MATCH_MP tmp)
  THEN ASM_CSQ_THEN (GEN_ALL (MATCH_MP EQ_IMP (SPEC_ALL is_valid_interface)))
    (STRIP_ASSUME_TAC o REWRITE_RULE[LET_DEFs])
  THEN ASM_CSQ_THEN PLANE_NON_EMPTY (STRIP_ASSUME_TAC 
    o REWRITE_RULE[GSYM MEMBER_NOT_EMPTY])
  THEN FIRST_ASSUM (fun x -> FIRST_X_ASSUM (C ASM_CSQ_THEN (MP_TAC o
    REWRITE_RULE[LET_DEFs]
    o GEN_REWRITE_RULE (RATOR_CONV o ONCE_DEPTH_CONV) [h_of_emf] 
    o CONJUNCT2 o SPEC `&0`) 
    o REWRITE_RULE[boundary_conditions;emf_add;LET_DEFs] o C MATCH_MP x))
  THEN ASSUM_LIST (MAP_EVERY (C ASM_CSQS_THEN (SINGLE REWRITE_TAC o CONJUNCT2
    o SPEC `&0` o REWRITE_RULE[emf_at_point_mul;eh_of_wave;FST;SND;EMF_EQ]))
    o mapfilter (REWRITE_RULE[FORALL_INTERFACE_THM;LET_DEFs;plane_of_interface;
      normal_of_interface;]
    o MATCH_MP PLANE_WAVE_WITH_PHASE_SHIFT_AT_INTERFACE))
  THEN ASM_CSQ_THEN (REWRITE_RULE[LET_DEFs] K_PLANE_PROJECTION_CONSERVATION)
    (CONJUNCTS_THEN (C ASM_CSQ_THEN (SINGLE REWRITE_TAC
      o SIMP_RULE[normal_of_interface;LET_DEFs])))
  THEN REWRITE_TAC[REAL_ARITH `x - y * &0 = x`;CX_ADD;COMPLEX_ADD_LDISTRIB;
    CEXP_ADD;GSYM CVECTOR_MUL_ASSOC;GSYM CVECTOR_ADD_LDISTRIB;
    CCROSS_RMUL;CVECTOR_MUL_LCANCEL;CEXP_NZ;magnitude_shifted_at_plane;
    eh_of_wave;emf_at_point_mul;GSYM CDOT_LADD;CCROSS_RADD]
  THEN ASM_CSQS_THEN (REWRITE_RULE[IMP_CONJ;CVECTOR_ARITH
    `--x=a%y <=> x = (--a)%y:complex^N`] (ONCE_REWRITE_RULE[CCROSS_SKEW]
    (REWRITE_RULE[LET_DEFs] H_CROSS_Z_WRT_E_IN_TE_MODE)))
    (SINGLE REWRITE_TAC 
      o REWRITE_RULE[LET_DEFs;plane_of_interface;normal_of_interface])
  THEN REWRITE_TAC[GSYM CX_NEG;real_div;REAL_NEG_LMUL;GSYM DOT_RNEG]
  THEN FIRST_ASSUM (ASSUME_TAC o ONCE_REWRITE_RULE[DOT_SYM]
    o REWRITE_RULE[DOT_RUNIT_EQ] o MATCH_MP SYMETRIC_VECTORS_PROJECTION_ON_AXIS
    o MATCH_MP UNIT_THM o CONJUNCT1 o MATCH_MP (GEQ_IMP is_normal_to_plane))
  THEN ASM_CSQS_THEN (SIMP_RULE[LET_DEFs] LAW_OF_REFLECTION) (fun x -> 
    POP_ASSUM (SINGLE REWRITE_TAC o C MATCH_MP (REWRITE_RULE
      [normal_of_interface;LET_DEFs] x)))
  THEN REWRITE_TAC[CX_MUL;DOT_RNEG;CX_NEG;COMPLEX_MUL_LNEG;CVECTOR_ARITH
    `a%(--(u*v))%x+c%(u*v)%y = d%(--(u'*v))%z
    <=> v%u%(a%x-c%y) = v%u'%d%z:complex^N`]
  THEN REWRITE_TAC[CVECTOR_MUL_LCANCEL;CX_INJ;REAL_INV_EQ_0;REAL_ENTIRE;
    MATCH_MP REAL_LT_IMP_NZ eta0;MATCH_MP REAL_LT_IMP_NZ k0]
  THEN REWRITE_TAC[unit;DOT_LMUL;CX_MUL;GSYM COMPLEX_MUL_ASSOC;
    COMPLEX_EQ_MUL_LCANCEL;CX_INJ;REAL_ENTIRE;REAL_INV_EQ_0;NORM_EQ_0]
  THEN ASM_CSQS_THEN NORMAL_OF_PLANE_NON_NULL (SINGLE REWRITE_TAC)
  THEN ASM_CSQS_THEN (REWRITE_RULE[IMP_CONJ;LET_DEFs] TE_MODE_PLANEWAVE_PROJ)
    (GEN_REWRITE_TAC (RATOR_CONV o ONCE_DEPTH_CONV) o CONJUNCTS)
  THEN REWRITE_TAC[CVECTOR_MUL_ASSOC;GSYM CVECTOR_SUB_RDISTRIB;
    CVECTOR_MUL_RCANCEL;VECTOR_TO_CVECTOR_ZERO_EQ;GSYM NORM_EQ_0]
  THEN FIRST_ASSUM (SINGLE REWRITE_TAC o MATCH_MP (REWRITE_RULE[is_mode_axis]
    (GEQ_IMP TE_MODE_THM)))
  THEN SIMP_TAC[REAL_ARITH `~(&1 = &0)`;CDOT_LMUL;COMPLEX_MUL_ASSOC]);;

let COMPLEX_MUL_LINV2 = 
prove
  (`!x y. ~(x=Cx(&0)) ==> inv x * (x * y) = y`,

  SIMP_TAC[COMPLEX_MUL_ASSOC;COMPLEX_MUL_LINV;COMPLEX_MUL_LID]);;

let COMPLEX_BALANCE_MUL_DIV = 
prove
  (`!x y z. ~(x=Cx(&0)) ==> (x*y=z <=> y=z/x)`,

 REPEAT STRIP_TAC THEN EQ_TAC THENL [
   DISCH_THEN (fun x -> ASM_SIMP_TAC[GSYM x;complex_div;SIMPLE_COMPLEX_ARITH
   `(x*y)*inv x=(x*inv x)*y:complex`;COMPLEX_MUL_RINV;COMPLEX_MUL_LID]);
   DISCH_THEN (fun x -> ASM_SIMP_TAC[x;COMPLEX_DIV_LMUL]);
 ]);;


let FRESNEL_REFLECTION_TE_MODE = 
prove
  (`!i emf_i emf_r emf_t.
  is_plane_wave_at_interface i emf_i emf_r emf_t /\ non_null_wave emf_r
  /\ non_null_wave emf_t /\ TE_mode i emf_i emf_r emf_t ==>
  let magnitude =
    \emf. FST (magnitude_shifted_at_plane (plane_of_interface i)
      (normal_of_interface i) emf) cdot (vector_to_cvector (TE_mode_axis i
        emf_i emf_r emf_t))
  in
  let theta = \emf. Cx(vector_angle (k_of_wave emf) (normal_of_interface i)) in
  let n1 = Cx(n1_of_interface i) in
  let n2 = Cx(n2_of_interface i) in
  magnitude emf_r = (n1 * ccos (theta emf_i) - n2 * ccos (theta emf_t)) / (n1 *
    ccos (theta emf_i) + n2 * ccos (theta emf_t)) * magnitude emf_i`,

  REPEAT STRIP_TAC THEN ASM_CSQS_THEN (REWRITE_RULE[IMP_CONJ]
    NON_PROJECTED_E_RELATION_IN_TE_MODE) MP_TAC
  THEN ASM_CSQS_THEN (REWRITE_RULE[IMP_CONJ] E_PRESERVED_IN_TE_MODE) (MP_TAC o
    GSYM)
  THEN LET_TAC THEN REWRITE_TAC[LET_DEFs] 
  THEN DISCH_THEN (fun x -> REWRITE_TAC[x;SIMPLE_COMPLEX_ARITH
    `x*(y-z)=t*(z+y) <=> (x+t)*z=(x-t)*y:complex`])
  THEN DISCH_THEN (MP_TAC o MATCH_MP (MATCH_MP (MESON[]
    `(!x y z. p x ==> (q x y z <=> r x y z)) ==> !x y z. q x y z ==> p x
    ==> r x y z`) COMPLEX_BALANCE_MUL_DIV))
  THEN ANTS_TAC THENL [
    REWRITE_TAC[GSYM CX_ADD;CX_INJ] THEN MATCH_MP_TAC REAL_POS_NZ
    THEN MATCH_MP_TAC REAL_LTE_ADD THEN REWRITE_TAC[DOT_LUNIT_POS;DOT_LUNIT_GE0]
    THEN ONCE_REWRITE_TAC[DOT_SYM] THEN REPEAT (POP_ASSUM MP_TAC)
    THEN SPEC_TAC (`i:interface`,`i:interface`)
    THEN SIMP_TAC[FORALL_INTERFACE_THM;is_plane_wave_at_interface;
      normal_of_interface;map_triple;LET_DEFs;real_ge;real_gt];
    SIMP_TAC[SIMPLE_COMPLEX_ARITH `(x*y)/z = (x/z)*y`]
    THEN DISCH_THEN (K ALL_TAC) THEN AP_THM_TAC THEN AP_TERM_TAC 
    THEN REWRITE_TAC[GSYM CX_ADD;GSYM CX_SUB;GSYM CX_COS;GSYM CX_MUL;
      GSYM CX_DIV;CX_INJ]
    THEN REPEAT (POP_ASSUM MP_TAC) THEN SPEC_TAC (`i:interface`,`i:interface`)
    THEN REWRITE_TAC[FORALL_INTERFACE_THM;normal_of_interface;
      plane_of_interface;LET_DEFs;n1_of_interface;n2_of_interface]
    THEN REPEAT STRIP_TAC
    THEN ASM_CSQ_THEN (GEQ_IMP is_plane_wave_at_interface) (STRIP_ASSUME_TAC
      o REWRITE_RULE[is_valid_interface;is_normal_to_plane;LET_DEFs;
        map_triple])
    THEN ASM_REWRITE_TAC[GSYM NORM_EQ_0]
    THEN ASM_SIMP_TAC [GEN_ALL (DISCH_ALL (GSYM (MATCH_MP REAL_LT_IMP_NE
      (UNDISCH_ALL (SPEC_ALL (MATCH_MP (REWRITE_RULE[IMP_CONJ] REAL_LT_MUL)
        k0))))))]
    THEN GEN_REWRITE_TAC (RATOR_CONV o ONCE_DEPTH_CONV) [VECTOR_ANGLE]
    THEN ASM_SIMP_TAC[UNIT_THM;REAL_MUL_LID]
    THEN REWRITE_TAC[GSYM REAL_MUL_ASSOC;GSYM REAL_ADD_LDISTRIB;
      GSYM REAL_SUB_LDISTRIB;real_div;REAL_INV_MUL]
    THEN REWRITE_TAC[REAL_ARITH `x*y*inv x*z=(x*inv x)*y*z:real`;
      MATCH_MP REAL_MUL_RINV (GSYM (MATCH_MP REAL_LT_IMP_NE k0));REAL_MUL_LID]
    THEN ASM_SIMP_TAC[unit;VECTOR_ANGLE_LMUL;REAL_INV_EQ_0;NORM_EQ_0;
      REAL_LE_INV_EQ;NORM_POS_LE;VECTOR_ANGLE_SYM]]);;

let FRESNEL_TRANSMISSION_TE_MODE = 
prove
  (`!i emf_i emf_r emf_t.
    is_plane_wave_at_interface i emf_i emf_r emf_t /\ non_null_wave emf_r 
    /\ non_null_wave emf_t /\ TE_mode i emf_i emf_r emf_t ==>
    let magnitude =
      \emf. FST (magnitude_shifted_at_plane (plane_of_interface i) 
        (normal_of_interface i) emf) cdot (vector_to_cvector (TE_mode_axis i
          emf_i emf_r emf_t))
    in
    let theta =
      \emf. Cx(vector_angle (k_of_wave emf) (normal_of_interface i))
    in
    let n1 = Cx(n1_of_interface i) in
    let n2 = Cx(n2_of_interface i) in
    magnitude emf_t = Cx(&2) * n1 * ccos (theta emf_i) / (n1 * ccos
      (theta emf_i) + n2 * ccos (theta emf_t)) * magnitude emf_i`,

  REPEAT STRIP_TAC THEN ASM_CSQS_THEN (REWRITE_RULE[IMP_CONJ] 
    E_PRESERVED_IN_TE_MODE) (MP_TAC o GSYM)
  THEN ASM_CSQS_THEN (REWRITE_RULE[IMP_CONJ] 
    NON_PROJECTED_E_RELATION_IN_TE_MODE) MP_TAC THEN LET_TAC
    THEN REWRITE_TAC[LET_DEFs] 
  THEN REWRITE_TAC[SIMPLE_COMPLEX_ARITH `x*(y-z)=t*u <=> x*z=x*y-t*u:complex`]
  THEN DISCH_THEN (MP_TAC o MATCH_MP (MATCH_MP (MESON[]
    `(!x y z. p x ==> (q x y z <=> r x y z)) ==> !x y z. q x y z ==> p x
      ==> r x y z`) COMPLEX_BALANCE_MUL_DIV))
  THEN MATCH_MP_TAC (MESON[] `p/\(p==>q==>r==>s) ==> ((p==>q)==>r==>s)`)
  THEN CONJ_TAC THENL [
    REWRITE_TAC[CX_INJ] THEN MATCH_MP_TAC REAL_POS_NZ
    THEN REWRITE_TAC[DOT_LUNIT_POS;DOT_LUNIT_GE0]
    THEN ONCE_REWRITE_TAC[DOT_SYM] THEN REPEAT (POP_ASSUM MP_TAC)
    THEN SPEC_TAC (`i:interface`,`i:interface`)
    THEN SIMP_TAC[FORALL_INTERFACE_THM;is_plane_wave_at_interface;
      normal_of_interface;map_triple;LET_DEFs;real_gt];
    DISCH_TAC THEN DISCH_THEN (fun x -> ASM_SIMP_TAC[x;SIMPLE_COMPLEX_ARITH
      `t = (x*i-y*t)/x+i <=> (y/x+Cx(&1))*t = (x/x+Cx(&1))*i:complex`;
      COMPLEX_DIV_REFL;SIMPLE_COMPLEX_ARITH `Cx(&1)+Cx(&1)=Cx(&2)`])
    THEN ASM_CSQ_THEN COMPLEX_DIV_REFL (SINGLE REWRITE_TAC o GSYM)
    THEN REWRITE_TAC[GSYM CX_ADD;GSYM CX_DIV;real_div;GSYM REAL_ADD_RDISTRIB]
    THEN DISCH_THEN (MP_TAC o MATCH_MP (MATCH_MP (MESON[]
      `(!x y z. p x ==> (q x y z <=> r x y z)) ==> (!x y z. q x y z ==> p x 
        ==> r x y z)`) COMPLEX_BALANCE_MUL_DIV))
    THEN ANTS_TAC THENL [
      ASM_REWRITE_TAC[CX_MUL;CX_INV;COMPLEX_ENTIRE;COMPLEX_INV_EQ_0] 
      THEN REWRITE_TAC[CX_INJ] THEN MATCH_MP_TAC REAL_POS_NZ 
      THEN MATCH_MP_TAC REAL_LET_ADD
      THEN REWRITE_TAC[DOT_LUNIT_POS;DOT_LUNIT_GE0]
      THEN ONCE_REWRITE_TAC[DOT_SYM] THEN REPEAT (POP_ASSUM MP_TAC)
      THEN SPEC_TAC (`i:interface`,`i:interface`)
      THEN SIMP_TAC[FORALL_INTERFACE_THM;is_plane_wave_at_interface;
        normal_of_interface;map_triple;LET_DEFs;real_ge;real_gt];
      SIMP_TAC[SIMPLE_COMPLEX_ARITH `(x*y)/z = x/z*y:complex`;GSYM CX_DIV;
        real_div;REAL_INV_MUL;REAL_INV_INV] THEN DISCH_THEN (K ALL_TAC)
      THEN REWRITE_TAC[GSYM CX_MUL;GSYM CX_COS;GSYM CX_ADD;GSYM CX_DIV;
        COMPLEX_MUL_ASSOC;CX_INJ]
      THEN AP_THM_TAC THEN AP_TERM_TAC
      THEN REWRITE_TAC[CX_INJ;REAL_ARITH `&2*x = (&2*y)*z <=> x=y*z`]
      THEN REPEAT (POP_ASSUM MP_TAC) THEN SPEC_TAC (`i:interface`,`i:interface`)
      THEN REWRITE_TAC[FORALL_INTERFACE_THM;normal_of_interface;
        plane_of_interface;LET_DEFs;n1_of_interface;n2_of_interface]
      THEN REPEAT STRIP_TAC
      THEN ASM_CSQ_THEN (GEQ_IMP is_plane_wave_at_interface) (STRIP_ASSUME_TAC
        o REWRITE_RULE[is_valid_interface;is_normal_to_plane;LET_DEFs;map_triple])
      THEN GEN_REWRITE_TAC (RATOR_CONV o ONCE_DEPTH_CONV) [VECTOR_ANGLE]
      THEN ASM_SIMP_TAC[UNIT_THM;REAL_MUL_LID]
      THEN REWRITE_TAC[GSYM REAL_MUL_ASSOC;GSYM REAL_ADD_LDISTRIB;
        GSYM REAL_SUB_LDISTRIB;real_div;REAL_INV_MUL]
      THEN REWRITE_TAC[REAL_ARITH `inv x*y*x*z=(x*inv x)*y*z:real`;
        MATCH_MP REAL_MUL_RINV (GSYM (MATCH_MP REAL_LT_IMP_NE k0));
          REAL_MUL_LID]
      THEN ASM_SIMP_TAC[unit;VECTOR_ANGLE_LMUL;REAL_INV_EQ_0;NORM_EQ_0;
        REAL_LE_INV_EQ;NORM_POS_LE]
      THEN GEN_REWRITE_TAC (RAND_CONV o ONCE_DEPTH_CONV) [VECTOR_ANGLE_SYM]
      THEN REAL_ARITH_TAC]]);;