(* Section 2.2.1 Divergence of geodesics--Computes two geodesics near each other*)

(* Lorentz inner prodcut *)

LDOT[VA_, VB_]:= VA[[1]]VB[[1]] + VA[[2]]VB[[2]] - VA[[3]]VB[[3]]

(*Lorentz metric matrix *)

IJ={{1, 0, 0}, {0, 1, 0}, {0, 0, -1}}

(* LCross is IJ.Cross *) 

(*Using Ratcliffe The 3.5. 7 * x, y vertex, u, v perpendicular vertex  u, v on a line L_1 given by x2=0, passing (-1, 0, 1) and (1, 0, 1) *) 
(* Vetex, P1, P2 at the perp angles *)

P1:={0, 0, 1} 
P2:={Sinh[l1], 0, Cosh[l1]}

(*Isometry acting on the line L_1*)

mA1[t_]:={{Cosh[t] , 0, Sinh[t]}, {0, 1, 0}, {Sinh[t], 0, Cosh[t]}}

(* first geodesic standard one on L_1*) 
eta1[t_] := Cosh[t]{0, 0,1} + Sinh[t]{1, 0, 0}

(* Geodesic zeta1 perpendicular to L_1 at P1 *)

zeta1[t_] := Cosh[t]{0,0,1} + Sinh[t]{0, 1, 0}

(* The geodesic zeta 2 perpendicular to L_1 at P2 *)

zeta2[t_] := mA1[l1].zeta1[t]

(* zeta[s, t] the geodesic passing eta1[s] perdicular to L_1 *)

zeta[s_, t_]:= mA1[s].zeta1[t]

(*Q1, Q2 vertex with angle alpha, beta Q1 on zeta1 distance E1, Q2 on zeta2 distance E2. First case both are in the same side of g1*) 
Q1:= {0, Sinh[E1], Cosh[E1]}

Q2=zeta2[E2]

(* Find a vector Lorentz perpendicualr to Q1 Q2. Determines the Line L_2 through Q1 and Q2*)

NQ1Q2:=IJ.Cross[Q1, Q2]

Simplify[NQ1Q2]

(* Find a point on the line perpendicualr to eta1(t) and meeting line L_2 containing Q1 and Q2 *)

Simplify[LDOT[NQ1Q2, zeta[t, eps]]]

Simplify[Solve[-Cosh[E2] Cosh[eps] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[eps] Sinh[l1]-Cosh[eps] Sinh[E2] Sinh[t])==0, eps, Reals]]

Simplify[(((Cosh[E1] Cosh[E2] Sinh[l1]+Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] Sinh[E2] Sinh[t])/(Cosh[E1] Cosh[E2] Sinh[l1]-Cosh[E2] Sinh[E1] Sinh[l1-t]-Cosh[E1] Sinh[E2] Sinh[t])))]

ExpandAll[((Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]+Sinh[E2] Sinh[t]))/(-Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]-Sinh[E2] Sinh[t])))]

(* The vector in the direction of the point on line L_2 and the line perp to eta1(t) is as follows *)

VV:= zeta[t,(1/2)* Log[(((Cosh[E1] Cosh[E2] Sinh[l1]+Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] Sinh[E2] Sinh[t])/(Cosh[E1] Cosh[E2] Sinh[l1]-Cosh[E2] Sinh[E1] Sinh[l1-t]-Cosh[E1] Sinh[E2] Sinh[t])))]]

Simplify[VV]

VVp[t_]:= {Cosh[(1/2) Log[((Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]+Sinh[E2] Sinh[t]))/(-Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]-Sinh[E2] Sinh[t])))]] Sinh[t],Sinh[(1/2) Log[((Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]+Sinh[E2] Sinh[t]))/(-Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]-Sinh[E2] Sinh[t])))]],Cosh[t] Cosh[(1/2) Log[((Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]+Sinh[E2] Sinh[t]))/(-Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]-Sinh[E2] Sinh[t])))]]}

(* verify that VV is one the hyperbola *)

Simplify[LDOT[VV, VV]]

((VV[[1]])^2)+((VV[[2]])^2)-((VV[[3]])^2)

(* We compute the distance from the perp plane to (0,1,0) to VV by taking the inner product and taking the inverse sinh. See Ratcliffe.P.75 *)

Simplify[LDOT[VV, {0, 1, 0}]]

(* Our distance function is given as follows *)

DistF[t_]= (1/2) Log[((Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]+Sinh[E2] Sinh[t]))/(-Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]-Sinh[E2] Sinh[t])))]

(* Some examples of the distance functions *)

Plot[{(1/2) Log[((Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]+Sinh[E2] Sinh[t]))/(-Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]-Sinh[E2] Sinh[t])))]}/.{l1-> 30, E1-> 0.1, E2-> 0.2}, {t, 0, 30}, PlotRange->{0, 0.2}]

Simplify[ExpandAll[(-2 Cosh[(1/2) (l1-2 t)] Sinh[E2] Sinh[(l1/2)]-Cosh[E1] Sinh[l1]) (2 Cosh[(1/2) (l1-2 t)] Sinh[E2] Sinh[(l1/2)]-Cosh[E1] Sinh[l1])]]

(* Now we look at E2 being on the other side -- This is not really needed, but we are checking*)

Q2p=zeta2[-E2]

(* Find a vector perpendicular to Q1 Q2p *)

NQ1Q2p:=IJ.Cross[Q1, Q2p]

Simplify[NQ1Q2p]

Simplify[LDOT[NQ1Q2p, zeta[t, eps]]]

Simplify[Solve[-Cosh[E2] Cosh[eps] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[eps] Sinh[l1]+Cosh[eps] Sinh[E2] Sinh[t])==0, eps, Reals]]

VVpp:= zeta[t,(1/2) Log[((Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]-Sinh[E2] Sinh[t]))/(-Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]+Sinh[E2] Sinh[t])))]]

Solve[((Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]-Sinh[E2] Sinh[t]))/(-Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]+Sinh[E2] Sinh[t])))==1, t, Reals]

Simplify[LDOT[VVpp, VVpp]]

(* We compute the distance from the perp plane to (0,1,0) to VV by taking the inner product and taking the inverse sinh. See Ratcliffe.P.75 *)

Simplify[LDOT[VVpp, {0, 1, 0}]]

DistFp[t_]:= (1/2) Log[((Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]-Sinh[E2] Sinh[t]))/(-Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]+Sinh[E2] Sinh[t])))]

(* Example of the distance *)

Plot[((1/2) Log[((Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]-Sinh[E2] Sinh[t]))/(-Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]+Sinh[E2] Sinh[t])))])/.{l1-> 30, E1-> 0.1, E2-> 0.2} , {t, 0, 30}, PlotRange->{-0.3, 0.3}]

Plot[(((Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]-Sinh[E2] Sinh[t]))/(-Cosh[E2] Sinh[E1] Sinh[l1-t]+Cosh[E1] (Cosh[E2] Sinh[l1]+Sinh[E2] Sinh[t]))))/.{l1-> 30, E1-> 0.2, E2-> 0.2} , {t, 0, 30}, PlotRange->{0, 4}]