(*We compute the components of integrals of eta the standard cusp form on various geodesic segments -- Proposition B.4*) 

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

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

(* The following map parameteized the Klein model from the upper half space model--computed by Hong Chan Kim in the reference *)

HKimMap[x_, y_] := {2x/(x^2+y^2+1), 1- 2/(x^2+ y^2 + 1), 1}

(*The begining point of the complete geodesic*)

VP:= HKimMap[0, 0]

VP

(* The end point of the complete geodesic ending at (0,0) and (0, -2R) in the upper half space model*)

VN:= HKimMap[-2R, 0]

Together[VN]

Simplify[VN.VN]

(* We define Lorentz Cross product of VN and VP *)

(* This was done by hand.. mathematica makes a strange error here *)

BV:=IJ.Cross[VN, VP]

Simplify[BV]

Simplify[LDOT[BV, BV]]

(* Normalize BV by taking the square roof of the above LDOT result *)

BVN:= Simplify[BV/(8R^2/(1+ 4R^2))]

BVN

(* Checking Normality*)

Simplify[LDOT[VP, BV]]

Simplify[LDOT[VN, BV]]

(* Parabolic coordinate $1$-form vector here. This is in the parabolic coordinates (See StandardFormII.nb *)

OurVect :={  mu x^2/2, + mu x, mu}

(* Matrix of change of coordinates to the usual one *)

ChangeM:=Transpose[{{0, -1/Sqrt[2], 1/Sqrt[2]}, {1, 0, 0}, {0, 1/Sqrt[2], 1/Sqrt[2]}}]

TeXForm[MatrixForm[ChangeM]]

(* This is the inner product matrix for parabolic coorinate *)

Mu:={{0, 0, -1}, {0, 1, 0}, {-1, 0, 0}}

(* Check the Lorentz inner product matrix becomes Mu under our coordinate chage *)

Simplify[Transpose[ChangeM].IJ. ChangeM]

(* Write our cusp $1$-form vector into our coorinates *)

Simplify[ChangeM.OurVect]

(* Write OurVector as a linear combination of the vectors *)

Simplify[a1 VN + a2 VP + a3 BVN]

(* Solving the get the coefficients a1, a2, a3 *)

Solve[{a3-((4 a1 R)/(1+4 (R^2)))==mu x,a1-a2+(a3/(2 R))-((2 a1)/(1+4 (R^2)))==((mu (1-((x^2)/2)))/(Sqrt[2])),a1+a2-(a3/(2 R))==((mu (2+(x^2)))/(2 (Sqrt[2])))}, {a1, a2, a3}]

Factor[Collect[{{a1->((mu (1+4 (R^2)))/(4 (Sqrt[2]) (R^2))),a2->((mu ((Sqrt[2])+4 R x+2 (Sqrt[2]) (R^2) (x^2)))/(8 (R^2))),a3->-((-(Sqrt[2]) mu-2 mu R x)/(2 R))}}, R]]

Collect[{{a1->((mu (1+4 (R^2)))/(4 (Sqrt[2]) (R^2))),a2->((mu ((Sqrt[2])+4 R x+2 (Sqrt[2]) (R^2) (x^2)))/(8 (R^2))),a3->((mu ((Sqrt[2])+2 R x))/(2 R))}}, x]

(* V_- component behavior *)

Integrate[ ((mu (1+4 (R^2)))/(4 (Sqrt[2]) (R^2))), {x, R - Sqrt[R^2-1], R+ Sqrt[R^2-1]}]

((mu (Sqrt[-1+(R^2)]) (1+4 (R^2)))/(2 (Sqrt[2]) (R^2)))/.{mu->1}

Plot[(((Sqrt[-1+(R^2)]) (1+4 (R^2)))/(2 (Sqrt[2]) (R^2))), {R, 1, 10}, PlotRange->{{1, 15},{0,15}}]

(*VP component integration *)

Collect[Integrate[(mu/(4 (Sqrt[2]) (R^2)))+((mu x)/(2 R))+((mu (x^2))/(2 (Sqrt[2]))),  {x, R - Sqrt[R^2-1], R+ Sqrt[R^2-1]}], R]

Factor[mu (Sqrt[-1+(R^2)])-((mu (Sqrt[-1+(R^2)]))/(3 (Sqrt[2])))+((mu (Sqrt[-1+(R^2)]))/(2 (Sqrt[2]) (R^2)))+(2/3) (Sqrt[2]) mu (R^2) (Sqrt[-1+(R^2)])]

Simplify[(((mu (Sqrt[-1+(R^2)]) (3 (Sqrt[2])+12 (R^2)-2 (Sqrt[2]) (R^2)+8 (Sqrt[2]) (R^4)))/(12 (R^2))))/.{b2->1, b3->1 , mu->1} ]

Plot[{(((Sqrt[-1+(R^2)]) (3 (Sqrt[2])-2 (-6+(Sqrt[2])) (R^2)+8 (Sqrt[2]) (R^4)))/(12 (R^2)))}, {R, 1, 10}]

(* Neutral direction integration *)

Integrate[((mu ((Sqrt[2])+2 R x))/(2 R)), {x, R - Sqrt[R^2-1], R+ Sqrt[R^2-1]}]

Factor[(((Sqrt[2]) mu (Sqrt[-1+(R^2)]))/R)+2 mu R (Sqrt[-1+(R^2)])]

Simplify[{((mu (Sqrt[-1+(R^2)]) ((Sqrt[2])+2 (R^2)))/R)}/.{mu-> 1}]

Plot[{(((Sqrt[-1+(R^2)]) ((Sqrt[2])+2 (R^2)))/R)}, {R, 1, 10}]

(* Last part for (4.45), we are using the invariance of the value under the parabolic transformations. Hence on r and R matters. We take the difference when r value changes to r'=r+eps', eps' leq 1.1 eps, Here eps can be negative *) 
(* We use R = Sqrt[r^2+1], r > 10, r < R = sqrt[r^2+1]< r+1, 1/R < 1/10, 1/R < r/100, eps < 1/8, eps < r/80,  r' > 9, R'> 9, r' < R' = sqrt[r'^2+1]< r'+1, 1/R' < 1/9,  1/R' < r/90 *)

(*We need to estimate the absolute value of*)

((((Sqrt[2]) r)/R)+2  R r) - ( (((Sqrt[2]) r')/(R'))+2  R' r')

(* This is =< Asolute value of *)

(Sqrt[2] r/R + R*r) - (Sqrt[2] r/ R' + Sqrt[2] eps'/ R' + R'*r + R'*eps')

(* Now use R' leq R+ Abs[eps'] since the slope of Sqrt[r^2+1] is less than 1, The above is =< Absolute value of *)

Abs[Sqrt[2] r (1/R - 1/R')] + Abs[(R-R')*r ]+Sqrt[2]* Abs[eps']/R' + (R+Abs[eps'])Abs[eps']

Abs[eps'] (R+Abs[eps'])+Abs[(Sqrt[2]) r ((1/R)-(1/(R')))]+Abs[r (R-R')]+(((Sqrt[2]) Abs[eps'])/(R'))

(* Convert to*)

Abs[eps' ](R+Abs[eps'])+(Sqrt[2]) r *Abs[((1/R)-(1/(R')))]+r*Abs[(R-R')]+Abs[(((Sqrt[2]) eps')/(R'))]

(* Turn (1/R - 1/R') = (R - R')/RR' leq Abs[eps']/RR' since R=Sqrt[r^2+1} is distance decreasing. Hence the above is less than or equal to*)

Abs[eps'] (R+Abs[eps'])+(Sqrt[2]) r*Abs[eps'/(R R')]+r*Abs[eps']+(Sqrt[2]) Abs[((eps')/(R'))]

(* Now use R leq r+1 *)

r*Abs[eps']+Abs[eps' ](r+1+ Abs[eps'])+(Sqrt[2]) Abs[((eps')/(R'))]+(Sqrt[2]) r*Abs[(( eps')/(R R'))]

(* Collect r* eps'*)

r*Abs[eps']*( 1 + 1+ Sqrt[2]/(R*R')) + Abs[eps'](1+  Abs[eps']+ 1/R')

N[Sqrt[2]]

(* Now use eps' < (1.1)/8 < (1.1)r/80, 1/(RR') < 1/(10*9), 1/R'< 1/9 < r/90. The above is less than or equal to *)

Abs[eps' ](2+ (1.42)/90) *r + Abs[eps'](r/10 + (1.1)*r/80 + r/90)

(* The following is what we can use as the factor for eps*)

N[2.14064*(1.1)]