(* ::Package:: *)

(* ::Title:: *)
(*Hyperbolic Toolbox Version 1.0*)

(* ::Subsubtitle:: *)
(*William Goldman*)
(*University of Maryland*)
(*http://egl.math.umd.edu/*)
(**)
(*Edited by Ryan Hoban*)
(*August  2007*)
(**)
(*Edited by Brian Baird*)
(*August 2010*)

(* ::Section:: *)
(*Introduction*)

(* ::Text:: *)
(*This toolbox is a set of Mathematica functions designed to:  *)
(**)
(*	\[Bullet] Draw hyperbolic geodesics and circles in the Poincar\[EAcute] unit disc as well as the Upper Half Plane model.*)
(*	*)
(*	\[Bullet] Compute and apply linear fractional transformations of the isometry group PSL(2,\[DoubleStruckCapitalR]).*)
(*	*)
(*	\[Bullet]  Map from the Upper Half Plane to the Unit Disc and vice versa.*)
(*	*)
(*	\[Bullet] Allow the user to compute distances, label points, and other useful tools for creating detailed and interactive projects in hyperbolic geometry.*)
(**)

(* ::Subsubtitle:: *)
(*To easily navigate this package, use the section drop-down box above.*)

(* ::Section:: *)
(*How to Use this Package*)

(* ::Text:: *)
(*Place this file in the same directory as your project notebook.  You can then load the package easily by inserting the following code at the beginning of your notebook:*)
(**)
(*SetDirectory[NotebookDirectory[]];*)
(*<< HyperbolicToolbox_v1.0.m*)
(**)
(*Then go to Cell\[RightArrow]Cell Properties and click on "Initialization Cell".  This will ensure the package will be loaded when you evaluate your notebook.*)

(* ::Section:: *)
(*Functions and their Usage*)

(* ::Text:: *)
(*For each function, the type of input required and output generated will be listed.  For example, (\[DoubleStruckCapitalH]\[RightArrow]\[DoubleStruckCapitalD], Graphics) would indicate the function takes objects in the Upper Half Plane to the Unit Disc and generates graphics objects, while  (\[DoubleStruckCapitalH]\[RightArrow]\[DoubleStruckCapitalD]) would indicate the function takes objects in the Upper Half Plane to the Unit Disc in the form of complex numbers.*)

(* ::Section:: *)
(*Graphics Functions*)

(* ::Text:: *)
(*Unless otherwise noted, all graphics functions should be used within the Graphics command.*)

(* ::Subsection:: *)
(*PGSegment[{z1,z2}]*)

(* ::Text:: *)
(*(\[DoubleStruckCapitalH] \[RightArrow] \[DoubleStruckCapitalD], Graphics) Given two points from the Upper Half Plane, PGSegment draws hyperbolic lines on the disc model.  An optional third argument specifies the minimum draw size.*)
(**)

(* ::Subsection:: *)
(*UHPSegment[{z1,z2}]*)

(* ::Text:: *)
(*(\[DoubleStruckCapitalH] \[RightArrow] \[DoubleStruckCapitalH], Graphics) Draws hyperbolic lines in the Upper Half Plane.*)
(**)

(* ::Subsection:: *)
(*PCircle[Center,Radius]*)

(* ::Text:: *)
(*(\[DoubleStruckCapitalH] \[RightArrow] \[DoubleStruckCapitalD], Graphics) Given a center in the Upper Half Plane and a radius, PCircle draws a hyperbolic circle on the Unit Disc.  PCircle[Center, Radius, True] will draw the hyperbolic center as well as the circle.*)
(**)

(* ::Subsection:: *)
(*UHPCircle[Center,Radius]*)

(* ::Text:: *)
(*(\[DoubleStruckCapitalH] \[RightArrow]\[DoubleStruckCapitalH], Graphics) Given a center in the Upper Half Plane and a radius, UHPCircle draws a hyperbolic circle in the Upper Half Plane.  UHPCircle[Center, Radius, True] will draw the hyperbolic center as well as the circle.*)
(**)

(* ::Subsection:: *)
(*UDPoint[z]*)

(* ::Text:: *)
(*(\[DoubleStruckCapitalH] \[RightArrow] \[DoubleStruckCapitalD], Graphics) Takes a point z in the Upper Half Plane and returns the Graphics primitive for this point in the Unit Disc.*)
(**)

(* ::Subsection:: *)
(*UDText[z,text]*)

(* ::Text:: *)
(*(\[DoubleStruckCapitalH], Text \[RightArrow] \[DoubleStruckCapitalD], Text) places text in the Unit Disc according to the point z in the Upper Half Plane.*)
(**)

(* ::Subsection:: *)
(*ShowP[List]*)

(* ::Text:: *)
(*Generates the Unit Disc and displays a list of Graphics primitives.  Use in place of Graphics when drawing in the Poincar\[EAcute] disc.*)
(**)

(* ::Subsection:: *)
(*Circumcircle*)

(* ::Text:: *)
(*Circumcircle draws the circle passing through three non-collinear points.  The input can be six real numbers, three ordered pairs, or three complex numbers.*)

(* ::Section:: *)
(*Mapping between the Upper Half Plane and the Unit Disc*)

(* ::Subsection:: *)
(*ToUnitDisc[z]*)

(* ::Text:: *)
(*(\[DoubleStruckCapitalH] \[RightArrow] \[DoubleStruckCapitalD]) maps the point z in the Upper Half Plane to the corresponding point in the Unit Disc.*)
(**)

(* ::Subsection:: *)
(*ToUHP*)

(* ::Text:: *)
(* (\[DoubleStruckCapitalD] \[RightArrow] \[DoubleStruckCapitalH]) or takes complex numbers or ordered pairs from the Unit Disk to the Upper Half Plane.*)
(**)

(* ::Subsection:: *)
(*ToR2[z]*)

(* ::Text:: *)
(*( \[DoubleStruckCapitalC] \[RightArrow] \[DoubleStruckCapitalR]\.b2 )  Returns an ordered pair from the real and imaginary parts of z.  Use this function when a function requires an ordered pair instead of a complex number (such as most Mathematica Graphics primitives).*)
(**)

(* ::Section:: *)
(*Linear Fractional Transformations*)

(* ::Subsection:: *)
(*ApplyLFT[M,z]*)

(* ::Text:: *)
(*(\[DoubleStruckCapitalH] \[RightArrow] \[DoubleStruckCapitalH]) Applies the linear fractional transformation corresponding to the 2x2 matrix M to the extended complex number z.  If z is a list (or list of lists) then ApplyLFT applies the linear fraction across the list.*)

(* ::Subsection:: *)
(*FindLFT[{z0, z1, zInf}, {w0, w1, wInf}]*)

(* ::Text:: *)
(*Returns a 2x2 matrix in PSL (2, \[DoubleStruckCapitalR]) whose linear fractional transformation takes z1 to w1, z2 to w2, and z3 to w3.  If the transformation is not in PSL (2, \[DoubleStruckCapitalR]), it will return ' ERROR'."*)

(* ::Subsection:: *)
(*InvolutionFixing[z1, z2]*)

(* ::Text:: *)
(*Returns the 2x2 matrix representing an involution fixing complex numbers z1, z2. It is normalized to have determinant - 1.*)

(* ::Subsection:: *)
(*ReflectIn[z1,z2]*)

(* ::Text:: *)
(*Reflects the complex number z (upper half - plane coordinates) in the Poincar\[EAcute] geodesic with endpoints z1, z2 (extended real numbers). This is the effect of applying an anti - involution (anti - homography of order two) which fixes the chain orthogonal to the unit circle (real axis in upper half - plane coordinates) whose endpoints on the unit circle correspond to extended real numbers z1, z2.*)

(* ::Subsection:: *)
(*SymmetryIn[z]*)

(* ::Text:: *)
(*Returns the involutary isometry of the Upper Half Plane fixing z.*)

(* ::Subsection:: *)
(*FixedPointsSet[M]*)

(* ::Text:: *)
(*Returns extended real numbers, representing the attracting fixed point, then the repelling fixed point of the linear fractional transformation represented by the 2x2 matrix M.*)

(* ::Subsection:: *)
(*Commutator[A,B]*)

(* ::Text:: *)
(*Returns the commutator of the matrices A and B.*)

(* ::Subsection:: *)
(*CrossRatioMatrix[{z0,z1,zInf}]*)

(* ::Text:: *)
(*Returns the matrix whose linear fractional transformation takes z0 to 0, z1 to 1, zInf to ComplexInfinity.*)

(* ::Subsection:: *)
(*HypElement[{x1,x2}]*)

(* ::Text:: *)
(*Returns some Hyperbolic Element in PSL (2, R) with ideal fixed points x1 and x2.*)

(* ::Section:: *)
(*Computing with Geodesics*)

(* ::Subsection:: *)
(*HyperbolicLength[z1,z2]*)

(* ::Text:: *)
(*(\[DoubleStruckCapitalH] \[RightArrow] \[DoubleStruckCapitalR]) Returns the length of the geodesic segment connecting z1 and z2.*)

(* ::Subsection:: *)
(*MidpointEndpointsPGeodesic[z1, z2]*)

(* ::Text:: *)
(*(\[DoubleStruckCapitalH] \[RightArrow] \[DoubleStruckCapitalH]) Returns the midpoint of the pair of endpoints of the Poincar\[EAcute] geodesic containing points z1 and z2 in the Upper Half Plane." *)

(* ::Subsection:: *)
(*EndpointsPGeodesic[z1,z2]*)

(* ::Text:: *)
(*(\[DoubleStruckCapitalH] \[RightArrow] \[DoubleStruckCapitalH]) Returns the endpoints of the Poincar\[EAcute] geodesic containing z1 and z2 in the Upper Half Plane.*)

(* ::Subsection:: *)
(*IntersectionPoint[{z1,z2},{w1,w2}]*)

(* ::Text:: *)
(*(\[DoubleStruckCapitalH] \[RightArrow] \[DoubleStruckCapitalH]) Returns the complex number that is the point of intersection of the geodesics whose Ideal End Points are {z1, z2} and {w1, w2}.*)

(* ::Subsection:: *)
(*PerpendicularGeodesic[{z1, z2},{w1, w2}] *)

(* ::Text:: *)
(*(\[DoubleStruckCapitalH] \[RightArrow] \[DoubleStruckCapitalH]) Returns the ideal endpoints of the geodesic perpendicular to the geodesics with endpoints {z1, z2} and {w1, w2}.*)

(* ::Subsection:: *)
(*PerpendicularBisector[z1,z2]*)

(* ::Text:: *)
(*(\[DoubleStruckCapitalH] \[RightArrow] \[DoubleStruckCapitalH]) Returns the ideal endpoints of the line that is the perpendicular bisector of the line segment connecting z1 and z2.*)

(* ::Section:: *)
(*Other Functions*)

(* ::Subsection:: *)
(*Circumcenter*)

(* ::Text:: *)
(*Returns the circumcenter of three noncollinear points.  When the input is six real numbers or three ordered pairs, Circumcenter outputs an ordered pair.  When the input is three complex numbers, Circumcenter outputs a complex number.*)

(* ::Subsection:: *)
(*Affine[List]*)

(* ::Text:: *)
(*Returns the inhomogeneous coordinates of the list of homogeneous coordinates of a point in projective space.*)

(* ::Subsection:: *)
(*Projective[z]*)

(* ::Text:: *)
(*Returns the homogeneous coordinates of a point in P^1 corresponding to the extended complex number z.*)

(* ::Subsection:: *)
(*Adj[M]*)

(* ::Text:: *)
(*Returns the adjugate of the 2x2 matrix M.  Adj[M] = Det[M] Inverse[M].*)

(* ::Section::Closed:: *)
(*Code for functions*)

(* ::Item:: *)
(*Turn off annoying error messages*)

Off[General::"spell1"] 
Off[Power::"infy"]
Off[General::"spell"]
Off[Simplify::"infd"]


(* ::Item:: *)
(*Drawing Geodesics, Circles, and Points*)

PGSegment::usage="PGSegment[{z1,z2}] (\[DoubleStruckCapitalH]\[RightArrow]\[DoubleStruckCapitalD], Graphics) Given two points from the Upper Half Plane, PGSegment draws hyperbolic lines on the disc model.  An optional third argument specifies the minimum draw size.";
PGSegment[{z1_,z2_}]:=PGSegment[{z1,z2},.003]
PGSegment[{z1_,z2_},Threshold_]:=Module[{Disc1,Disc2,x1,y1,x2,y2,Output},
(*Move the coordinates to the disc model.*)
Disc1=ToUnitDisc[z1];
Disc2=ToUnitDisc[z2];

(*Convert to x,y coordinates for faster calculation.*)
x1=Re[Disc1];
y1=Im[Disc1];
x2=Re[Disc2];
y2=Im[Disc2];

(*If the coordinates are close to being antipodal, we wish to draw a straight line instead of a circular arc.*)
If[Chop[y1*x2-x1*y2]==0,Output=Line[{{x1,y1},{x2,y2}}],
(*Otherwise, we draw an arc of a circle tangent to the Unit Disc.*)
Output=DrawArc[x1,y1,x2,y2,Threshold]];
Return[Output]
]
DrawArc[x1_,y1_,x2_,y2_,Threshold_]:=Module[{CenterX,CenterY,Radius,Angle1,Angle2,BigAngle,LittleAngle,Output},

(*Calculates the center of the tangent circles using the formulas for circumcenter and the inversion of a point.  Simplified for the sake of brevity and speed.*)
CenterX=((1+x1^2) y2+y1^2 y2-y1 (1+x2^2+y2^2))/(-2 x2 y1+2 x1 y2);
CenterY=(x2+x1^2 x2+x2 y1^2-x1 (1+x2^2+y2^2))/(2 x2 y1-2 x1 y2);

(*Calculates the radius of the tangent circle*)
Radius=Norm[{x1-CenterX,y1-CenterY}];

(*Translates the angles for purposes of measurement, and converts angles in the 3rd and 4th quadrants to positive values greater than \[DoubledPi].*)
If[Chop[y1-CenterY]>=0,Angle1=Abs[ArcTan[x1-CenterX,y1-CenterY]],Angle1=ArcTan[x1-CenterX,y1-CenterY]+2*Pi];
If[Chop[y2-CenterY]>=0,Angle2=Abs[ArcTan[x2-CenterX,y2-CenterY]],Angle2=ArcTan[x2-CenterX,y2-CenterY]+2*Pi];

(*Orders the angles for ease of use.*)
If[Angle1<Angle2,LittleAngle=Angle1;BigAngle=Angle2;,LittleAngle=Angle2;BigAngle=Angle1;];

(*Places conditions for drawing the arc given the resulting angles.  Because Mathematica is finicky with input coordinates, several cases must be considered to keep the final angle less than \[DoubledPi].*)
If[Radius<Threshold,Output={}];
If[Radius>Threshold&&LittleAngle!=0&&BigAngle-LittleAngle<Pi,Output=Circle[{CenterX,CenterY},Radius,{LittleAngle,BigAngle}]];
If[Radius>Threshold&&LittleAngle!=0&&BigAngle-LittleAngle>Pi,Output=Circle[{CenterX,CenterY},Radius,{BigAngle-2*Pi,LittleAngle}]];
If[Radius>Threshold&&Chop[LittleAngle]==0&&BigAngle<Pi,Output=Circle[{CenterX,CenterY},Radius,{LittleAngle,BigAngle}]];
If[Radius>Threshold&&Chop[LittleAngle]==0&&BigAngle>Pi,Output=Circle[{CenterX,CenterY},Radius,{BigAngle,2*Pi}]];
(*The Threshold value is used to avoid drawing arcs that are too small to see.*)

Return[Output]
]

UHPSegment::usage="UHPSegment[{z1,z2}] (\[DoubleStruckCapitalH]\[RightArrow]\[DoubleStruckCapitalH],Graphics) draws hyperbolic lines in the Upper Half Plane.";
UHPSegment[{z1_,z2_}]:=Module[
{Output},
If[Abs[Re[z1-z2]]<.05,Output=Line[{ToR2[z1],ToR2[z2]}],Output=UHPArc[z1,z2]];
Return[Output]
]
UHPSegment[{z1_,ComplexInfinity}]:=Line[{ToR2[z1],{Re[z1],100000}}];
UHPSegment[{ComplexInfinity,z2_}]:=Line[{ToR2[z2],{Re[z2],100000}}];
UHPArc[z1_,z2_]:=Module[
{x1,y1,x2,y2,Center,Radius,Angle1,Angle2,Output},
x1=Re[z1];
y1=Im[z1];
x2=Re[z2];
y2=Im[z2];

Center =(x1^2-x2^2+y1^2-y2^2)/(2 (x1-x2));

Radius=Norm[{x1-Center,y1}];

Angle1=N[ArcTan[x1-Center,y1]];
Angle2=N[ArcTan[x2-Center,y2]];

If[Angle1>Angle2,Output={Circle[{Center,0},Radius,{Angle2,Angle1}],Point[{Center,0}]},Output={Circle[{Center,0},Radius,{Angle1,Angle2}],Point[{Center,0}]}];

Return[Output]
]

PCircle::usage="PCircle[Center,Radius] (\[DoubleStruckCapitalH]\[RightArrow]\[DoubleStruckCapitalD],Graphics) Given a center in the Upper Half Plane and a radius, PCircle draws a hyperbolic circle on the Unit Disc.  PCircle[Center,Radius,True] will draw the hyperbolic center as well as the circle.";
PCircle[Center_,Radius_,ShowCenter_]:=Module[
{P1,DiscP1,P2,P3,T,TP1,TP2,TP3,DiscCenter,Output},
(*First we assume the center of the circle in the Upper Half Plane is I.  Two points lying on the circle are easily found on the line from 0 to ComplexInfinty.*)
P1=E^Radius*I;
P2=E^-Radius*I;

(*We find the third point by mapping one point lying on the circle to the unit disc.  When mapped to the Unit Disc, the center of the circle will be the origin, which is also the Euclidean center.  Using the center and Euclidean radius, we can find a third point on the circle and map it back to the Upper Half Plane.*)
DiscP1=ToUnitDisc[P1];
P3=ToUHP[Norm[DiscP1]];

(*Now we must find the linear fractional transformation taking I to the center we specified.*)
T=FindLFT[{0,I,ComplexInfinity},{Re[Center],Center,ComplexInfinity}];

(*Now we transform the points and send them to the Unit Disc*)
TP1=ToUnitDisc[ApplyLFT[T,P1]];
TP2=ToUnitDisc[ApplyLFT[T,P2]];
TP3=ToUnitDisc[ApplyLFT[T,P3]];
DiscCenter=ToUnitDisc[Center];

(*If Showcenter=True, a point will be displayed on the Unit Disc representing the hyperbolic center.*)
If[ShowCenter==True,Output={Circumcircle[TP1,TP2,TP3],Point[ToR2[DiscCenter]]},Output=Circumcircle[TP1,TP2,TP3]];
Return[Output]
]
PCircle[Center_,Radius_]:=PCircle[Center,Radius,False]

UHPCircle::usage="UHPCircle[Center,Radius] (\[DoubleStruckCapitalH]\[RightArrow]\[DoubleStruckCapitalH],Graphics) Given a center in the Upper Half Plane and a radius, UHPCircle draws a hyperbolic circle in the Upper Half Plane.  UHPCircle[Center,Radius,True] will draw the hyperbolic center as well as the circle.";
UHPCircle[Center_,Radius_,ShowCenter_]:=Module[{P1,DiscP1,P2,P3,T,TP1,TP2,TP3,Output},
(*First we assume the center of the circle in the Upper Half Plane is I.  Two points lying on the circle are easily found on the line from 0 to ComplexInfinty.*)
P1=E^Radius*I;
P2=E^-Radius*I;

(*We find the third point by mapping one point lying on the circle to the unit disc.  When mapped to the Unit Disc, the center of the circle will be the origin, which is also the Euclidean center.  Using the center and Euclidean radius, we can find a third point on the circle and map it back to the Upper Half Plane.*)
DiscP1=ToUnitDisc[P1];
P3=ToUHP[Norm[DiscP1]];

(*Now we must find the linear fractional transformation taking I to the center we specified.*)
T=FindLFT[{0,I,ComplexInfinity},{Re[Center],Center,ComplexInfinity}];

(*Now we transform the three points lying on the circle*)
TP1=ApplyLFT[T,P1];
TP2=ApplyLFT[T,P2];
TP3=ApplyLFT[T,P3];

(*If Showcenter=True, a point will be displayed in the Upper Half Plane representing the hyperbolic center.*)
If[ShowCenter==True,Output={Circumcircle[TP1,TP2,TP3],Point[ToR2[Center]]},Output=Circumcircle[TP1,TP2,TP3]];
Return[Output]
]
UHPCircle[Center_,Radius_]:=UHPCircle[Center,Radius,False]

UDPoint::usage="UDPoint[z] (\[DoubleStruckCapitalH]\[RightArrow]\[DoubleStruckCapitalD],Graphics) takes a point z in the Upper Half Plane and returns the Graphics primitive for this point in the Unit Disc.";
UDPoint[z_]:=Module[{udpoint,x,y},
udpoint=ToUnitDisc[z];
x=Re[udpoint] //N;
y=Im[udpoint] //N;
Point[{x,y}]
];

UDText::usage="UDText[text,z] (\[DoubleStruckCapitalH],Text\[RightArrow]\[DoubleStruckCapitalD],Text) places text in the Unit Disc according to the point z in the Upper Half Plane.";
UDText[text_,z_]:=Module[{udpoint,x,y},
udpoint=ToUnitDisc[z];
x=Re[udpoint];
y=Im[udpoint];
Text[text,{x,y}]
];

Circumcircle::usage="Circumcircle draws the circle passing through three non-collinear points.  The input can be six real numbers, three ordered pairs, or three complex numbers.";
Circumcircle[ax_,ay_,bx_,by_,cx_,cy_]:=Module[{c,r,output},
c=Circumcenter[ax,ay,bx,by,cx,cy];
r=Norm[Circumcenter[ax,ay,bx,by,cx,cy]-{ax,ay}];
output=Circle[c,r];
Return[output]]
Circumcircle[z1_,z2_,z3_]:=Circumcircle[Re[z1],Im[z1],Re[z2],Im[z2],Re[z3],Im[z3]]
Circumcircle[{ax_,ay_},{bx_,by_},{cx_,cy_}]:=Circumcircle[ax,ay,bx,by,cx,cy]

UnitCircle={Circle[{0,0},1]};
ShowP::usage="ShowP[List] generates a graphic object in the Poincar\[EAcute] Unit Disc from a list of graphics primitives.";
ShowP[list_]:=Show[Graphics[Append[list,UnitCircle]],PlotRange->{{-1.1,1.1},{-1.1,1.1}},AspectRatio->1]


(* ::Item:: *)
(*Mapping between the Upper Half Plane and the Unit Disc*)

ToUnitDisc::usage="ToUnitDisc[z] (\[DoubleStruckCapitalH]\[RightArrow]\[DoubleStruckCapitalD]) maps the point z in the Upper Half Plane to the corresponding point in the Unit Disc.";
ToUnitDisc[ComplexInfinity]=I;
ToUnitDisc[-I]=ComplexInfinity;
ToUnitDisc[z_]:= (z-I)/(-I z + 1)

ToCircle::usage="ToCircle[r] is to maintain backwards compatibility only.";
ToCircle[x_]:=ToUnitDisc[x]

ToUHP::usage="To UHP[z] (\[DoubleStruckCapitalD]\[RightArrow]\[DoubleStruckCapitalH]) or takes complex numbers or ordered pairs from the Unit Disk to the Upper Half Plane.";
ToUHP[z_]:=ApplyLFT[({
 {1/Sqrt[2], I/Sqrt[2]},
 {I/Sqrt[2], 1/Sqrt[2]}
}),z]
ToUHP[{a_,b_}]:=ToUHP[a+b*I]

ToR2::usage="ToR2[z] (\[DoubleStruckCapitalC]\[RightArrow]\[DoubleStruckCapitalR]\.b2) returns an ordered pair in \[DoubleStruckCapitalR]\.b2 from the real and imaginary parts of z.";
ToR2[z_]:={Re[z],Im[z]}


(* ::Item:: *)
(*Linear Fractional Transformations (Homographies)*)

ApplyLFT::usage="ApplyLFT[M,z] (\[DoubleStruckCapitalH]\[RightArrow]\[DoubleStruckCapitalH]) applies the linear fractional transformation corresponding to the 2x2 matrix M to the extended complex number z.  If z is a list (or list of lists) then ApplyLFT applies the linear fraction across the list.";
ApplyLFT[{{a_,b_},{c_,d_}},z_]:=N[(a*z+b)/(c*z+d)]
ApplyLFT[{{a_,b_},{c_,d_}},ComplexInfinity]:=N[a/c]
ApplyLFT[M_,l_List]:=Map[ApplyLFT[M,#]&,l]

FindLFT::usage="FindLFT[{z0,z1,zInf},{w0,w1,wInf}] Returns a 2x2 matrix in PSL(2,\[DoubleStruckCapitalR]) whose linear fractional transformation takes z1 to w1, z2 to w2, and z3 to w3.  If the transformation is not in PSL(2,\[DoubleStruckCapitalR]), it will return 'ERROR'.";
FindLFT[{z1_,z2_,z3_},{w1_,w2_,w3_}]:=Module[{A,B,M1,M2,Output},
A=CrossRatioMatrix[{z1,z2,z3}];
B=CrossRatioMatrix[{w1,w2,w3}];
M1=Adj[B].A;
M2=Simplify[M1/Sqrt[Det[M1]]];
If[Norm[Im[M2]]<.0002,Output=Re[M2]];
If[Norm[Re[M2]]<.0002,Output=Im[M2]];
If[Norm[Re[M2]]>.0002&&Norm[Im[M2]]>.0002,Output="ERROR"];
Return[Output]
]
(*To maintain backwards compatibility with previous versions.*)
LFTtaking[{z1_,z2_,z3_},{w1_,w2_,w3_}]:=FindLFT[{z1,z2,z3},{w1,w2,w3}]

(* ::Subitem:: *)
(*Commutators*)

Commutator::usage="Commutator[A,B] returns the commutator of the matrices A and B";
Commutator[A_,B_]:=A.B.Inverse[A].Inverse[B]

(* ::Subitem:: *)
(*Involutions*)

InvolutionFixing::usage="InvolutionFixing[z1,z2] Returns the 2x2 matrix representing an involution fixing
complex numbers z1, z2. It is normalized to have determinant -1.";
InvolutionFixing[z1_,z2_]:=1/(z1 - z2) { {z1 + z2, -2 z1 z2},{2, -(z1 + z2)}}
InvolutionFixing[z_,ComplexInfinity]:={{-1,2 z},{0,1}}
InvolutionFixing[ComplexInfinity,z_]:={{-1,2 z},{0,1}}

ReflectIn::usage="ReflectIn[{z1,z2},z] reflects the complex number z (upper half-plane coordinates)
in the Poincar\[EAcute] geodesic with endpoints z1, z2 (extended real numbers). This is the effect of applying
an anti-involution (anti-homography of order two) which fixes the chain orthogonal to the unit circle
(real axis in upper half-plane coordinates) whose endpoints on the unit circle correspond to
extended real numbers z1, z2.";
ReflectIn[{z1_,z2_},z_]:=ApplyLFT[InvolutionFixing[z1,z2],Conjugate[z]]

(* ::Subitem:: *)
(*Point-symmetries*)

SymmetryIn::usage="SymmetryIn[z] returns the involutary isometry of the Upper Half Plane fixing z.";
SymmetryIn[z_]:=Module[{x=Re[z],y=Im[z]},{{x/y,-(x^2/y)-y},{1/y,-(x/y)}}]

(* ::Subitem:: *)
(*Fixed Points*)

FixedPointSet::usage="FixedPointSet[M] Returns extended real numbers, representing
the attracting fixed point, then the repelling fixed point of the linear fractional
transformation represented by the 2x2 matrix M.";
FixedPointSet[M_]:=N[Map[Affine[#]&,Eigenvectors[M]]] /; Abs[Tr[M]]!=2;
FixedPointSet[M_]:=N[Affine[Eigenvectors[M][[1]]]] /;Abs[Tr[M]]==2;

(* ::Subitem:: *)
(*Cross Ratios*)

CrossRatio::usage="CrossRatio[a,b,c,d] returns the image of c under the collineation taking a to 0, b to 1 and d to \[Infinity].";
CrossRatio[a_,b_,c_,d_]:=(c-a)(b-d)/((c-d)(b-a))
CrossRatio[ComplexInfinity,b_,c_,d_]:=(b-d)/(c-d)
CrossRatio[a_,ComplexInfinity,c_,d_]:=(c-a)/(c-d) 
CrossRatio[a_,b_,ComplexInfinity,d_]:=(b-d)/(b-a)
CrossRatio[a_,b_,c_,ComplexInfinity]:=(c-a)/(b-a)

CrossRatioMatrix::usage="CrossRatioMatrix[{z0,z1,zInf}] Returns the matrix whose linear fractional transformation takes z0 to 0, z1 to 1, zInf to ComplexInfinity.";
CrossRatioMatrix[{z0_,z1_,zInf_}]:={{(z1-zInf),(z1-zInf)(-z0)},{(z1-z0),(z1-z0)(-zInf)}}
CrossRatioMatrix[{z0_,z1_,ComplexInfinity}]:={{1,-z0},{0,z1-z0}}
CrossRatioMatrix[{z0_,ComplexInfinity,zInf_}]:={{1,-z0},{1,-zInf}}
CrossRatioMatrix[{ComplexInfinity,z1_,zInf_}]:={{0,-(z1-zInf)},{-1,zInf}}


(* ::Subitem:: *)
(*Hyperbolic Elements*)

HypElement::usage="HypElement[{x1,x2}] Returns some Hyperbolic Element in PSL(2,R) with ideal fixed points x1 and x2.";
HypElement[{x1_,x2_}]:=Module[{A},A=LFTtaking[{x1,x2,Max[x1,x2]+1},{x1,x2,ComplexInfinity}] ;
A 
];
HypElement[{x1_,ComplexInfinity}]:=Module[{A},A=LFTtaking[{x1,ComplexInfinity,x1+1},{x1,ComplexInfinity,x1+2}];
A
];
HypElement[{ComplexInfinity,x2_}]:=Module[{A},A=LFTtaking[{ComplexInfinity,x2,x2+1},{ComplexInfinity,x2,x2+2}];
A
];


(* ::Item::Closed:: *)
(*Computing with Geodesics*)

(* ::Subitem:: *)
(*Distance*)

HyperbolicLength::usage="HyperbolicLength[z1,z2] (\[DoubleStruckCapitalH]\[RightArrow]\[DoubleStruckCapitalR]) Returns the length of the geodesic segment connecting z1 and z2.";
HyperbolicLength[z1_,z2_]:=Log[(Abs[z1-Conjugate[z2]]+Abs[z1-z2])/(Abs[z1-Conjugate[z2]]-Abs[z1-z2])];


(* ::Subitem:: *)
(*Endpoints and Midpints of Geodesics*)

MidpointEndpointsPGeodesic::usage="MidpointEndpointsPGeodesic[z1,z2] (\[DoubleStruckCapitalH]\[RightArrow]\[DoubleStruckCapitalH]) Returns the midpoint of the pair of endpoints of the Poincar\[EAcute] geodesic containing points z1 and z2 in the Upper Half Plane." ;
MidpointEndpointsPGeodesic[z1_,z2_]:=Re[(z1 + z2)/2 - I/2 (Im[z1+z2]/Re[z1-z2] (z1-z2))]

EndpointsPGeodesic::usage="EndpointsPGeodesic[z1,z2] (\[DoubleStruckCapitalH]\[RightArrow]\[DoubleStruckCapitalH]) Returns the endpoints of the Poincar\[EAcute] geodesic containing z1 and z2 in the Upper Half Plane.";
EndpointsPGeodesic[z1_,z2_]:=Module[{distance,midpoint=MidpointEndpointsPGeodesic[z1,z2]},
distance=Abs[z1-midpoint];
{midpoint-distance,midpoint+distance}]
EndpointsPGeodesic[z1_,z2_]:={Re[z1],ComplexInfinity}/;(Re[z1]==Re[z2])
EndpointsPGeodesic[z_,ComplexInfinity]:={Re[z],ComplexInfinity}
EndpointsPGeodesic[ComplexInfinity,z_]:={Re[z],ComplexInfinity}

(* ::Subitem:: *)
(*Intersections of Geodesics*)

IntersectionPoint::usage="IntersectionPoint[{z1,z2},{w1,w2}] (\[DoubleStruckCapitalH]\[RightArrow]\[DoubleStruckCapitalH]) Returns the complex number that is the point of intersection of the geodesics whose Ideal End Points are {z1,z2} and {w1,w2}.";
IntersectionPoint[{p1_,p2_},{q1_,q2_}]:=Module[{A,B,X,Y,fixed,intersectionpt},
A=InvolutionFixing[p1,p2];
B=InvolutionFixing[q1,q2];
X=HypElement[{p1,p2}];
Y=HypElement[{q1,q2}];
If[Tr[Commutator[X,Y]]<2,
fixed=FixedPointSet[A.B];
(*Print["These Geodesics Intersect"];*)
intersectionpt=If[Im[fixed[[1]]]>0,fixed[[1]],fixed[[2]]];
Return[intersectionpt];
];
If[Tr[Commutator[X,Y]]>2,
(*Print["These Geodesics do not intersect, they are ultraparallel"];*)
Return[{}];
];
If[Tr[Commutator[X,Y]]==2,
fixed=FixedPointSet[A.B];
(*Print["These geodesics are asymptotic and intersect in an ideal point"];*)
intersectionpt=fixed;
Return[ intersectionpt];
];
];
IntersectionPoint[{{p1_,p2_},{q1_,q2_}}]:=IntersectionPoint[{p1,p2},{q1,q2}];


(* ::Subitem:: *)
(*Orthogonal Geodesics*)

PerpendicularGeodesic::usage="PerpendicularGeodesic[{z1,z2},{w1,w2}] (\[DoubleStruckCapitalH]\[RightArrow]\[DoubleStruckCapitalH]) Returns the ideal endpoints of the geodesic perpendicular to the geodesics with endpoints {z1,z2} and {w1,w2}.";
PerpendicularGeodesic[{p1_,p2_},{q1_,q2_}]:=Module[{A,B,X,Y,out},
A=InvolutionFixing[p1,p2];
B=InvolutionFixing[q1,q2];
X=HypElement[{p1,p2}];
Y=HypElement[{q1,q2}];
If[Tr[Commutator[X,Y]]>2,out=FixedPointSet[A.B] ;Return[out];,Null];
If[Tr[Commutator[X,Y]]<=2,Print["There is no common perpendicular, these geodesics intersect"];,Null];
];
PerpendicularGeodesic[{{p1_,p2_},{q1_,q2_}}]:=PerpendicularGeodesic[{p1,p2},{q1,q2}];

PerpendicularBisector::usage="PerpendicularBisector[z1,z2] (\[DoubleStruckCapitalH]\[RightArrow]\[DoubleStruckCapitalH]) Returns the ideal endpoints of the line that is the perpendicular bisector of the line segment connecting z1 and z2.";
PerpendicularBisector[z_,w_]:=Module[{dist,endpt1,endpt2,T,S,midpoint,out},
dist=HyperbolicLength[z,w]/2;
{endpt1,endpt2}=EndpointsPGeodesic[z,w];
T={{Exp[dist/2],0},{0,Exp[-dist/2]}}.{{Cos[Pi/4],-Sin[Pi/4]},{Sin[Pi/4],Cos[Pi/4]}};
S=LFTtaking[{endpt1,z,endpt2},{0,I,ComplexInfinity}];
S=Chop[S/Sqrt[Det[S]] //N];
midpoint=ApplyLFT[Inverse[S].T.S,z];

(*S may have things in the incorrect order, we may have to switch the attracting and repelling fixed points*)
If[Abs[HyperbolicLength[midpoint,z]-HyperbolicLength[midpoint,w]]>10^-12,
S=LFTtaking[{endpt2,z,endpt1},{0,I,ComplexInfinity}];
S=Chop[S/Sqrt[Det[S]] //N];
];

out=Map[ApplyLFT[Inverse[S].T.S,#]&,{endpt1,endpt2}];
(*Replace numbers that are super duper large with ComplexInfinity*)
If[Abs[out[[1]]]>10^12,
out[[1]]=ComplexInfinity;
];
If[Abs[out[[2]]]>10^12,
out[[2]]=ComplexInfinity;
];
out=Chop[out];
Return[out];
];

(* ::Item:: *)
(*Other Useful Functions*)

Circumcenter::usage="Circumcenter returns the circumcenter of three noncollinear points.  When the input is six real numbers or three ordered pairs, Circumcenter outputs an ordered pair.  When the input is three complex numbers, Circumcenter outputs a complex number.";
Circumcenter[ax_,ay_,bx_,by_,cx_,cy_]:={
((ax^2+ay^2) (by-cy)+(bx^2+by^2) (cy-ay)+(cx^2+cy^2) (ay-by))/(2 (ax (by-cy)+bx (cy-ay)+cx (ay-by))),((ax^2+ay^2) (cx-bx)+(bx^2+by^2) (ax-cx)+(cx^2+cy^2) (bx-ax))/(2 (ax (by-cy)+bx (cy-ay)+cx (ay-by)))
}
Circumcenter[{ax_,ay_},{bx_,by_},{cx_,cy_}]:=Circumcenter[ax,ay,bx,by,cx,cy]
Circumcenter[z1_,z2_,z3_]:=Module[{ztemp},
ztemp=Circumcenter[Re[z1],Im[z1],Re[z2],Im[z2],Re[z3],Im[z3]];
Return[ztemp[[1]]+ztemp[[2]] I]
]

Affine::usage="Affine[List] Returns the inhomogeneous coordinates of the list of homogeneous coordinates of a point in projective space.";
Affine[List_]:=First[List]/Last[List]
Affine[{z_,0}]:=ComplexInfinity

Projective::usage="Projective[z] Returns the homogeneous coordinates of a point in P^1 corresponding to the extended complex number z.";
Projective[z_]:={z,1}
Projective[ComplexInfinity]:={1,0}

Adj::usage="Adj[M] Returns the adjugate of the 2x2 matrix M.  Adj[M]=Det[M]Inverse[M].";
Adj[{{a_,b_},{c_,d_}}]:={{d,-b},{-c,a}}