(* ========================================================================= *)
(*                                                                           *)
(*                        Geometrical Optics Library                         *)
(*                                                                           *)
(*        (c) Copyright, Muhammad Umair Siddique, Vincent Aravantinos, 2012  *)
(*                       Hardware Verification Group,                        *)
(*                       Concordia University                                *)
(*                                                                           *)
(*     Contact: <muh_sidd@ece.concordia.ca>, <vincent@ece.concordia.ca>      *)
(*                                                                           *)
(* Last update: Feb 16, 2012                                                 *)
(*                                                                           *)
(* ========================================================================= *)


needs "Multivariate/make_complex.ml";;
needs "utils.ml";;
(* ------------------------------------------------------------------------- *)
(* Preliminaries                                                             *)
(* ------------------------------------------------------------------------- *)

let LET_PAIR = prove
  (`!P t. (let x,y = t in P x y) = P (FST t) (SND t)`,
  ONCE_REWRITE_TAC[GSYM PAIR] THEN PURE_REWRITE_TAC[LET_DEF;LET_END_DEF;FST;SND]
  THEN REWRITE_TAC[]);;

let ALL_REVERSE = prove
  (`!l P. ALL P (REVERSE l) = ALL P l`,
  LIST_INDUCT_TAC THEN ASM_REWRITE_TAC[REVERSE;ALL;ALL_APPEND] THEN ITAUT_TAC);;

let BOUNDED_VECTOR_COMPONENT = prove
  (`!v:real^N. (?B:real^N. !i. 1 <= i /\ i <= dimindex(:N) ==> v$i <= B$i)
      <=> (?b. !i. 1 <= i /\ i <= dimindex(:N) ==> v$i <= b)`,
  REPEAT (STRIP_TAC ORELSE EQ_TAC) THENL [
    MESON_TAC[REAL_LE_TRANS;REAL_ABS_LE;COMPONENT_LE_INFNORM];
    EXISTS_TAC `b % vec 1:real^N`
    THEN ASM_REWRITE_TAC[VEC_COMPONENT;VECTOR_MUL_COMPONENT;REAL_MUL_RID]
  ]);;

let BOUNDED_MATRIX_BOUNDED_VECTOR_MUL = prove
  (`!M:num->real^M^N.
      (?b. !n i j. 1 <= i /\ i <= dimindex(:N) /\ 1 <= j /\ j <= dimindex(:M)
          ==> abs ((M n)$i$j) <= b)
      ==> (!v:real^M.
        ?b. !n i. 1 <= i /\ i <= dimindex(:N) ==> abs ((M n ** v)$i) <= b)`,
  REPEAT STRIP_TAC
  THEN SIMP_TAC[MATRIX_VECTOR_MUL_COMPONENT;dot]
  THEN Pa.EXISTS_TAC "sum (1..dimindex(:M)) (\i. b * abs (v$i))"
  THEN REPEAT GEN_TAC THEN DISCH_TAC THEN MATCH_MP_TAC SUM_ABS_LE
  THEN REWRITE_TAC[FINITE_NUMSEG;IN_NUMSEG;REAL_ABS_MUL]
  THEN GEN_TAC THEN Pa.ASM_CASES_TAC "v$x = &0"
  THEN ASM_REWRITE_TAC[REAL_ABS_0;REAL_MUL_RZERO;REAL_LE_REFL]
  THEN MP_REWRITE_TAC REAL_LE_RMUL_EQ
  THEN ASM_MESON_TAC[REAL_LT_LE;REAL_ABS_ZERO;REAL_ABS_POS]);;


(* ------------------------------------------------------------------------- *)
(* System formalization                                                      *)
(* ------------------------------------------------------------------------- *)

new_type_abbrev ("free_space",`:real#real`);;

let FORALL_FREE_SPACE_THM = prove
  (`!P. (!fs:free_space. P fs) <=> !n d. P (n,d)`,
  REWRITE_TAC[FORALL_PAIR_THM]);;

let EXISTS_FREE_SPACE_THM = prove
  (`!P. (?fs:free_space. P fs) <=> ?n d. P (n,d)`,
  REWRITE_TAC[EXISTS_PAIR_THM]);;

let interface_IND,interface_REC = define_type
  "interface = plane | spherical real";;

let interface_kind_IND,interface_kind_REC = define_type
  "interface_kind = transmitted | reflected";;

new_type_abbrev("optical_component",`:free_space#interface#interface_kind`);; 

let FORALL_OPTICAL_COMPONENT_THM = prove
  (`!P. (!c:optical_component. P c) <=> !fs i ik. P (fs,i,ik)`,
  REWRITE_TAC[FORALL_PAIR_THM]);;

let EXISTS_OPTICAL_COMPONENT_THM = prove
  (`!P. (?c:optical_component. P c) <=> ?fs i ik. P (fs,i,ik)`,
  REWRITE_TAC[EXISTS_PAIR_THM]);;

new_type_abbrev("optical_system",`:optical_component list#free_space`);;

let FORALL_OPTICAL_SYSTEM_THM = prove
  (`!P. (!os:optical_system. P os) <=> !cs fs. P (cs,fs)`,
  REWRITE_TAC[FORALL_PAIR_THM]);;

let is_valid_free_space = define
  `is_valid_free_space ((n,d):free_space) <=> &0 < n /\ &0 <= d`;;

let is_valid_interface = define
  `(is_valid_interface (plane) <=> T) /\
   (is_valid_interface (spherical R) <=> ~(&0 = R))`;;

let is_valid_optical_component = new_definition
  `is_valid_optical_component (fs,i,ik) <=>
    is_valid_free_space fs /\ is_valid_interface i`;;

let is_valid_optical_system = new_definition
  `is_valid_optical_system (cs,fs) <=>
    ALL is_valid_optical_component cs /\ is_valid_free_space fs`;;

let head_index = define
  `head_index ([],(n,d)) = n /\ 
  head_index (CONS ((n,d),i) cs, (nt,dt)) = n`;;

let HEAD_INDEX_VALID_OPTICAL_SYSTEM = prove
  (`!os:optical_system. is_valid_optical_system os ==> &0 < head_index os`,
  REWRITE_TAC[FORALL_OPTICAL_SYSTEM_THM;is_valid_optical_system]
  THEN MATCH_MP_TAC list_INDUCT
  THEN SIMP_TAC[FORALL_OPTICAL_COMPONENT_THM;FORALL_FREE_SPACE_THM;
    head_index;is_valid_free_space;ALL;is_valid_optical_component]);;

let HEAD_INDEX_APPEND = prove
  (`!cs1:optical_component list cs2 fs a. 
      head_index (APPEND cs1 cs2, fs)
        = head_index (cs1, (head_index (cs2,fs),a))`,
    MATCH_MP_TAC list_INDUCT
    THEN REWRITE_TAC[APPEND;head_index;FORALL_OPTICAL_COMPONENT_THM;
      FORALL_FREE_SPACE_THM]);;

let HEAD_INDEX_CONS_ALT = prove
  (`!fs i cs:optical_component list fs'.
      head_index (CONS (fs,i) cs,fs')
        = head_index ([]:optical_component list,fs)`,
  REWRITE_TAC[FORALL_FREE_SPACE_THM;head_index]);;

let HEAD_INDEX_CONS_EQ = prove
  (`!n1 n2 d1 d2 i1 i2 cs1 cs2 fs1 fs2.
      head_index (CONS ((n1,d1),i1) cs1,fs1)
        = head_index (CONS ((n2,d2),i2) cs2,fs2)
          <=> n1 = n2`,
  REWRITE_TAC[FORALL_PAIR_THM;head_index;PAIR_EQ]);;
  
let HEAD_INDEX_CONS_EQ_ALT = prove
  (`!c1 c2 cs1 cs2 fs1 fs2.
      c1 = c2 ==> head_index (CONS c1 cs1,fs1) = head_index (CONS c2 cs2,fs2)`,
  SIMP_TAC[FORALL_PAIR_THM;head_index;PAIR_EQ]);;

(* ------------------------------------------------------------------------- *)
(* Ray formalization                                                         *)
(* ------------------------------------------------------------------------- *)

new_type_abbrev("single_ray",`:real#real`);;

let FORALL_SINGLE_RAY_THM = prove
  (`!P. (!sr:single_ray. P sr) <=> !y theta. P (y,theta)`,
  REWRITE_TAC[FORALL_PAIR_THM]);;

new_type_abbrev("ray",
  `:single_ray # single_ray # (single_ray # single_ray) list`);;

let FORALL_RAY_THM = prove
  (`!P. (!r:ray. P r) <=> !sr1 sr2 srs. P (sr1,sr2,srs)`,
  REWRITE_TAC[FORALL_PAIR_THM]);;

let is_valid_ray_in_free_space = define 
  `is_valid_ray_in_free_space (y0,theta0) (y1,theta1) ((n0,d):free_space) <=>
    y1 = y0 + d * theta0 /\ theta0 = theta1`;;

let is_valid_ray_at_interface = define 
  `(is_valid_ray_at_interface (y0,theta0) (y1,theta1) n0 n1 plane transmitted
    <=> y1 = y0 /\ n0 * theta0 =  n1 * theta1) /\
  (is_valid_ray_at_interface (y0,theta0) (y1,theta1) n0 n1 (spherical R)
    transmitted <=> 
      let phi_i = theta0 + y1 / R and phi_t = theta1 + y1 / R in
      y1 = y0 /\ n0 * phi_i = n1 * phi_t) /\
  (is_valid_ray_at_interface (y0,theta0) (y1,theta1) n0 n1 plane reflected <=>
    y1 = y0 /\ n0 * theta0 = n0 * theta1) /\
  (is_valid_ray_at_interface (y0,theta0) (y1,theta1) n0 n1 (spherical R)
    reflected <=> 
      let phi = y1 / R - theta0  in
      y1 = y0 /\ theta1 = --(theta0 + &2 * phi))`;;

let last_single_ray = define
  `last_single_ray ((sr1,sr2,[]):ray) = sr2 /\
  last_single_ray ((sr1,sr2,CONS sr3 t):ray) = SND (LAST (CONS sr3 t))`;;

let fst_single_ray = define
  `fst_single_ray ((sr1,sr2,rs):ray )= sr1`;; 


let EXISTS_IS_VALID_RAY_IN_SYSTEM = prove(
  `?f. !sr1 sr2 h h' fs cs rs i ik y0 theta0 y1 theta1 y2 theta2 y3 theta3 n d n' d'.
  (f (sr1,sr2,[]) (CONS h cs,fs) <=> F) /\
  (f (sr1,sr2, CONS h' rs) ([],fs) <=> F) /\
  (f ((y0,theta0),(y1,theta1), []) ([],(n,d)) <=>
    is_valid_ray_in_free_space (y0,theta0) (y1,theta1) (n,d)) /\
  (f ((y0,theta0),(y1,theta1), CONS ((y2,theta2),(y3,theta3)) rs) 
    ((CONS ((n',d'),i,ik) cs),(n,d)) <=>
  (is_valid_ray_in_free_space (y0,theta0) (y1,theta1) (n',d') /\ 
  (is_valid_ray_at_interface (y1,theta1) (y2,theta2) n' (head_index(cs,(n,d))) i ik))
  /\ f ((y2,theta2),(y3,theta3),rs) (cs,(n,d)))`,
  SUBGOAL_THEN `!y0 theta0 y1 theta1 n d. ?f1. (f1 ([]:(free_space#interface#interface_kind)list) = is_valid_ray_in_free_space (y0,theta0) (y1,theta1) (n,d)
    /\ !h t. f1 (CONS h t) = F)` (STRIP_ASSUME_TAC o REWRITE_RULE[SKOLEM_THM])
  THENL [
    REPEAT GEN_TAC THEN MATCH_ACCEPT_TAC list_RECURSION;
    SUBGOAL_THEN 
      `!y2 theta2 self y0 theta0 y1 theta1 (fs:free_space). ?f2. (f2 [] = F /\ !h t.
      f2 (CONS h t) <=>
      let (n,d),i,ik = h in
      is_valid_ray_in_free_space (y0,theta0) (y1,theta1) (n,d) /\ 
      is_valid_ray_at_interface (y1,theta1) (y2,theta2) n (head_index(t,fs)) i ik /\
      self t)`
    (STRIP_ASSUME_TAC o REWRITE_RULE[SKOLEM_THM])
    THENL [
      REPEAT GEN_TAC THEN MATCH_ACCEPT_TAC list_RECURSION;
      SUBGOAL_THEN
        `!n:real d:real. ?g. (g [] = (\y0:real theta0:real y1:real theta1:real. f1 y0 theta0 y1 theta1 n d) /\ !h t. g (CONS h t) =
        \y0:real theta0:real y1:real theta1:real.
        let (y2,theta2),(y3,theta3) = (h:single_ray#single_ray) in
        (f2 y2 theta2 (g t y2 theta2 y3 theta3) y0 theta0 y1 theta1 (n,d)):(free_space#interface#interface_kind)list ->bool)`
        (STRIP_ASSUME_TAC o REWRITE_RULE [SKOLEM_THM])
      THENL [
        REPEAT GEN_TAC THEN MATCH_ACCEPT_TAC list_RECURSION;
        EXISTS_TAC `\(((y0,theta0),(y1,theta1),rs):ray) ((cs,(n,d)):optical_system). (g n d rs y0 theta0 y1 theta1 cs):bool`
        THEN ASM_REWRITE_TAC[LET_DEF;LET_END_DEF;CONJ_ACI;FORALL_PAIR_THM]
      ]
    ]
  ]);;

let is_valid_ray_in_system = new_specification ["is_valid_ray_in_system"] EXISTS_IS_VALID_RAY_IN_SYSTEM;;


(* ------------------------------------------------------------------------- *)
(* 2x2 Matrices                                                              *)
(* ------------------------------------------------------------------------- *)

let mat_pow = define
  `mat_pow (M:real^N^N) (0:num) :real^N^N = mat 1 /\ 
   mat_pow M (SUC n) = M ** mat_pow M n`;;

overload_interface("pow",`mat_pow:real^N^N->num->real^N^N`);;

let MAT_POW_ALT = prove
  (`!n M. mat_pow M 0 = mat 1 /\ mat_pow M (SUC n) = mat_pow M n ** M`,
  REWRITE_TAC[mat_pow] THEN INDUCT_TAC
  THEN REWRITE_TAC[mat_pow;MATRIX_MUL_LID;MATRIX_MUL_RID]
  THEN ASM_REWRITE_TAC[GSYM MATRIX_MUL_ASSOC]);;

let mat2x2 = new_definition
  `mat2x2 A B C D :A^2^2 = vector[ vector[A;B]; vector[C;D] ]`;;

let COMMON_TAC xs =
  SIMP_TAC (xs @ [mat2x2;CART_EQ;VECTOR_2;DIMINDEX_2;FORALL_2;DOT_2;SUM_2]);;

let MAT2X2_VECTOR_MUL = prove
  (`!A B C D X Y.
    mat2x2 A B C D ** vector [X;Y] = vector[A*X+B*Y;C*X+D*Y]`,
  COMMON_TAC [MATRIX_VECTOR_MUL_COMPONENT]);;

let MAT2X2_VECTOR_MUL_GEN = prove
  (`!(M:real^2^2) (X:real) (Y:real).
    M ** vector [X;Y] = vector[(M$1$1)*X+ (M$1$2)*Y;(M$2$1*X)+ (M$2$2)*Y]`,
  COMMON_TAC [MATRIX_VECTOR_MUL_COMPONENT]);;

let MAT2X2_EQ = prove
  (`!A B C D P Q R S. mat2x2 A B C D = mat2x2 P Q R S <=>
    A = P /\ B = Q /\ C = R /\ D = S`,
  COMMON_TAC [] THEN CONV_TAC TAUT);;  

let VECTOR2_EQ = prove
  (`!x y z t. vector [x;y] :A^2 = vector [z;t] <=> x=z /\ y=t`,
  COMMON_TAC[]);;

let MAT2X2_MUL = prove
  (`!A B C D P Q R S.
    mat2x2 A B C D ** mat2x2 P Q R S =
      mat2x2 (A*P + B*R) (A*Q + B*S) (C*P + D*R) (C*Q + D*S)`,
  COMMON_TAC[matrix_mul;MATRIX_VECTOR_MUL_COMPONENT;LAMBDA_BETA]);;

let MAT2X2_MAT1 = prove
  (`mat2x2 (&1) (&0) (&0) (&1) = mat 1`,
  COMMON_TAC[mat;LAMBDA_BETA] THEN CONV_TAC NUM_REDUCE_CONV);;

let MAT2X2_LID = prove
  (`!m:real^2^2. mat2x2 (&1) (&0) (&0) (&1) ** m = m`,
  REWRITE_TAC[MAT2X2_MAT1;MATRIX_MUL_LID]);;

let MAT2X2_RID = prove
  (`!m:real^2^2. m ** mat2x2 (&1) (&0) (&0) (&1) = m`,
  REWRITE_TAC[MAT2X2_MAT1;MATRIX_MUL_RID]);;

let MAT2X2_ID = CONJ MAT2X2_LID MAT2X2_RID;;

let MAT2X2_SCALAR_MUL = prove
  (`!A B C D k. k %% mat2x2 A B C D = mat2x2 (k*A) (k*B) (k*C) (k*D)`,
  COMMON_TAC [MATRIX_CMUL_COMPONENT]);;

let MAT2X2_SCALAR_MUL_ASSOC = prove
  (`!A B C D L M N P k.
    mat2x2 L M N P ** (k %% mat2x2 A B C D)
      = k %% (mat2x2 L M N P ** mat2x2 A B C D)`,
  SIMP_TAC[MAT2X2_SCALAR_MUL;MAT2X2_MUL;MAT2X2_EQ] THEN CONV_TAC REAL_FIELD);;

let MAT2X2_EQ_SCALAR_LMUL = prove
  (`!A B C D P Q R S. ~(k = &0) ==>
    (k %% mat2x2 A B C D = k %% mat2x2 P Q R S
      <=> A = P /\ B = Q /\ C = R /\ D = S)`,
  COMMON_TAC[MAT2X2_SCALAR_MUL] THEN CONV_TAC REAL_FIELD);;

let MAT2X2_DET = prove
  (`!A B C D. det (mat2x2 A B C D) = A * D - B * C`,
  COMMON_TAC[DET_2]);;

let FORALL_MATRIX_2X2 = prove
  (`!P. (!M:real^2^2. P M) <=> ! a b c d.P (mat2x2 a b c d )`,
 REWRITE_TAC[FORALL_VECTOR_2;mat2x2] );;

let MAT2X2_VECTOR_MUL_ALT = prove 
 (`!(M:real^2^2) (x:real) y x' y'.  (vector[x';y'] = M ** vector[x;y]) =    
                   (x' =  (M ** vector[x;y])$1 /\ y' = (M ** vector[x;y])$2)`,
COMMON_TAC [ FORALL_MATRIX_2X2;VECTOR2_EQ;MAT2X2_VECTOR_MUL]);;


let MAT2_MUL = prove (
 ` !(M:real^2^2). M**M =  mat2x2((M$1$1 * M$1$1) +  (M$1$2 * M$2$1))
                                ((M$1$1 * M$1$2) +  (M$1$2 * M$2$2))
                                ((M$2$1 * M$1$1) +  (M$2$2 * M$2$1))
                                ((M$2$1 * M$1$2) +  (M$2$2 * M$2$2))`,
 COMMON_TAC[matrix_mul;MATRIX_VECTOR_MUL_COMPONENT;LAMBDA_BETA]);;

 
let MAT_POW_ADD = prove (
`!(M:real^a^a) m n. (M pow (m + n)) = ((M pow m)**(M pow n))`,
GEN_TAC THEN INDUCT_TAC THEN 
REWRITE_TAC[mat_pow;ADD_CLAUSES;MATRIX_MUL_LID] THEN
REWRITE_TAC[mat_pow;ADD_CLAUSES;MATRIX_MUL_LID] THEN
REWRITE_TAC[GSYM MATRIX_MUL_ASSOC] THEN ASM_SIMP_TAC[]);;


let MAT2X2_POW2_POW = prove (
 ` !N. ((mat2x2 A B C D) pow 2) pow N = (mat2x2 A B C D) pow (2*N)`,
INDUCT_TAC THEN REWRITE_TAC[mat_pow;MULT_0] THEN 
REWRITE_TAC[mat_pow;ADD1;LEFT_ADD_DISTRIB;MULT_CLAUSES;MAT_POW_ADD] THEN 
POP_ASSUM MP_TAC THEN SIMP_TAC[EQ_SYM_EQ]);;

let MAT2X2_POW2_POW_GEN = prove 
(` !n (M:real^k^k). (M pow 2) pow n = M pow (2*n)`,
INDUCT_TAC THEN REWRITE_TAC[mat_pow;MULT_0] THEN 
REWRITE_TAC[mat_pow;ADD1;LEFT_ADD_DISTRIB;MULT_CLAUSES;MAT_POW_ADD] THEN 
POP_ASSUM MP_TAC THEN SIMP_TAC[EQ_SYM_EQ]);;

let MAT_POW2 = prove (`
  !(M:real^a^a). M pow 2 = M**M`,
REWRITE_TAC[mat_pow;TWO;ONE] THEN 
REWRITE_TAC[GSYM ONE; MATRIX_MUL_RID]);;

let TRACE_MAT2X2 = prove( `!A B C D. trace (mat2x2 A B C D) = A+D`,
 SIMP_TAC[SUM_2;trace;mat2x2;CART_EQ;DIMINDEX_2;FORALL_2;VECTOR_2]);;

let DET_MAT2X2 = prove( `!A B C D. det (mat2x2 A B C D) = A*D - B*C`,
 SIMP_TAC[SUM_2;DET_2;mat2x2;CART_EQ;DIMINDEX_2;FORALL_2;VECTOR_2]);;


(* ------------------------------------------------------------------------- *)
(* Matrix formalism                                                          *)
(* ------------------------------------------------------------------------- *)

let free_space_matrix = new_definition
  `free_space_matrix (n,d) = mat2x2 (&1) d (&0) (&1)`;;

let interface_matrix = define 
  `interface_matrix n0 n1 plane transmitted = mat2x2 (&1) (&0) (&0) (n0/n1) /\
  interface_matrix n0 n1 (spherical R) transmitted 
    = mat2x2 (&1) (&0) ((n0-n1)/n1 * &1/R) (n0/n1)/\
  interface_matrix n0 n1 plane reflected = mat2x2 (&1) (&0) (&0) (&1) /\
  interface_matrix n0 n1 (spherical R) reflected
    = mat2x2 (&1) (&0) (-- &2/R) (&1)`;;


(*-------------------------------------------------------------------------
---Verification of matrix relationships for free space and interfaces------
-- ** PI= Plane Interface *** SI = Spherical Interface ** FS = Free Space--
--------------------------------------------------------------------------*)

let common_prove t = prove (t,
  SIMP_TAC[free_space_matrix;FORALL_FREE_SPACE_THM;MAT2X2_VECTOR_MUL;VECTOR2_EQ;
    LET_END_DEF;is_valid_ray_in_free_space;REAL_MUL_LZERO;REAL_MUL_LID;LET_DEF;
    REAL_ADD_LID;interface_matrix;is_valid_ray_at_interface;REAL_ADD_RID]
  THEN CONV_TAC REAL_FIELD);;
 
let FS_MATRIX = common_prove
  `!fs y0 theta0 y1 theta1.
    is_valid_free_space fs
    /\ is_valid_ray_in_free_space (y0,theta0) (y1,theta1) fs
    ==> vector[y1;theta1] = free_space_matrix fs ** vector[y0;theta0]`;;

let PI_MATRIX_TRANSMITTED = common_prove
  `!n0 n1 y0 theta0 y1 theta1.
    is_valid_ray_at_interface (y0,theta0) (y1,theta1) n0 n1 plane transmitted
    /\ &0 < n0 /\ &0 < n1
      ==> vector [y1;theta1]
        = interface_matrix n0 n1 plane transmitted ** vector [y0; theta0]`;;

let PI_MATRIX_REFLECTED = common_prove
  `!n0 y0 theta0 y1 theta1.
    is_valid_ray_at_interface (y0,theta0) (y1,theta1) n0 n1 plane reflected
    /\ &0 < n0 /\ &0 < n1
      ==> vector [y1;theta1]
        = interface_matrix n0 n1 plane reflected ** vector [y0;theta0]`;;

let PI_MATRIX = prove
  (`!ik n0 n1 y0 theta0 y1 theta1.
     is_valid_ray_at_interface (y0,theta0) (y1,theta1) n0 n1 plane ik
     /\ &0 < n0 /\ &0 < n1
       ==> vector [y1;theta1]
        = interface_matrix n0 n1 plane ik ** vector [y0;theta0]`,
  MATCH_MP_TAC interface_kind_IND
  THEN SIMP_TAC[PI_MATRIX_TRANSMITTED;PI_MATRIX_REFLECTED]);;

let SI_MATRIX_TRANSMITTED = common_prove
  `!n0 n1 y0 theta0 y1 theta1 R.
   is_valid_interface (spherical R) /\ &0 < n0 /\  &0 < n1
   /\ is_valid_ray_at_interface (y0,theta0) (y1,theta1) n0 n1 (spherical R)
    transmitted
     ==> vector [y1;theta1]
      = interface_matrix n0 n1 (spherical R) transmitted ** vector [y0;theta0]`;;

let SI_MATRIX_REFLECTED = common_prove
  `!n0 y0 theta0 y1 theta1 R.
   is_valid_interface (spherical R) /\ &0 < n0
   /\ is_valid_ray_at_interface (y0,theta0) (y1,theta1) n0 n1 (spherical R)
    reflected
     ==> vector [y1;theta1]
      = interface_matrix n0 n1 (spherical R) reflected ** vector [y0; theta0]`;;

let SI_MATRIX = prove
 (`!ik n0 n1 y0 theta0 y1 theta1.
  is_valid_interface (spherical R)
  /\ is_valid_ray_at_interface (y0,theta0) (y1,theta1) n0 n1 (spherical R) ik
  /\ &0 < n0 /\ &0 < n1
   ==> vector [y1; theta1] = interface_matrix n0 n1 (spherical R) ik ** vector [y0; theta0]`,
  MATCH_MP_TAC interface_kind_IND
  THEN SIMP_TAC[SI_MATRIX_TRANSMITTED;SI_MATRIX_REFLECTED ]);;


(*-------General Relationship for any Interface i --------------------*)


let INTERFACE_MATRIX = prove
  (`!n0 n1 y0 theta0 y1 theta1 i ik.
  is_valid_interface i /\ &0 < n0 /\ &0 < n1 /\
  is_valid_ray_at_interface (y0,theta0) (y1,theta1) n0 n1 i ik
    ==> vector [y1;theta1] = interface_matrix n0 n1 i ik ** vector [y0;theta0]`,
  REPEAT (MATCH_MP_TAC interface_IND ORELSE GEN_TAC)
  THEN MESON_TAC[PI_MATRIX;SI_MATRIX]);;


(*-----------------------------------------------*)

let system_composition = define 
  `system_composition ([],fs) = free_space_matrix fs /\ 
   system_composition (CONS ((n',d'),i,ik) cs,n,d) =
     system_composition (cs,n,d)
     ** interface_matrix n' (head_index (cs,n,d)) i ik
     ** free_space_matrix (n',d')`;;

let IS_VALID_OPTICAL_SYSTEM_REC = prove
  (`!fs fs' i ik t.
    is_valid_optical_system ([],fs) = is_valid_free_space fs /\
    (is_valid_optical_system (CONS (fs,i,ik) t,fs') <=> is_valid_free_space fs
    /\ is_valid_interface i /\ is_valid_optical_system (t,fs'))`,
  REWRITE_TAC[is_valid_optical_system;ALL;is_valid_optical_component]
  THEN MESON_TAC[]);;

let LAST_SINGLE_RAY_REC = prove
  (`!sr1 sr2 sr3 sr4 t.
    last_single_ray (sr1,sr2,[]) = sr2 /\
    last_single_ray (sr1,sr2,CONS (sr3,sr4) t) = last_single_ray (sr3,sr4,t)`,
  REPEAT (LIST_INDUCT_TAC ORELSE GEN_TAC)
  THEN SIMP_TAC[last_single_ray;LAST;NOT_CONS_NIL]);;

let IS_VALID_OPTICAL_SYSTEM_HEAD_INDEX = prove
  (`!os:optical_system. is_valid_optical_system os ==> &0 < head_index os`,
  REWRITE_TAC[FORALL_OPTICAL_SYSTEM_THM] THEN MATCH_MP_TAC list_INDUCT
  THEN SIMP_TAC[FORALL_FREE_SPACE_THM;IS_VALID_OPTICAL_SYSTEM_REC;ALL;
    is_valid_free_space;head_index;FORALL_OPTICAL_COMPONENT_THM]);;

let SYSTEM_MATRIX = prove
  (`!sys ray.
    is_valid_optical_system sys /\ is_valid_ray_in_system ray sys ==> 
      let ((y0,theta0),(y1,theta1),rs) = ray in
      let (y_n,theta_n) = last_single_ray ray in
        vector [y_n;theta_n] = system_composition sys ** vector [y0;theta0]`,
  REWRITE_TAC[FORALL_OPTICAL_SYSTEM_THM;FORALL_RAY_THM]
  THEN MATCH_MP_TAC list_INDUCT THEN CONJ_TAC
  THEN REWRITE_TAC[FORALL_SINGLE_RAY_THM;FORALL_OPTICAL_COMPONENT_THM] 
  THENL [ALL_TAC;REPEAT GEN_TAC THEN DISCH_TAC]
  THEN REPEAT (MATCH_MP_TAC list_INDUCT ORELSE GEN_TAC) THEN CONJ_TAC
  THEN REWRITE_TAC[IS_VALID_OPTICAL_SYSTEM_REC;is_valid_ray_in_system;ALL;
    system_composition;head_index]
  THENL [
    REWRITE_TAC[last_single_ray;LET_DEF;LET_END_DEF;FS_MATRIX];
    let lemma = prove(`!P. (!a. P a) <=> !y1 th1 y2 th2. P((y1,th1),y2,th2)`,
      REWRITE_TAC[FORALL_PAIR_THM]) in
    REWRITE_TAC[lemma;is_valid_ray_in_system;GSYM MATRIX_VECTOR_MUL_ASSOC;
      LAST_SINGLE_RAY_REC;is_valid_free_space]
    THEN REPEAT STRIP_TAC THEN REPEAT LET_TAC
    THEN ASSUM_LIST (REWRITE_TAC o map (REWRITE_RULE[PAIR_EQ] o GSYM))
    THEN REPEAT_GTCL IMP_RES_THEN (wrap REWRITE_TAC)
      (REWRITE_RULE[IMP_CONJ;is_valid_free_space;FORALL_FREE_SPACE_THM]
        (GSYM FS_MATRIX))
    THEN IMP_RES_THEN ASSUME_TAC IS_VALID_OPTICAL_SYSTEM_HEAD_INDEX
    THEN REPEAT_GTCL IMP_RES_THEN (wrap REWRITE_TAC)
      (REWRITE_RULE[IMP_CONJ] (GSYM INTERFACE_MATRIX))
    THEN EVERY_ASSUM (REPEAT_GTCL IMP_RES_THEN MP_TAC o REWRITE_RULE[IMP_CONJ])
    THEN ASM_REWRITE_TAC[LET_DEF;LET_END_DEF] THEN SIMP_TAC[]
  ]);;