(* ========================================================================= *)
(*               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>                     *)
(*                                                                           *)
(* Extentions of the vector library.                                         *)
(* ========================================================================= *)


prioritize_vector ();;

(** Additions to euclidean space *)

let [ORTHOGONAL_RZERO;ORTHOGONAL_RMUL;ORTHOGONAL_RNEG;ORTHOGONAL_RADD;
  ORTHOGONAL_RSUB;ORTHOGONAL_LZERO;ORTHOGONAL_LMUL;ORTHOGONAL_LNEG;
  ORTHOGONAL_LADD;ORTHOGONAL_LSUB] = 
  CONJUNCTS ORTHOGONAL_CLAUSES;;

let NON_VEC0 = 
prove
  (`!x:real^N. ~(x=vec 0) <=> ?i. 1 <= i /\ i <= dimindex(:N) /\ ~(x$i = &0)`,

  GEN_TAC THEN EQ_TAC THEN SIMP_TAC[vec;CART_EQ;LAMBDA_BETA]
  THEN ASM_MESON_TAC[]);;

let NON_VEC0_DIM3 = 
prove
  (`!x:real^3. ~(x=vec 0) <=> ~(x$1= &0) \/ ~(x$2= &0) \/ ~(x$3= &0)`,

  REWRITE_TAC[NON_VEC0;DIMINDEX_3]
  THEN MESON_TAC[ARITH_RULE `1 <= i /\ i <= 3 <=> i = 1 \/ i = 2 \/ i = 3`]);;

let ORTHOGONAL_SQR_NORM = 
prove
  (`!x y:real^N. orthogonal (x-y) y ==> x dot x = (x-y) dot (x-y) + y dot y`,

  REWRITE_TAC[orthogonal] THEN REPEAT STRIP_TAC
  THEN GEN_REWRITE_TAC (RATOR_CONV o ONCE_DEPTH_CONV) 
    [VECTOR_ARITH `x dot x = (x-y) dot (x-y) + y dot y + &2 * (x-y) dot y`]
  THEN ASM_REWRITE_TAC[REAL_MUL_RZERO;REAL_ADD_RID]);;

let COLLINEAR_VEC0 = 
prove
  (`!s:real^N->bool. collinear s /\ vec 0 IN s ==> ?u. !x. x IN s
    ==> ?c. x = c %u`,

  REWRITE_TAC[collinear] THEN REPEAT (STRIP_TAC ORELSE EXISTS_TAC `u:real^N`)
  THEN ASM_MESON_TAC[VECTOR_SUB_RZERO]);;

let CROSS_NORM_ORTHOGONAL = 
prove
  (`!x y. orthogonal x y ==> norm (x cross y) = norm x * norm y`,

  ONCE_REWRITE_TAC[GSYM (REWRITE_RULE[GSYM REAL_ABS_REFL] NORM_POS_LE)]
  THEN SIMP_TAC[GSYM REAL_ABS_MUL;REAL_EQ_SQUARE_ABS;GSYM NORM_CROSS_DOT;
    orthogonal;REAL_POW_2]
  THEN REAL_ARITH_TAC);;


let COLLINEAR_DEPENDENT = 
prove
  (`!x y:real^N. ~(x=y) /\ collinear {vec 0,x,y} ==> dependent {x,y}`,

  REWRITE_TAC[collinear;dependent;IN_INSERT;IN_SING] THEN REPEAT STRIP_TAC
  THEN ASM_CASES_TAC `x=vec 0:real^N`
  THENL ON_FIRST_GOAL (ASM_MESON_TAC[SPAN_0])
  THEN ASM_CASES_TAC `y=vec 0:real^N`
  THENL ON_FIRST_GOAL (ASM_MESON_TAC[SPAN_0])
  THEN SUBGOAL_THEN `?c. x = c % y :real^N` STRIP_ASSUME_TAC THENL [
    FIRST_ASSUM (fun x -> MAP_EVERY (STRIP_ASSUME_TAC o
      REWRITE_RULE[IN_INSERT;IN_SING;VECTOR_SUB_RZERO] o C SPECL x) [
        [`x:real^N`;`vec 0:real^N`]; [`y:real^N`;`vec 0:real^N`]])
    THEN ASM_REWRITE_TAC[VECTOR_MUL_ASSOC;VECTOR_MUL_RCANCEL]
    THEN SUBGOAL_THEN `~(c'= &0)` ASSUME_TAC
    THENL [ ASM_MESON_TAC[VECTOR_MUL_LZERO]; EXISTS_TAC `c/c':real`
    THEN ASM_SIMP_TAC[REAL_DIV_RMUL] ];
    SUBGOAL_TAC "" `~(y=c%y:real^N)` [ASM_MESON_TAC[]]
    THEN EXISTS_TAC `c%y:real^N` 
    THEN ASM_REWRITE_TAC[DELETE_INSERT;EMPTY_DELETE]
    THEN MESON_TAC[SUBSPACE_MUL;SUBSET;SPAN_INC;IN_SING;SUBSPACE_SPAN]
  ]);;

let ORTHOGONAL_DIFFERENT = 
prove
  (`!x y:real^N. orthogonal x y /\ ~(x=vec 0) /\ ~(y=vec 0) ==> ~(x=y)`,

  REWRITE_TAC[orthogonal] THEN MESON_TAC[DOT_EQ_0]);;

let ORTHOGONAL_NON_COLLINEAR = 
prove
  (`!x y:real^N. orthogonal x y /\ ~(x=vec 0) /\ ~(y=vec 0) 
    ==> ~(collinear {vec 0,x,y})`,

  REPEAT STRIP_TAC
  THEN SUBGOAL_THEN `independent {x,y:real^N}` ASSUME_TAC THENL [
    MATCH_MP_TAC PAIRWISE_ORTHOGONAL_INDEPENDENT
    THEN REWRITE_TAC[pairwise;IN_INSERT;IN_SING]
    THEN ASM_MESON_TAC[ORTHOGONAL_SYM];
    POP_ASSUM (MP_TAC o REWRITE_RULE[independent])
    THEN REWRITE_TAC[] THEN MATCH_MP_TAC COLLINEAR_DEPENDENT
    THEN ASM_SIMP_TAC[ORTHOGONAL_DIFFERENT]
  ]);;

(* NORMALIZING VECTORS *)

let unit = new_definition `unit v = inv(norm v) % v`;;

let UNIT_THM = 
prove
  (`!v. ~(v = vec 0) ==> norm(unit v) = &1`,

  SIMP_TAC[unit;NORM_MUL;REAL_ABS_INV;REAL_ABS_NORM;GSYM NORM_EQ_0;
    REAL_MUL_LINV]);;

let UNIT_INTRO = 
prove
  (`!v. ~(v = vec 0) ==> v = norm v % unit v`,

  SIMP_TAC[unit;VECTOR_MUL_ASSOC;GSYM NORM_EQ_0;REAL_MUL_RINV;
    VECTOR_MUL_LID]);;

let UNIT_ELIM = GSYM UNIT_INTRO;;

let UNIT_ZERO = 
prove
  (`unit(vec 0) = vec 0`,
 REWRITE_TAC[unit;VECTOR_MUL_RZERO]);;

let UNIT_EQ_ZERO = 
prove
  (`!v. unit v = vec 0 <=> v = vec 0`,

  REWRITE_TAC[unit;VECTOR_MUL_EQ_0;REAL_INV_EQ_0;NORM_EQ_0]);;

let UNIT_UNIT = 
prove
  (`!v. ~(v = vec 0) ==> unit (unit v) = unit v`,

  SIMP_TAC[unit;UNIT_THM;REAL_INV_1;VECTOR_MUL_LID]);;

let UNIT_DOT_UNIT_SELF = 
prove
  (`!v. ~(v = vec 0) ==> unit v dot unit v = &1`,

  SIMP_TAC[GSYM NORM_EQ_1;UNIT_THM]);;

let DOT_UNIT_SELF = 
prove
  (`!v. ~(v=vec 0) ==> v dot unit v = norm v`,

  SIMP_TAC[MESON[UNIT_INTRO] `!v. ~(v = vec 0) ==>
    v dot unit v = (norm v % unit v) dot unit v`;DOT_LMUL;UNIT_DOT_UNIT_SELF;
      REAL_MUL_RID]);;

let ORTHOGONAL_UNIT = 
prove
  (`!x y:real^N. (orthogonal x (unit y) <=> orthogonal x y)
    /\ (orthogonal (unit x) y <=> orthogonal x y)`,

  let common h =
    ASM_CASES_TAC h
    THENL [
      RULE_ASSUM_TAC (REWRITE_RULE[NORM_EQ_0]) THEN ASM_REWRITE_TAC[UNIT_ZERO];
      ASM_SIMP_TAC[unit;orthogonal;DOT_RMUL;DOT_LMUL;REAL_ENTIRE;REAL_INV_EQ_0]
    ]
  in
  REPEAT GEN_TAC THEN CONJ_TAC
  THENL [ common `norm (y:real^N) = &0`; common `norm (x:real^N) = &0` ]);;

let ORTHOGONAL_RUNIT = 
prove
  (`!x y:real^N. orthogonal x (unit y) <=> orthogonal x y`,

  MESON_TAC[ORTHOGONAL_UNIT]);;

let ORTHOGONAL_LUNIT = 
prove
  (`!x y:real^N. orthogonal (unit x) y <=> orthogonal x y`,

  MESON_TAC[ORTHOGONAL_UNIT]);;

let ORTHOGONAL_LRUNIT = 
prove
  (`!x y:real^N. orthogonal (unit x) (unit y) <=> orthogonal x y`,

  MESON_TAC[ORTHOGONAL_UNIT]);;

let COLLINEAR_UNIT = 
prove
  (`!x y:real^N. collinear {vec 0,x,y} <=> collinear {vec 0,unit x,unit y}`,

  REPEAT (STRIP_TAC ORELSE EQ_TAC) 
  THEN ASM_CSQ_THEN (REWRITE_RULE[GSYM IMP_IMP] COLLINEAR_VEC0) 
    (STRIP_ASSUME_TAC o REWRITE_RULE[IN_INSERT;IN_SING])
  THEN ASM_CSQ_THEN ~remove:true (MESON[] `(!x. x=a \/ x=b \/ x=c ==> p x) ==>
    p b /\ p c`) STRIP_ASSUME_TAC
  THENL [
    REWRITE_TAC[unit];
    MAP_EVERY ASM_CASES_TAC [`x=vec 0:real^N`;`y=vec 0:real^N`]
    THEN REPEAT (FIRST_X_ASSUM SUBST_VAR_TAC)
    THEN REWRITE_TAC[INSERT_AC;COLLINEAR_SING;COLLINEAR_2]
    THEN SUBGOAL_THEN `collinear {x:real^N, y, vec 0}
      <=> collinear {norm x % unit x, norm y % unit y, vec 0}` 
      (SINGLE REWRITE_TAC)
    THENL ON_FIRST_GOAL (ASM_SIMP_TAC[GSYM UNIT_INTRO])
  ]
  THEN REWRITE_TAC[collinear;IN_INSERT;IN_SING]
  THEN EXISTS_TAC `u:real^N` THEN REPEAT STRIP_TAC
  THEN ASM_REWRITE_TAC[VECTOR_MUL_ASSOC;VECTOR_SUB_RZERO;VECTOR_SUB_LZERO;
    VECTOR_NEG_0;GSYM VECTOR_MUL_LNEG;GSYM VECTOR_SUB_RDISTRIB]
  THEN MESON_TAC[VECTOR_MUL_LZERO]);;

let DOT_RUNIT_POS,DOT_LUNIT_POS = CONJ_PAIR(prove
  (`(!x y. &0 < x dot (unit y) <=> &0 < x dot y)
  /\ (!x y. &0 < (unit x) dot y <=> &0 < x dot y)`,
  REWRITE_TAC[unit;DOT_RMUL;DOT_LMUL;REAL_MUL_POS_LT;REAL_LT_INV_EQ;
    NORM_POS_LT;REWRITE_RULE[REAL_INV_NEG;REAL_NEG_GT0] 
    (SPEC `--x:real` REAL_LT_INV_EQ);
    REWRITE_RULE[REWRITE_RULE[MESON[]`(p<=> ~q) <=> (q<=> ~p)`] real_lt] 
    NORM_POS_LE]
  THEN MESON_TAC[DOT_RZERO;DOT_LZERO;REAL_LT_REFL]));;

let DOT_RUNIT_EQ0,DOT_LUNIT_EQ0 = CONJ_PAIR(prove
  (`(!x y. x dot (unit y) = &0 <=> x dot y = &0) 
  /\ (!x y. (unit x) dot y = &0 <=> x dot y = &0)`,
  REWRITE_TAC[GSYM orthogonal;ORTHOGONAL_LUNIT;ORTHOGONAL_RUNIT]));;

let DOT_RUNIT_LT0,DOT_LUNIT_LT0 = CONJ_PAIR(prove
  (`(!x y. x dot (unit y) < &0 <=> x dot y < &0) 
  /\ (!x y. (unit x) dot y < &0 <=> x dot y < &0)`,
  REWRITE_TAC[real_lt;MESON[] `(~p <=> ~q) <=> (p <=> q)`] 
  THEN REWRITE_TAC[REAL_LE_LT]
  THEN MESON_TAC[DOT_RUNIT_EQ0;DOT_LUNIT_EQ0;DOT_RUNIT_POS;DOT_LUNIT_POS]));;

let DOT_RUNIT_LE0,DOT_LUNIT_LE0 = CONJ_PAIR(prove
  (`(!x y. x dot (unit y) <= &0 <=> x dot y <= &0)
  /\ (!x y. (unit x) dot y <= &0 <=> x dot y <= &0)`,
  MESON_TAC[REAL_LE_LT;DOT_RUNIT_LT0;DOT_LUNIT_LT0;DOT_RUNIT_EQ0;
    DOT_LUNIT_EQ0]));;

let DOT_RUNIT_GE0,DOT_LUNIT_GE0 = CONJ_PAIR(prove
  (`(!x y. &0 <= x dot (unit y) <=> &0 <= x dot y)
  /\ (!x y. &0 <= (unit x) dot y <=> &0 <= x dot y)`,
  REWRITE_TAC[REAL_LE_LT]
  THEN MESON_TAC[DOT_RUNIT_POS;DOT_LUNIT_POS;DOT_RUNIT_EQ0;DOT_LUNIT_EQ0]));;

let DOT_RUNIT_EQ = 
prove
  (`!x y z:real^N. x dot (unit z) = y dot (unit z) <=> x dot z = y dot z`,

  REPEAT GEN_TAC THEN ASM_CASES_TAC `z=vec 0:real^N`
  THEN ASM_REWRITE_TAC[unit;DOT_RMUL;REAL_EQ_MUL_LCANCEL;REAL_INV_EQ_0;
    NORM_EQ_0;DOT_RZERO]);;

let DOT_LUNIT_EQ = 
prove
  (`!x y z:real^N. (unit z) dot x = (unit z) dot y <=> z dot x = z dot y`,

  REPEAT GEN_TAC THEN ASM_CASES_TAC `z=vec 0:real^N`
  THEN ASM_REWRITE_TAC[unit;DOT_LMUL;REAL_EQ_MUL_LCANCEL;REAL_INV_EQ_0;
    NORM_EQ_0;DOT_LZERO]);;


(* CONNECTION BETWEEN SPANS AND PLANES *)

new_type_abbrev("point",`:real^3`);;
new_type_abbrev("plane",`:point->bool`);;

let PLANE_NON_EMPTY = 
prove
  (`!p:plane. plane p ==> ~(p={})`,

  REWRITE_TAC[plane;GSYM MEMBER_NOT_EMPTY] THEN REPEAT STRIP_TAC 
  THEN EXISTS_TAC `u:point` 
  THEN ASM_REWRITE_TAC[] THEN MATCH_MP_TAC (REWRITE_RULE[SUBSET] HULL_SUBSET) 
  THEN SET_TAC[]);;

let is_plane_basis = new_definition
  `is_plane_basis s p <=>
    ?u v:real^3. s = {u,v} /\ ~(collinear {vec 0,u,v}) /\ !pt:point. pt IN p 
      ==> !pt':point. pt' IN p <=> ?a b:real. pt' = pt + a % u + b % v`;;

(* That theorem is (one of) the last which makes use of the affine definition of `plane` *)
let EXISTS_PLANE_BASIS = 
prove
  (`!p. plane p ==> ?b. is_plane_basis b p`,

  GEN_TAC
  THEN DISCH_THEN (fun x -> 
    STRIP_ASSUME_TAC (REWRITE_RULE[plane;AFFINE_HULL_3] x)
    THEN STRIP_ASSUME_TAC (MATCH_MP (REWRITE_RULE[GSYM MEMBER_NOT_EMPTY] 
      PLANE_NON_EMPTY) x))
  THEN POP_ASSUM MP_TAC THEN ASM_REWRITE_TAC[is_plane_basis] THEN DISCH_TAC
  THEN MAP_EVERY EXISTS_TAC [`{v-u,w-u:real^3}`;`v-u:real^3`; `w-u:real^3`]
  THEN ASM_REWRITE_TAC[] THEN REPEAT STRIP_TAC 
  THENL [
    ASM_MESON_TAC[INSERT_AC;COLLINEAR_3];
    POP_ASSUM MP_TAC THEN REWRITE_TAC[IN_ELIM_THM] THEN REPEAT STRIP_TAC
    THEN EQ_TAC THEN STRIP_TAC
    THENL [
      RULE_ASSUM_TAC (REWRITE_RULE[REAL_ARITH `x+y+z= &1 <=> x= &1-y-z`])
      THEN MAP_EVERY EXISTS_TAC [`v''-v':real`; `w''-w':real`] 
      THEN ASM_REWRITE_TAC[] THEN VECTOR_ARITH_TAC;
      MAP_EVERY EXISTS_TAC [`u'-a-b:real`;`v'+a:real`;`w'+b:real`] 
      THEN CONJ_TAC
      THENL [
        ASM_ARITH_TAC;
        REPLICATE_TAC 2 (POP_ASSUM MP_TAC) THEN VECTOR_ARITH_TAC]
    ]]);;


let BASIS_NON_NULL = 
prove
  (`!p. plane p ==> !b. is_plane_basis b p ==> !u. u IN b ==> ~(u=vec 0)`,

  REWRITE_TAC[is_plane_basis] THEN REPEAT STRIP_TAC
  THEN SUBGOAL_THEN `(u':real^3) IN {u,v}` MP_TAC
  THENL ON_FIRST_GOAL (ASM_MESON_TAC[])
  THEN REWRITE_TAC[IN_INSERT;IN_SING]
  THEN ASM_MESON_TAC[INSERT_AC;COLLINEAR_2]);;

let PAIR_SETS_EQ = 
prove
  (`!x y z t. {x,y} = {z,t} ==> (x=z /\ y=t) \/ (x=t /\ y=z)`,

  SET_TAC[]);;

let NON_COLLINEAR3_IMPLIES_DIFFERENT = 
prove
  (`!x y z. ~(collinear {x,y,z}) ==> ~(x=y) /\ ~(y=z) /\ ~(x=z)`,

  MESON_TAC[COLLINEAR_2;INSERT_AC]);;

let BASIS_DIFFERENT = 
prove
  (`!p. plane p ==> !u v. is_plane_basis {u,v} p ==> ~(u=v)`,

  REWRITE_TAC[is_plane_basis] THEN REPEAT STRIP_TAC
  THEN ASM_MESON_TAC[PAIR_SETS_EQ;NON_COLLINEAR3_IMPLIES_DIFFERENT]);;

let UNIT_OF_BASIS_IS_BASIS = 
prove
  (`!p. plane p ==> !u v. is_plane_basis {u,v} p ==> is_plane_basis {unit u,unit v} p`,

  REPEAT STRIP_TAC THEN REWRITE_TAC[is_plane_basis]
  THEN MAP_EVERY EXISTS_TAC [`unit u:real^3`; `unit v:real^3`]
  THEN REWRITE_TAC[]
  THEN POP_ASSUM 
    (DISTRIB [ASSUME_TAC;STRIP_ASSUME_TAC o REWRITE_RULE[is_plane_basis]])
  THEN ASM_REWRITE_TAC[GSYM COLLINEAR_UNIT]
  THEN ASM_CSQ_THEN ~remove:true PAIR_SETS_EQ STRIP_ASSUME_TAC
  THEN REPEAT (FIRST_X_ASSUM SUBST_VAR_TAC)
  THEN ASM_CSQ_THEN BASIS_NON_NULL (C ASM_CSQ_THEN (STRIP_ASSUME_TAC 
    o REWRITE_RULE[IN_INSERT;IN_SING;
      MESON[]`(!u. u=a \/ u=b ==> p u) <=> p a /\ p b`]))
  THEN GEN_TAC 
  THEN DISCH_THEN (fun x -> FIRST_X_ASSUM (ASSUME_TAC o C MATCH_MP x))
  THEN ASM_REWRITE_TAC[unit;VECTOR_MUL_ASSOC]
  THEN REPEAT (STRIP_TAC ORELSE EQ_TAC) THEN ASM_REWRITE_TAC[]
  THENL [
    MAP_EVERY EXISTS_TAC [`a*norm(u':real^3)`;`b*norm(v':real^3)`] 
    THEN RULE_ASSUM_TAC (REWRITE_RULE[GSYM NORM_EQ_0])
    THEN ASM_SIMP_TAC[GSYM REAL_MUL_ASSOC;REAL_MUL_RINV] THEN VECTOR_ARITH_TAC;
    MESON_TAC[];
    MAP_EVERY EXISTS_TAC [`b*norm(v':real^3)`;`a*norm(u':real^3)`] 
    THEN RULE_ASSUM_TAC (REWRITE_RULE[GSYM NORM_EQ_0])
    THEN ASM_SIMP_TAC[GSYM REAL_MUL_ASSOC;REAL_MUL_RINV] THEN VECTOR_ARITH_TAC;
    MESON_TAC[VECTOR_ADD_AC];
  ]);;

(* This theorem makes the connection between the (affine) notion of plane and the (vector) notion of span.
 * From now on, we can generally rely only on this result to reason about the spanning set of a plane.
 *)
let PLANE_AS_SPAN = 
prove
  (`!p. plane p ==> !b. is_plane_basis b p ==> !pt:point. pt IN p 
    ==> p = { pt + y | y IN span b }`,

  REWRITE_TAC[is_plane_basis] THEN REPEAT STRIP_TAC
  THEN ASM_REWRITE_TAC[EXTENSION;IN_ELIM_THM;SPAN_2;UNIV]
  THEN ASM_MESON_TAC[]);;

let PLANE_SPAN = 
prove
  (`!p. plane p ==> !b. is_plane_basis b p ==> !pt:point. pt IN p
    ==> { pt' - pt | pt' IN p } = span b`,

  REPEAT STRIP_TAC THEN ASM_CSQ_THEN PLANE_AS_SPAN (C ASM_CSQ_THEN (C
    ASM_CSQ_THEN (SINGLE (GEN_REWRITE_TAC (RATOR_CONV o ONCE_DEPTH_CONV)))))
  THEN REWRITE_TAC[IN_ELIM_THM;EXTENSION] THEN GEN_TAC THEN EQ_TAC
  THENL [
    REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[VECTOR_ARITH `(x + y)-x = y:real^N`];
    DISCH_TAC THEN EXISTS_TAC `x+pt:real^3` THEN REWRITE_TAC (map VECTOR_ARITH
      [`(x + y) - y = x:real^N`; `x+z=z+y <=> x=y:real^N`]) THEN ASM_MESON_TAC[]
  ]);;

let ALL_BASIS_SAME_SPAN = 
prove
  (`!p. plane p ==> !b1. is_plane_basis b1 p ==> !b2. is_plane_basis b2 p
    ==> span b1 = span b2`,

  REPEAT STRIP_TAC
  THEN ASM_CSQ_THEN PLANE_SPAN (fun x ->
    ASM_CSQ_THEN ~match_:false (SPEC `b1:real^3->bool` x) ASSUME_TAC
      THEN ASM_CSQ_THEN (SPEC `b2:real^3->bool` x) ASSUME_TAC)
  THEN ASM_MESON_TAC[PLANE_NON_EMPTY;MEMBER_NOT_EMPTY]);;

let PLANE_SUBSPACE = 
prove
  (`!p. plane p ==> !pt:point. pt IN p ==> subspace { pt' - pt | pt' IN p }`,

  REPEAT STRIP_TAC THEN ASM_CSQ_THEN EXISTS_PLANE_BASIS STRIP_ASSUME_TAC
  THEN ASM_SIMP_TAC[PLANE_SPAN;SUBSPACE_SPAN]);;

let PLANE_SING = 
prove
  (`!x:real^N. ~(plane{x})`,

  REWRITE_TAC[plane;NOT_EXISTS_THM;MESON[] `~(A/\B) <=> (A ==> ~B)`]
  THEN REPEAT STRIP_TAC 
  THEN ASM_CSQ_THEN NON_COLLINEAR3_IMPLIES_DIFFERENT ASSUME_TAC
  THEN MATCH_MP_TAC (SET_RULE `~(u=v:real^N) /\ u IN {x} /\ v IN {x} ==> F`) 
  THEN ASM_REWRITE_TAC[]
  THEN CONJ_TAC THEN MATCH_MP_TAC HULL_INC THEN SET_TAC[]);;

let NON_COLLINEAR_INDEPENDENT = 
prove
  (`!x y:real^N. ~(collinear {vec 0,x,y}) ==> independent{x,y}`,

  MAP_EVERY (SINGLE ONCE_REWRITE_TAC) [independent;CONTRAPOS_THM;NOT_CLAUSES]
  THEN REPEAT GEN_TAC
  THEN ASM_CASES_TAC `x=y:real^N` THENL [
    ASM_REWRITE_TAC[INSERT_AC;COLLINEAR_2];
    REWRITE_TAC[dependent;IN_INSERT;IN_SING] THEN STRIP_TAC THENL
    let common main_var exist_list =
      POP_ASSUM MP_TAC
      THEN ASM_REWRITE_TAC[DELETE_INSERT;EMPTY_DELETE;SPAN_SING;IN_ELIM_THM;
        UNIV]
      THEN STRIP_TAC
      THEN REWRITE_TAC[collinear;IN_INSERT;IN_SING] THEN EXISTS_TAC main_var
      THEN REPEAT STRIP_TAC
      THENL map (fun t ->
        EXISTS_TAC t THEN ASM_REWRITE_TAC[] THEN VECTOR_ARITH_TAC) exist_list
    in
    [ common `y:real^N`
        [`&0`;`--u:real`;`-- &1`;`u:real`;`&0`;`u- &1`;`&1`;`&1-u`;`&0`];
      common `x:real^N`
        [`&0`;`-- &1`;`--u:real`;`&1`;`&0`;`&1-u`;`u:real`;`u- &1`;`&0`]
    ]]);;

let DIM_OF_PLANE_SPAN = 
prove
  (`!p:plane. plane p ==> !b. is_plane_basis b p ==> dim (span b) = 2`,

  REWRITE_TAC[is_plane_basis] THEN REPEAT STRIP_TAC 
  THEN ASM_CSQ_THEN NON_COLLINEAR_INDEPENDENT ASSUME_TAC
  THEN ASM_CSQ_THEN NON_COLLINEAR3_IMPLIES_DIFFERENT ASSUME_TAC
  THEN SUBGOAL_THEN `{u,v:real^3} HAS_SIZE (dim (span {u,v}))` MP_TAC
  THENL [
    MATCH_MP_TAC BASIS_HAS_SIZE_DIM THEN ASM_REWRITE_TAC[];
    ASM_SIMP_TAC[HAS_SIZE;FINITE_INSERT;FINITE_EMPTY;CARD_CLAUSES;IN_INSERT;
    IN_SING;NOT_IN_EMPTY] THEN ARITH_TAC
  ]);;

let DIM_OF_PLANE_SUBSPACE = 
prove
  (`!p. plane p ==> !pt:point. pt IN p ==> dim { pt' - pt | pt' IN p } = 2`,

  MESON_TAC[PLANE_SPAN;DIM_OF_PLANE_SPAN;EXISTS_PLANE_BASIS]);;

let PLANE_SPAN_DECOMPOS = 
prove
  (`!p. plane p ==> !u v. is_plane_basis {u,v} p ==>
    !x. x IN span {u,v} <=> ?a b. x = a % u + b % v`,

  REPEAT STRIP_TAC THEN POP_ASSUM (SINGLE REWRITE_TAC o GSYM)
  THEN REWRITE_TAC[SPAN_2;IN_ELIM_THM;UNIV]);;

let units = new_definition `units = { x | norm x = &1 }`;;

let IN_UNITS = 
prove
  (`!x. x IN units ==> norm x = &1`,

  REWRITE_TAC[units;IN_ELIM_THM]);;

let is_normal_to_plane = new_definition
  `is_normal_to_plane (n:real^3) (p:plane) <=> ~(n=vec 0)
  /\ !b. is_plane_basis b p ==> !v. v IN (span b) ==> orthogonal v n`;;

let IS_NORMAL_TO_PLANE_UNIT = 
prove
  (`!p n. is_normal_to_plane (unit n) p <=> is_normal_to_plane n p`,

  REWRITE_TAC[is_normal_to_plane;ORTHOGONAL_RUNIT;UNIT_EQ_ZERO]);;

let EXISTS_NORMAL_OF_PLANE = 
prove
  (`!p:plane. plane p ==> ?n:real^3. is_normal_to_plane n p`,

  REWRITE_TAC[is_normal_to_plane] THEN REPEAT STRIP_TAC
  THEN ASM_CSQ_THEN EXISTS_PLANE_BASIS (CHOOSE_THEN (fun x ->
    MAP_EVERY STRIP_ASSUME_TAC [x;REWRITE_RULE[is_plane_basis] x]))
  THEN EXISTS_TAC `u cross v` THEN ASM_REWRITE_TAC[CROSS_EQ_0] 
  THEN REPEAT STRIP_TAC THEN POP_ASSUM MP_TAC
  THEN ASM_CSQ_THEN ALL_BASIS_SAME_SPAN (C (ASM_CSQ_THEN ~remove:true) (C
    ASM_CSQ_THEN (SINGLE REWRITE_TAC))) 
  THEN ASM_REWRITE_TAC[SPAN_2;IN_ELIM_THM;UNIV]
  THEN REPEAT STRIP_TAC
  THEN ASM_REWRITE_TAC[orthogonal;DOT_LADD;DOT_LMUL;DOT_CROSS_SELF]
  THEN SIMPLE_COMPLEX_ARITH_TAC);;

let NORMAL_OF_PLANE_IS_NORMAL_TO_BASIS = 
prove
  (`!p. plane p ==> !b. is_plane_basis b p ==> !u. u IN b ==>
    !n. is_normal_to_plane n p ==> orthogonal u n`,

  MESON_TAC[is_normal_to_plane;SPAN_SUPERSET;]);;

let NORMAL_OF_PLANE_NON_NULL = 
prove
  (`!p. plane p ==> !n. is_normal_to_plane n p ==> ~(n=vec 0)`,

  MESON_TAC[is_normal_to_plane]);;

let NORMAL_OF_PLANE_IS_ORTHOGONAL_TO_SEGMENT = 
prove
  (`!p. plane p ==> !n. is_normal_to_plane n p 
    ==> !pt1 pt2. pt1 IN p /\ pt2 IN p ==> orthogonal (pt1-pt2) n`,

  REWRITE_TAC[is_normal_to_plane] THEN REPEAT STRIP_TAC 
  THEN ASM_CSQ_THEN EXISTS_PLANE_BASIS STRIP_ASSUME_TAC
  THEN POP_ASSUM
    (fun x -> FIRST_X_ASSUM (MATCH_MP_TAC o C MATCH_MP x) THEN ASSUME_TAC x)
  THEN FIRST_X_ASSUM (STRIP_ASSUME_TAC o CONV_RULE (REWR_CONV is_plane_basis))
  THEN POP_ASSUM (C (ASM_CSQ_THEN ~remove:true) (fun x -> 
    FIRST_X_ASSUM (STRIP_ASSUME_TAC o CONV_RULE (REWR_CONV x))))
  THEN ASM_REWRITE_TAC[SPAN_2;UNIV;IN_ELIM_THM;VECTOR_ARITH
    `(x+y+z)-x = y+z:real^N`]
  THEN MESON_TAC[]);;

let DIM_OF_NORMAL_SPAN = 
prove
  (`!p:plane. plane p ==> !n. is_normal_to_plane n p ==> dim (span {n}) = 1`,

  REWRITE_TAC[is_normal_to_plane] THEN REPEAT STRIP_TAC 
  THEN SUBGOAL_THEN `{n:real^3} HAS_SIZE dim (span {n})` MP_TAC 
  THENL [
    MATCH_MP_TAC BASIS_HAS_SIZE_DIM 
    THEN REWRITE_TAC[independent;dependent;NOT_EXISTS_THM;
      MESON[] `~(p/\q) <=> p ==> ~q`;IN_SING]
    THEN GEN_TAC THEN DISCH_TAC
    THEN ASM_REWRITE_TAC[DELETE_INSERT;EMPTY_DELETE;SPAN_EMPTY;IN_SING];
    ASM_SIMP_TAC[HAS_SIZE;FINITE_INSERT;FINITE_EMPTY;CARD_CLAUSES;IN_INSERT;
      IN_SING;NOT_IN_EMPTY]
    THEN ARITH_TAC]);;

let DIM_OF_ORTHOGONAL = 
prove
  (`!p. plane p ==> !b. is_plane_basis b p 
    ==> dim {z:real^3 | !x. x IN span b ==> orthogonal z x} = 1`,

  REPEAT STRIP_TAC
  THEN SUBGOAL_THEN
    `(:real^3) = { x+y 
      | x IN span b /\ y IN {z:real^3 | !x. x IN span b ==> orthogonal z x}}`
  (ASSUME_TAC o REWRITE_RULE[DIM_UNIV;DIMINDEX_3] o 
    AP_TERM `dim:(real^3->bool)->num`) 
  THENL ON_FIRST_GOAL 
    (REWRITE_TAC[GSYM SUBSET_ANTISYM_EQ;SUBSET;UNIV;IN;IN_ELIM_THM] 
    THEN MESON_TAC[ORTHOGONAL_SUBSPACE_DECOMP_EXISTS;IN])
  THEN SUBGOAL_THEN
    `span b
      INTER {z:real^3 | !x. x IN span b ==> orthogonal z x} SUBSET {vec 0}` 
    (ASSUME_TAC o GSYM o REWRITE_RULE[GSYM DIM_EQ_0])
  THENL ON_FIRST_GOAL 
    (REWRITE_TAC[INTER;SUBSET;IN_SING;IN_ELIM_THM]
    THEN MESON_TAC[ORTHOGONAL_REFL])
  THEN ASM_CSQ_THEN DIM_OF_PLANE_SPAN (fun x ->
    ASM_SIMP_TAC[ARITH_RULE `x=1 <=> 3+0=2+x`;GSYM x])
  THEN MATCH_MP_TAC DIM_SUMS_INTER
  THEN REWRITE_TAC[SUBSPACE_ORTHOGONAL_TO_VECTORS;ORTHOGONAL_SYM;
    SUBSPACE_SPAN]
  );;

let NORMAL_SPAN_SUBSET_ORTHOGONAL_SUBSPACE = 
prove
  (`!p. plane p ==> !b. is_plane_basis b p ==> !n. is_normal_to_plane n p 
    ==> span {n} SUBSET {z:real^3 | !x. x IN span b ==> orthogonal z x}`,

  REWRITE_TAC[is_normal_to_plane] THEN REPEAT STRIP_TAC 
  THEN POP_ASSUM (C ASM_CSQ_THEN ASSUME_TAC)
  THEN REWRITE_TAC[SUBSET;SPAN_SING;IN_ELIM_THM;UNIV] 
  THEN ASM_MESON_TAC[ORTHOGONAL_CLAUSES;ORTHOGONAL_SYM]);;

let ORTHOGONAL_SUBSPACE_IS_NORMAL_SPAN = 
prove
  (`!p. plane p ==> !b. is_plane_basis b p ==> !n. is_normal_to_plane n p 
    ==> {z:real^3 | !x. x IN span b ==> orthogonal z x} = span {n}`,

  MAP_EVERY (SINGLE ONCE_REWRITE_TAC) [
    ORTHOGONAL_SYM; MESON[SPAN_SPAN] `span{x}=span(span{x})`; 
    GSYM (REWRITE_RULE[GSYM SPAN_EQ_SELF] SUBSPACE_ORTHOGONAL_TO_VECTORS)]
  THEN REPEAT STRIP_TAC THEN MATCH_MP_TAC EQ_SYM 
  THEN MATCH_MP_TAC DIM_EQ_SPAN THEN ONCE_REWRITE_TAC[ORTHOGONAL_SYM]
  THEN CONJ_TAC
  THENL ON_FIRST_GOAL (ASM_MESON_TAC[NORMAL_SPAN_SUBSET_ORTHOGONAL_SUBSPACE])
  THEN MAP_EVERY (fun x ->
    ASM_CSQ_THEN x (C ASM_CSQ_THEN (SINGLE REWRITE_TAC))) 
    [DIM_OF_ORTHOGONAL;DIM_OF_NORMAL_SPAN]
  THEN ARITH_TAC);;

let PLANE_DECOMP = 
prove
  (`!p. plane p ==> !b. is_plane_basis b p ==> !n. is_normal_to_plane n p 
    ==> !v. ?w x. v = w + x % n /\ w IN span b`,

  REPEAT STRIP_TAC 
  THEN MP_TAC (ISPEC `b:real^3->bool` ORTHOGONAL_SUBSPACE_DECOMP)
  THEN ASM_CSQ_THEN ORTHOGONAL_SUBSPACE_IS_NORMAL_SPAN 
    (C ASM_CSQ_THEN (C ASM_CSQ_THEN (SINGLE REWRITE_TAC)))
  THEN REWRITE_TAC[SPAN_SING;IN_ELIM_THM;EXISTS_UNIQUE_DEF;EXISTS_PAIRED_THM]
  THEN MESON_TAC[]);;

let NORMAL_ORTHOGONAL_IN_PLANE_SPAN = 
prove
  (`!p. plane p ==> !b. is_plane_basis b p ==> !n. is_normal_to_plane n p 
    ==> !v. orthogonal v n <=> v IN span b`,

  let EQ_IMPLY =
    MATCH_MP (MESON[] `(!x y. p x y <=> q x y) ==> (!x y. p x y ==> q x y)`)
  in
  REPEAT STRIP_TAC
  THEN ASM_CSQ_THEN PLANE_DECOMP (C ASM_CSQ_THEN (C ASM_CSQ_THEN 
    (STRIP_ASSUME_TAC o SPEC `v:real^3`)))
  THEN ASM_CSQ_THEN (EQ_IMPLY is_normal_to_plane) STRIP_ASSUME_TAC 
  THEN POP_ASSUM (C ASM_CSQ_THEN ASSUME_TAC) THEN EQ_TAC THEN ASM_REWRITE_TAC[]
  THEN ASM_CSQ_THEN (EQ_IMPLY is_plane_basis) STRIP_ASSUME_TAC 
  THEN RULE_ASSUM_TAC (REWRITE_RULE[orthogonal])
  THEN ASM_SIMP_TAC[orthogonal;DOT_LADD;DOT_LMUL;REAL_ADD_LID;REAL_ENTIRE;
    DOT_EQ_0;VECTOR_MUL_LZERO;VECTOR_ADD_RID] THEN ASM_MESON_TAC[]);;

let PLANE_DECOMP_DOT_NORMAL = 
prove
  (`!p. plane p ==> !b. is_plane_basis b p ==> !n. is_normal_to_plane n p
    ==> !v. ?w. w IN span b /\ v = w + (v dot (unit n)) % (unit n)`,

  REPEAT STRIP_TAC
  THEN FIRST_ASSUM (ASSUME_TAC o CONJUNCT1 o CONV_RULE (REWR_CONV 
    is_normal_to_plane))
  THEN ASM_CSQ_THEN (GSYM NORMAL_ORTHOGONAL_IN_PLANE_SPAN) (C ASM_CSQ_THEN 
    (C ASM_CSQ_THEN ASSUME_TAC))
  THEN EXISTS_TAC `v - (v dot unit n) % unit n :real^3` 
  THEN ASM_REWRITE_TAC[VECTOR_ARITH `x=x-y+z <=> y=z:real^N`] 
  THEN ONCE_REWRITE_TAC[GSYM ORTHOGONAL_RUNIT]
  THEN REWRITE_TAC[orthogonal;DOT_LSUB;DOT_LMUL] 
  THEN ASM_MESON_TAC[UNIT_DOT_UNIT_SELF;REAL_MUL_RID;REAL_SUB_0]);;

let SUB_UNIT_NORMAL_IS_ORTHOGONAL_TO_NORMAL = 
prove
  (`!p. plane p ==> !n. is_normal_to_plane n p 
    ==> !x. orthogonal (x - (x dot unit n) % unit n) n`,

  REPEAT STRIP_TAC THEN ASM_CSQ_THEN EXISTS_PLANE_BASIS STRIP_ASSUME_TAC
  THEN ASM_CSQ_THEN PLANE_DECOMP_DOT_NORMAL 
    (C ASM_CSQ_THEN (C ASM_CSQ_THEN ASSUME_TAC))
  THEN RULE_ASSUM_TAC (REWRITE_RULE[VECTOR_ARITH `x=y+z <=> x-z=y`]) 
  THEN ASM_MESON_TAC[NORMAL_ORTHOGONAL_IN_PLANE_SPAN]);;

let is_orthogonal_plane_basis = new_definition
  `is_orthogonal_plane_basis b p 
    <=> is_plane_basis b p /\ pairwise orthogonal b`;;

let EXISTS_ORTHOGONAL_PLANE_BASIS = 
prove
  (`!p. plane p ==> ?b. is_orthogonal_plane_basis b p`,

  REWRITE_TAC[is_orthogonal_plane_basis;is_plane_basis] THEN REPEAT STRIP_TAC
  THEN ASM_CSQ_THEN EXISTS_PLANE_BASIS (CHOOSE_THEN (fun x -> ASSUME_TAC x
    THEN STRIP_ASSUME_TAC (REWRITE_RULE[is_plane_basis] x)))
  THEN ASM_CSQ_THEN EXISTS_NORMAL_OF_PLANE (CHOOSE_THEN (fun x ->
    let x' = REWRITE_RULE[is_normal_to_plane] x in
    ASSUME_TAC x THEN ASM_CSQ_THEN (CONJUNCT2 x') ASSUME_TAC
    THEN ASSUME_TAC (CONJUNCT1 x')))
  THEN SUBGOAL_TAC "" `~((u:real^3)=vec 0)` [ASM SET_TAC[BASIS_NON_NULL]]
  THEN SUBGOAL_TAC "" `orthogonal u (n:real^3)` [
    FIRST_ASSUM MATCH_MP_TAC THEN MATCH_MP_TAC SPAN_SUPERSET 
    THEN ASM SET_TAC [SPAN_SUPERSET]
  ]
  THEN SUBGOAL_TAC "" `~(u cross n = vec 0)` [
    ASM_SIMP_TAC[CROSS_EQ_0;ORTHOGONAL_NON_COLLINEAR] ]
  THEN EXISTS_TAC `{u, u cross n}` THEN CONJ_TAC 
  THENL [
    MAP_EVERY EXISTS_TAC [`u:real^3`;`u cross n`] THEN REWRITE_TAC[] 
    THEN CONJ_TAC THENL ON_FIRST_GOAL 
      (MATCH_MP_TAC ORTHOGONAL_NON_COLLINEAR 
      THEN ASM_REWRITE_TAC[orthogonal;DOT_CROSS_SELF])
    THEN SUBGOAL_THEN `?a b. u cross n = a % u + b % v` STRIP_ASSUME_TAC 
    THENL ON_FIRST_GOAL
      (ASM_MESON_TAC[GSYM PLANE_SPAN_DECOMPOS;
        GSYM NORMAL_ORTHOGONAL_IN_PLANE_SPAN;orthogonal;DOT_CROSS_SELF])
    THEN SUBGOAL_THEN `~(b'= &0)` (fun x -> POP_ASSUM MP_TAC THEN ASSUME_TAC x) 
    THENL ON_FIRST_GOAL 
      (ASM_CASES_TAC `a= &0` THENL [
        ASM_MESON_TAC[VECTOR_MUL_LZERO;VECTOR_ADD_LID];
        DISCH_THEN (fun x -> POP_ASSUM (fun y -> POP_ASSUM (MP_TAC 
          o REWRITE_RULE[x;VECTOR_MUL_LZERO;VECTOR_ADD_RID]) 
          THEN ASSUME_TAC y))
        THEN DISCH_THEN (MP_TAC o AP_TERM `(dot) (u:real^3)`) 
        THEN REWRITE_TAC[DOT_CROSS_SELF;DOT_RMUL] 
        THEN ONCE_REWRITE_TAC[EQ_SYM_EQ]
        THEN ASM_REWRITE_TAC[REAL_ENTIRE;DOT_EQ_0]
      ])
    THEN DISCH_THEN (fun x ->
      ASM_SIMP_TAC[x;VECTOR_ARITH `a%x+b%(c%x+d%y)=(a+b*c)%x+(b*d)%y`])
    THEN REPEAT STRIP_TAC THEN EQ_TAC THENL [
      REPEAT STRIP_TAC 
      THEN MAP_EVERY EXISTS_TAC [`a'-b''/b'*a:real`;`b''/b':real`]
      THEN ASM_SIMP_TAC[REAL_SUB_ADD;VECTOR_ARITH `x+y=x+z <=> y=z`;
        REAL_DIV_RMUL];
      MESON_TAC[];
    ];
    REWRITE_TAC[pairwise;orthogonal] THEN SET_TAC [DOT_CROSS_SELF]
  ]);;

let ORTHOGONAL_PLANE_BASIS_AS_PAIR = 
prove
  (`!b p. is_orthogonal_plane_basis b p ==> ?u v. is_orthogonal_plane_basis {u,v} p`,

  MESON_TAC[is_orthogonal_plane_basis;is_plane_basis]);;

let BASIS_OF_PLANE_ORTHONORMAL_DECOMPOS = 
prove
  (`!p. plane p ==> !u v. is_orthogonal_plane_basis {u,v} p ==>
  !x. x IN span {u,v}
    <=> x = (x dot unit u) % unit u + (x dot unit v) % unit v`,

  REWRITE_TAC [is_orthogonal_plane_basis] THEN REPEAT STRIP_TAC
  THEN ASM_CSQ_THEN PLANE_SPAN_DECOMPOS (C ASM_CSQ_THEN (SINGLE REWRITE_TAC))
  THEN EQ_TAC THENL [
    REPEAT STRIP_TAC
    THEN SUBGOAL_TAC "" `~(u = vec 0:real^3) /\ ~(v = vec 0:real^3)` 
      [ASM_MESON_TAC[BASIS_NON_NULL;IN_INSERT]]
    THEN SUBGOAL_THEN
      `orthogonal (u:real^3) (unit v) /\ orthogonal v (unit u)` 
      (ASSUME_TAC o REWRITE_RULE[orthogonal]) 
    THENL ON_FIRST_GOAL 
      (RULE_ASSUM_TAC (REWRITE_RULE[pairwise;orthogonal])
      THEN REWRITE_TAC[ORTHOGONAL_UNIT;orthogonal]
      THEN SUBGOAL_TAC "" `~(u=v:real^3)` [ASM_MESON_TAC[BASIS_DIFFERENT]]
      THEN CONJ_TAC THEN FIRST_ASSUM MATCH_MP_TAC THEN ASSUM_LIST SET_TAC)
    THEN ASM_SIMP_TAC[DOT_LADD;DOT_LMUL;REAL_MUL_RZERO;REAL_ADD_RID;
      REAL_ADD_LID;DOT_UNIT_SELF;GSYM VECTOR_MUL_ASSOC;GSYM UNIT_INTRO];
    REWRITE_TAC[unit;VECTOR_MUL_ASSOC] THEN MESON_TAC[]
  ]);;

let PLANE_DECOMP_DOT = 
prove
  (`!p:plane. plane p ==> !u v. is_orthogonal_plane_basis {u,v} p 
    ==> !n. is_normal_to_plane n p ==> !x:real^3.
      x = (x dot unit u) % unit u + (x dot unit v) % unit v
        + (x dot unit n) % unit n`,

  REPEAT STRIP_TAC
  THEN FIRST_ASSUM (STRIP_ASSUME_TAC o CONV_RULE (REWR_CONV 
    is_orthogonal_plane_basis))
  THEN ASM_CSQ_THEN PLANE_DECOMP_DOT_NORMAL 
    (C ASM_CSQ_THEN (C ASM_CSQ_THEN (STRIP_ASSUME_TAC o SPEC `x:real^3`)))
  THEN SUBGOAL_THEN
    `(x:real^3) dot unit u = w dot unit u /\ x dot unit v = w dot unit v`
    ASSUME_TAC
  THENL ON_FIRST_GOAL
  begin
    SUBGOAL_THEN `(u:real^3) IN {u,v} /\ v IN {u,v}` STRIP_ASSUME_TAC 
    THENL ON_FIRST_GOAL (ASSUM_LIST SET_TAC)
    THEN ASM_CSQ_THEN NORMAL_OF_PLANE_IS_NORMAL_TO_BASIS 
      (C ASM_CSQ_THEN ASSUME_TAC)
    THEN POP_ASSUM (fun x -> ASM_CSQ_THEN ~remove:true x 
      (C ASM_CSQ_THEN ASSUME_TAC) 
      THEN ASM_CSQ_THEN x (C ASM_CSQ_THEN ASSUME_TAC))
    THEN RULE_ASSUM_TAC (REWRITE_RULE[orthogonal] 
      o ONCE_REWRITE_RULE[GSYM ORTHOGONAL_LRUNIT]
      o ONCE_REWRITE_RULE[ORTHOGONAL_SYM])
    THEN ASM ONCE_REWRITE_TAC[]
    THEN ASM_REWRITE_TAC[DOT_LADD;DOT_LMUL;REAL_MUL_RZERO;REAL_ADD_RID]
  end
  THEN ASM_CSQ_THEN (REWRITE_RULE[GSYM IMP_IMP] 
    BASIS_OF_PLANE_ORTHONORMAL_DECOMPOS) (C ASM_CSQ_THEN (ASSUME_TAC
      o SPEC `w:real^3`))
  THEN ASSUM_LIST (GEN_REWRITE_TAC (RATOR_CONV o DEPTH_CONV))
  THEN ASM_REWRITE_TAC[VECTOR_ARITH `x+y=z+t+y <=> x=z+t`]
  THEN ASM_MESON_TAC[]);;

let UNIT_OF_ORTHOGONAL_BASIS_IS_ORTHOGONAL_BASIS = 
prove
  (`!p. plane p ==> !u v. is_orthogonal_plane_basis {u,v} p 
    ==> is_orthogonal_plane_basis {unit u,unit v} p`,

  SIMP_TAC[is_orthogonal_plane_basis;UNIT_OF_BASIS_IS_BASIS;pairwise;IN_INSERT;
    IN_SING]
  THEN REPEAT STRIP_TAC THEN ASM_REWRITE_TAC[]
  THEN ASM_MESON_TAC[ORTHOGONAL_LRUNIT]);;

let EXISTS_UNIT_ORTHOGONAL_PLANE_BASIS = 
prove
  (`!p. plane p ==> ?u v. is_orthogonal_plane_basis {unit u,unit v} p`,

  REPEAT STRIP_TAC
  THEN ASM_CSQ_THEN EXISTS_ORTHOGONAL_PLANE_BASIS (CHOOSE_THEN
    (DISTRIB (map (fun xs -> STRIP_ASSUME_TAC o REWRITE_RULE xs) [
      [];
      [is_orthogonal_plane_basis];
      [is_orthogonal_plane_basis;is_plane_basis]])))
  THEN FIRST_X_ASSUM SUBST_VAR_TAC 
  THEN ASM_MESON_TAC[UNIT_OF_ORTHOGONAL_BASIS_IS_ORTHOGONAL_BASIS]);;

let FORALL_PLANE_THM = 
prove
  (`!p. plane p ==>
    ?u v. is_orthogonal_plane_basis {unit u,unit v} p /\ !P pt0. pt0 IN p 
      ==> ((!pt. pt IN p ==> P pt)
        <=> (!a b. P (pt0 + a % unit u + b % unit v)))`,

  REPEAT STRIP_TAC 
  THEN ASM_CSQ_THEN EXISTS_UNIT_ORTHOGONAL_PLANE_BASIS (CHOOSE_THEN 
    (CHOOSE_THEN (DISTRIB (map (fun xs -> STRIP_ASSUME_TAC o REWRITE_RULE xs) [
      [];
      [is_orthogonal_plane_basis];
      [is_orthogonal_plane_basis;is_plane_basis]]))))
  THEN MAP_EVERY EXISTS_TAC [`u:real^3`;`v:real^3`]
  THEN ASM_REWRITE_TAC[] THEN REPEAT STRIP_TAC
  THEN ASM_CSQ_THEN PAIR_SETS_EQ STRIP_ASSUME_TAC
  THEN REPEAT (FIRST_X_ASSUM SUBST_VAR_TAC)
  THEN ASM_CSQ_THEN PLANE_AS_SPAN (C ASM_CSQ_THEN (C ASM_CSQ_THEN (ASSUME_TAC
    o REWRITE_RULE[SPAN_2;UNIV;IN_ELIM_THM]))) 
  THEN ASSUM_LIST SET_TAC);;

let FORALL_PLANE_THM_2 = 
prove
  (`!p. plane p ==>
    ?u v pt0. is_orthogonal_plane_basis {unit u,unit v} p /\ pt0 IN p
    /\ !P. (!pt. pt IN p ==> P pt)
      <=> (!a b. P (pt0 + a % unit u + b % unit v))`,

  REPEAT STRIP_TAC THEN ASM_CSQ_THEN FORALL_PLANE_THM STRIP_ASSUME_TAC
  THEN ASM_CSQ_THEN PLANE_NON_EMPTY (STRIP_ASSUME_TAC o REWRITE_RULE[GSYM MEMBER_NOT_EMPTY]) THEN ASM_MESON_TAC[]);;

(** Is u the projection of v on the hyperplane represented by the non-null vector w? *)
let is_projection_on_hyperplane = new_definition
  `is_projection_on_hyperplane u v (w:real^N)
  <=> norm w = &1 ==> u = v - (v dot w) % w`;;

let projection_on_hyperplane = new_definition
  `projection_on_hyperplane v w = @u. is_projection_on_hyperplane u v w`;;

let PROJECTION_ON_HYPERPLANE_THM = 
prove
  (`!v w. norm w = &1 
    ==> projection_on_hyperplane v w = v - (v dot w) % w`,

  REWRITE_TAC[projection_on_hyperplane;is_projection_on_hyperplane]
  THEN SELECT_ELIM_TAC THEN MESON_TAC[]);;

let PROJECTION_ON_HYPERPLANE_THM_LNEG = 
prove
  (`!v w. norm w = &1 
    ==> projection_on_hyperplane (--v) w = --(projection_on_hyperplane v w)`,

  SIMP_TAC[PROJECTION_ON_HYPERPLANE_THM;VECTOR_NEG_SUB;DOT_LNEG;
    VECTOR_MUL_LNEG]
  THEN REPEAT STRIP_TAC THEN VECTOR_ARITH_TAC);;

let PROJECTION_ON_HYPERPLANE_DECOMPOS = 
prove
  (`!v w. norm w = &1 ==> projection_on_hyperplane v w + (v dot w) % w = v`,
 
  SIMP_TAC[projection_on_hyperplane;is_projection_on_hyperplane]
  THEN REPEAT STRIP_TAC THEN VECTOR_ARITH_TAC);;

let symetric_vectors_wrt = new_definition
  `symetric_vectors_wrt u v (w:real^3) 
    <=> norm w = &1 ==> u + v = (&2 * u dot w) % w`;;

let SYMETRIC_VECTORS_WRT_SYM = 
prove
  (`!u v w. symetric_vectors_wrt u v w <=> symetric_vectors_wrt v u w`,

  REWRITE_TAC[symetric_vectors_wrt] THEN REPEAT GEN_TAC THEN EQ_TAC
  THEN DISCH_THEN (fun x -> DISCH_THEN (DISTRIB (map ((o) ASSUME_TAC)
    [MP x; REWRITE_RULE[NORM_EQ_1]])))
  THEN ASSUM_LIST (GEN_REWRITE_TAC (RAND_CONV o DEPTH_CONV)
    o map (REWRITE_RULE[VECTOR_ARITH `x+y=z <=> y=z-x:real^N`]))
  THEN ASM_REWRITE_TAC[DOT_LSUB;DOT_LMUL;REAL_ARITH `(&2 * x) * &1 - x = x`]
  THEN ONCE_REWRITE_TAC[VECTOR_ADD_SYM] THEN ASM_REWRITE_TAC[]);;

let SYMETRIC_VECTORS_PROJECTION_ON_HYPERPLANE = 
prove
  (`!u v w. norm w = &1 ==> symetric_vectors_wrt u v w 
    ==> --(projection_on_hyperplane v w) = projection_on_hyperplane u w`,

  ONCE_REWRITE_TAC[MATCH_MP (MESON[] `(p <=> q) ==> (p <=> p /\ q)`) 
    (SPEC_ALL SYMETRIC_VECTORS_WRT_SYM)]
  THEN SIMP_TAC[symetric_vectors_wrt;projection_on_hyperplane;
    is_projection_on_hyperplane;
    VECTOR_ARITH `x = (&2 * y) % z <=> y%z = inv(&2) % x`]
  THEN REPEAT STRIP_TAC THEN VECTOR_ARITH_TAC);;

let SYMETRIC_VECTORS_PROJECTION_ON_AXIS = 
prove
  (`!u v w. norm w = &1 ==> symetric_vectors_wrt u v w ==> u dot w = v dot w`,

  REPEAT STRIP_TAC
  THEN ASM_CSQ_THEN (GEQ_IMP SYMETRIC_VECTORS_WRT_SYM) ASSUME_TAC 
  THEN REPEAT (POP_ASSUM MP_TAC) THEN SIMP_TAC[symetric_vectors_wrt]
  THEN REWRITE_TAC[IMP_IMP;GSYM CONJ_ASSOC;REAL_EQ_MUL_LCANCEL;
    VECTOR_ARITH `u+v=w /\ v+u=x <=> u+v=x /\ x=w`;VECTOR_MUL_RCANCEL]
  THEN REPEAT STRIP_TAC THEN ASM_SIMP_TAC[DOT_RZERO] THEN POP_ASSUM MP_TAC
  THEN REAL_ARITH_TAC);;


(** Additions to affine spaces *)

let AFFINE_HULL_3_ZERO = 
prove
  (`affine hull {vec 0, a, b} = {u % a + v % b | u IN (:real) /\ v IN(:real)}`,

  REWRITE_TAC[AFFINE_HULL_3;UNIV;EXTENSION;IN_ELIM_THM;VECTOR_MUL_RZERO;
    VECTOR_ADD_LID;REAL_ARITH `x+y+z=t <=> x=t-y-z:real`]
  THEN MESON_TAC[]);;

let LSUB_PLANE_EQ_RSUB_PLANE = 
prove
  (`!p:plane. plane p ==> !x. x IN p ==> {y - x | y IN p} = {x - y | y IN p}`,

  let COMMON_TAC t =
    EXISTS_TAC t THEN CONJ_TAC THEN TRY VECTOR_ARITH_TAC
    THEN SUBGOAL_TAC "" `(x':real^3) IN span b` [ASM SET_TAC []]
  in
  REWRITE_TAC[GSPEC;SETSPEC;FUN_EQ_THM] THEN REPEAT STRIP_TAC
  THEN ASM_CSQ_THEN EXISTS_PLANE_BASIS STRIP_ASSUME_TAC
  THEN ASM_CSQS_THEN PLANE_AS_SPAN ASSUME_TAC THEN EQ_TAC THENL [
    ASM_CSQS_THEN PLANE_SPAN ASSUME_TAC THEN STRIP_TAC 
    THEN COMMON_TAC `x-x':point` THEN REWRITE_TAC[VECTOR_SUB]
    THEN POP_ASSUM (ASSUME_TAC o MATCH_MP SPAN_NEG) 
    THEN ASSUM_LIST (FIRST o List.rev_map (CHANGED_TAC o SINGLE REWRITE_TAC)) 
    THEN ASSUM_LIST SET_TAC;
    STRIP_TAC THEN ASM_CSQS_THEN PLANE_SPAN ASSUME_TAC 
    THEN COMMON_TAC `x+x':point` THEN ASSUM_LIST SET_TAC]);;

let EXISTS_MULTIPLE_OF_NORMAL_IN_PLANE = 
prove
  (`!p:plane. plane p ==> !n. is_normal_to_plane n p ==> ?a. a % unit n IN p`,

  REPEAT STRIP_TAC
  THEN ASM_CSQ_THEN PLANE_NON_EMPTY (STRIP_ASSUME_TAC 
    o REWRITE_RULE[GSYM MEMBER_NOT_EMPTY])
  THEN ASM_CSQ_THEN EXISTS_PLANE_BASIS STRIP_ASSUME_TAC
  THEN ASM_CSQS_THEN (GSYM PLANE_SPAN) (fun x -> 
    ASM_CSQS_THEN PLANE_DECOMP_DOT_NORMAL 
    (MP_TAC o SPEC `x:point` o REWRITE_RULE[x]))
  THEN ASM_SIMP_TAC[LSUB_PLANE_EQ_RSUB_PLANE] THEN REWRITE_TAC[IN_ELIM_THM]
  THEN STRIP_TAC THEN EXISTS_TAC `(x:point) dot unit n`
  THEN MATCH_MP_TAC (MESON[] `(pt':point) IN p /\ pt' = y ==> y IN p`) 
  THEN ASM_REWRITE_TAC[] THEN REPLICATE_TAC 2 (POP_ASSUM MP_TAC)
  THEN VECTOR_ARITH_TAC);;

(* Additional functions *)
let map_triple = new_definition `map_triple (f:A->B) (x1,x2,x3) = (f x1,f x2,f x3)`;;