(* ::Package:: *)

(* ::Title:: *)
(*Turbo Tessellate Version 1.1*)

(* ::Subsubtitle:: *)
(*Ammar Husain*)
(*Ravi Joseph*)
(*Brian Baird*)
(**)
(*Experimental Geometry Lab*)
(*University of Maryland*)
(*January 2011*)
(**)
(*http://egl.math.umd.edu/*)

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

(* ::Text:: *)
(*The purpose of this function is to perform tessellations of the once-punctured torus (regular and flared) and display the results on the Poincar\[EAcute] unit disc model.  Turbo Tessellate uses a depth-first recursion method that applies only those transformations that will provide visible results.  This allows the user to create highly detailed images in a short amount of time, and even tessellate in real-time as the input variables are changed.*)

(* ::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[]];*)
(*<< TurboTessellate_v1.1.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:: *)
(*Tessellating Hyperbolic Space with Ideal Quadrilaterals*)

(* ::Subsection:: *)
(*Input*)

(* ::Text:: *)
(*The required input for a basic tessellation is:*)
(**)
(*TurboTessellate[A,B]*)
(**)
(*Where A and B are two matrices that form a subgroup of PSL(2, \[DoubleStruckCapitalR]).  When the commutator of A and B is parabolic, a complete tiling of the hyperbolic plane will occur.  When the commutator is not parabolic, an incomplete tiling of the hyperbolic plane will occur.*)
(**)

(* ::Subsection:: *)
(*Threshold Parameters*)

(* ::Text:: *)
(*TurboTessellate allows for two additional parameters, Threshold and MaxWordLength.  Threshold sets the minimum size of the tessellation, while MaxWordLength specifies the maximum depth of the tessellation.*)
(**)
(*When using a specified Threshold and MaxWordLength, the required input is:*)
(**)
(*TurboTessellate[A,B,Threshold,MaxWordLength]*)
(**)
(*The minimum recommended value for Threshold is .01.  The range for MaxWordLength is 1 to 100.  By default these parameters are set to .1 and 50, respectively.*)
(**)

(* ::Section:: *)
(*Tessellating Hyperbolic Space with Fundamental Domains That are Not Ideal Quadrilaterals*)

(* ::Text:: *)
(*When the commutator of A and B is hyperbolic, tessellating hyperbolic space with an ideal quadrilateral will result in an incomplete tiling.  However, when the commutator is hyperbolic, a complete tiling can be achieved by using a region in hyperbolic space bounded by four disjoint sides with eight distinct endpoints on the boundary at infinity.  Such a region is not compact and will meet the boundary at infinity in four arcs.  This region will be a fundamental domain for a once-punctured torus with a flared end (i.e. infinite area).*)

(* ::Subsection:: *)
(*Input*)

(* ::Text:: *)
(*In this case, TurboTessellate requires an additional input: FirstPoint.  Each value of FirstPoint (range -1 to 1) corresponds to a different tiling of the flared torus for the same marked group, with the values -1 and 1 corresponding to incomplete tilings.*)
(**)
(*The required input for this tessellation is:*)
(**)
(*TurboTessellate[A, B, FirstPoint]*)
(**)

(* ::Subsection:: *)
(*Threshold Parameters*)

(* ::Text:: *)
(*Threshold and MaxWordLength can also be set using this method.  The required input is:*)
(**)
(*TurboTessellate[A, B, FirstPoint, Threshold, MaxWordLength]*)
(**)
(*The minimum value for Threshold is .01.  The range for MaxWordLength is 1 to 100.  By default these parameters are set to .1 and 50, respectively.*)

(* ::Section:: *)
(*Optimizing TurboTessellate Based on Desired Output*)

(* ::Subsection:: *)
(*Dynamic Rendering*)

(* ::Text:: *)
(*The rendering of a tessellation can be a computationally intensive task in Mathematica.  However, the method used by TurboTessellate allows for detailed tessellations to be drawn in real time using the Dynamic or Manipulate commands.*)
(**)
(*When using TurboTessellate to dynamically render a tessellation, keep in mind the total number of objects drawn will limit the performance.  For this reason, it is generally recommended to set the Threshold parameter above .2 and MaxWordLength below 25 when using TurboTessellate dynamically.*)

(* ::Subsection:: *)
(*Static Renders*)

(* ::Text:: *)
(*To create a highly detailed tessellation, set Threshold to .01 and MaxWordLength to 100.*)

(* ::Subsection:: *)
(*Adjusting Colors and Image Size*)

(* ::Text:: *)
(*TurboTessellate has two additional parameters to allow the user to configure the colors of the geodesics, as well as the resulting image size.  ColorsList determines the values relating to Mathematica's RGBColor graphics directive, while Size determines the pixel dimensions of the resulting graphics object.  Both parameters must be specified at the same time along with Threshold and MaxWordLength, such as:*)
(**)
(*TurboTessellate[A, B, Threshold, MaxWordLength, ColorList, Size]*)
(**)
(*or*)
(**)
(*TurboTessellate[A, B, FirstPoint, Threshold, MaxWordLength, ColorList, Size]*)
(**)
(**)
(*ColorList should be formatted in the following manner: {{R1,G1,B1},{R2,G2,B2},{R3,G3,B3}}.  Size should be an integer.  "DefaultColors" can be substituted for ColorList if the user only wishes to change the image size.*)

(* ::Section::Closed:: *)
(*Code*)

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

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


(* ::Item:: *)
(*Function Interface*)

TurboTessellate::usage="TurboTessellate[A,B,Threshold,MaxWordLength] Draws the tessellation of the quadrilateral given by the LFTs A and B.  Threshold is the minimum chord length required to continue applying transformations into a given domain, MaxWordLength is the maximum word length allowed during the tessellation.  If unspecificied, these values will default to .1 and 50, respectively";

TurboTessellate[A_,B_,Threshold_,MaxWordLength_,ColorsList_,Size_]:=Module[{qAB,qAb,qab,qaB,Output},
{AGeodesics={};
BGeodesics={};
Diagonals={};

(*Find the trace of the commutator for cutoff purposes*)
CommutatorTrace=Tr[TTCommutator[A,B]];

(*Find the vertices of the ideal quadrilateral to be used*)
qAB=ParabolicFP[TTCommutator[A,B]];
qAb=TTApplyLFT[Inverse[B],qAB];
qab=TTApplyLFT[Inverse[A].Inverse[B],qAB];
qaB=TTApplyLFT[Inverse[A],qAB];

Tessellation[0,IdentityMatrix[2],N[A],N[B],qAb, qab,qaB, qAB,Threshold,MaxWordLength]};

Output=Graphics[Join[{RGBColor[ColorsList[[3]]]},AGeodesics,{RGBColor[ColorsList[[1]]]},BGeodesics,{Opacity[.5],RGBColor[ColorsList[[2]]]},Diagonals,{Opacity[1],Thick,Gray,Circle[{0,0},1],Thick,White,Circle[{0,0},1.0039]}],ImageSize->{Size,Size}]
];

(*Generate a tiling of the plane with two generators A and B with a hyperbolic commutator. 
FirstPoint is a value between -1 and 1. Each value between -1 and 1 will correspond to a different tiling for the same marked group,
with -1 and 1 corresponding to incomplete structures.*)
TurboTessellate[A_,B_,FirstPoint_,Threshold_,MaxWordLength_,ColorsList_,Size_]:=Module[{PointList,qa1, qa2,qA1, qA2,qb1, qb2,qB1 ,qB2,Output},
{AGeodesics={};
BGeodesics={};

(*Find the trace of the commutator for cutoff purposes*)
CommutatorTrace=Tr[TTCommutator[A,B]];

(*Calculate endpoints of the fundamental domain to be used*)
PointList= ReturnEndPointsFlared[A,B,FirstPoint];
qa1=PointList[[1]];
qa2=PointList[[2]];
qA1=PointList[[3]];
qA2=PointList[[4]];
qb1=PointList[[5]];
qb2=PointList[[6]];
qB1=PointList[[7]];
qB2=PointList[[8]];

If[Abs[FirstPoint]==1,Tessellation[0,IdentityMatrix[2],N[A],N[B],qA2,qa2,qa1,qA1,Threshold,MaxWordLength],
TessellationFlared[0,IdentityMatrix[2],N[A],N[B],qa1, qa2,qA1, qA2,qb1, qb2,qB1 ,qB2,Threshold,MaxWordLength]];

};

Output=Graphics[Join[{RGBColor[ColorsList[[1]]]},AGeodesics,{RGBColor[ColorsList[[2]]]},BGeodesics,{Thick,Gray,Circle[{0,0},1],Thick,White,Circle[{0,0},1.0039]}],ImageSize->{Size,Size}]
];

(*If additional variables are not specified, defaults are used*)
TurboTessellate[A_,B_]:=TurboTessellate[A,B,DefaultThreshold,DefaultWordLength,DefaultColors,DefaultSize]
TurboTessellate[A_,B_,Threshold_,MaxWordLength_]:=TurboTessellate[A,B,Threshold,MaxWordLength,DefaultColors,DefaultSize]
TurboTessellate[A_,B_,FirstPoint_]:=TurboTessellate[A,B,FirstPoint,DefaultThreshold,DefaultWordLength,DefaultColors,DefaultSize]
TurboTessellate[A_,B_,FirstPoint_,Threshold_,MaxWordLength_]:=TurboTessellate[A,B,FirstPoint,Threshold,MaxWordLength,DefaultColors,DefaultSize]


(* ::Item:: *)
(*Basic Functions Needed For Routine*)

(*Set defaults for function*)
DefaultThreshold=.1;
DefaultWordLength=50;
DefaultColors={{0,.8,0},{0,0,1},{1,0,0}};
DefaultSize=600;

(*Only circular arcs with an Euclidean radius above this threshold will be drawn.*)
DefaultTTSegmentThreshold=.003;

ChordLenth::usage="The chord length between two points on the boundary is used to determine if a given transformation will be applied.";
ChordLength[q1_,q2_]:=Norm[TTtoUnitDisc[q1]-TTtoUnitDisc[q2]]

(*Finds the fixed point of the commutator*)
ParabolicFP[M_]:=If[Chop[Part[Part[M,2],1]]==0,ComplexInfinity,(Part[Part[M,1],1]-Part[Part[M,2],2])/(2*Part[Part[M,2],1])]

(*Returns the ideal fixed points of the commutator.*)
IdealFPs[M_] := If[Tr[M]!=2,TTFixedPointSet[M],{ParabolicFP[M],ParabolicFP[M]}]


(* ::Item:: *)
(*Normal Routine*)

Tessellation[lastDrawn_,CurrentMatrix_, A_,B_, qAb_, qab_,qaB_, qAB_,Threshold_,MaxWordLength_]:=Module[
{ChildrenToKeep={},CurrentWordLength},{
CurrentWordLength=MaxWordLength-1;

ChildrenToKeep=KeepWorthy[lastDrawn, qAb, qab,qaB, qAB,CurrentMatrix,Threshold];

If[CurrentWordLength<=0,ChildrenToKeep={}];

If[MemberQ[ChildrenToKeep,1],
Tessellation[1,CurrentMatrix.A,A,B, qAb, qab,qaB, qAB,Threshold,CurrentWordLength]];
If[ MemberQ[ChildrenToKeep,2],
Tessellation[2,CurrentMatrix.B,A,B, qAb, qab,qaB, qAB,Threshold,CurrentWordLength]];
If[MemberQ[ChildrenToKeep,-1],
Tessellation[-1,CurrentMatrix.Inverse[A],A,B, qAb, qab,qaB, qAB,Threshold,CurrentWordLength]];
If[MemberQ[ChildrenToKeep,-2],
Tessellation[-2,CurrentMatrix.Inverse[B],A,B, qAb, qab,qaB, qAB,Threshold,CurrentWordLength]];
}]

KeepWorthy::usage="KeepWorthy tests the Euclidean chord length of the ideal points of hyperbolic geodesics to determine if they are larger than the set threshold.  If they are not, no further transformations will procede into that domain.";
KeepWorthy[lastDrawn_,qAb_, qab_,qaB_, qAB_,CurrentMatrix_,Threshold_]:=
Module[{q1,q2,q3,q4,r12,r23,r34,r41,r13,r24,Output1={},Output2={}},
q1=TTApplyLFT[CurrentMatrix,qAb];
q2=TTApplyLFT[CurrentMatrix,qab];
q3=TTApplyLFT[CurrentMatrix,qaB];
q4=TTApplyLFT[CurrentMatrix,qAB];

r12=ChordLength[q1,q2];
r23=ChordLength[q2,q3];
r34=ChordLength[q3,q4];
r41=ChordLength[q4,q1];
r13=ChordLength[q1,q3];
r24=ChordLength[q2,q4];

GetQuadrilateral[q1,q2,q3,q4,lastDrawn];

If[r12>Threshold&&lastDrawn!=2,Output2=Append[Output2,-2]];
If[r23>Threshold&&lastDrawn!=1,Output2=Append[Output2,-1]];
If[r34>Threshold&&lastDrawn!=-2,Output2=Append[Output2,2]];
If[r41>Threshold&&lastDrawn!=-1,Output2=Append[Output2,1]];

If[CommutatorTrace!=-2&&(r13<Threshold/2||r24<Threshold/2),Output2={}];

Return[Output2];

];

GetQuadrilateral[qAb_, qab_,qaB_, qAB_,lastDrawn_]:=Module[{geoA={},geoB={},geoa={},geob={},diagonal},
If[lastDrawn!=-1,
geoA={TTSegment[{qAB,qAb},DefaultTTSegmentThreshold]}];
If[lastDrawn!=1,
geoa={TTSegment[{qaB,qab},DefaultTTSegmentThreshold]}];
If[lastDrawn!=2,
geob={TTSegment[{qAb,qab},DefaultTTSegmentThreshold]}];
If[lastDrawn!=-2,
geoB={TTSegment[{qAB,qaB},DefaultTTSegmentThreshold]}];
diagonal={TTSegment[{qAB,qab},DefaultTTSegmentThreshold]};

AGeodesics=Join[AGeodesics,geoA,geoa];
BGeodesics=Join[BGeodesics,geoB,geob];
Diagonals=Append[Diagonals,diagonal];
];


(* ::Item:: *)
(*Hyperbolic Commutator Routine*)

TessellationFlared[lastDrawn_,currentMatrix_, A_,B_, qa1_, qa2_,qA1_, qA2_,qb1_, qb2_,qB1_, qB2_,Threshold_,MaxWordLength_]:=Module[

{ChildrenToKeep={},CurrentWordLength},{
CurrentWordLength=MaxWordLength-1;

ChildrenToKeep=KeepWorthyFlared[lastDrawn, qa1, qa2,qA1, qA2,qb1, qb2,qB1, qB2,currentMatrix,Threshold];

If[CurrentWordLength<=0,ChildrenToKeep={}];

If[MemberQ[ChildrenToKeep,1],
TessellationFlared[1,currentMatrix.A,A,B, qa1, qa2,qA1, qA2,qb1, qb2,qB1, qB2,Threshold,CurrentWordLength]];
If[ MemberQ[ChildrenToKeep,2],
TessellationFlared[2,currentMatrix.B,A,B,qa1, qa2,qA1, qA2,qb1, qb2,qB1, qB2,Threshold,CurrentWordLength]];
If[MemberQ[ChildrenToKeep,-1],
TessellationFlared[-1,currentMatrix.Inverse[A],A,B, qa1, qa2,qA1, qA2,qb1, qb2,qB1, qB2,Threshold,CurrentWordLength]];
If[MemberQ[ChildrenToKeep,-2],
TessellationFlared[-2,currentMatrix.Inverse[B],A,B, qa1, qa2,qA1, qA2,qb1, qb2,qB1, qB2,Threshold,CurrentWordLength]];
}]

KeepWorthyFlared[lastDrawn_,qa1_, qa2_,qA1_, qA2_,qb1_, qb2_,qB1_, qB2_,CurrentMatrix_,Threshold_]:=
Module[{q1,q2,q3,q4,q5,q6,q7,q8,r1,r2,r3,r4,Output1={},Output2={}},

q1=TTApplyLFT[CurrentMatrix,qb2];
q2=TTApplyLFT[CurrentMatrix,qb1];
q3=TTApplyLFT[CurrentMatrix,qA1];
q4=TTApplyLFT[CurrentMatrix,qA2];
q5=TTApplyLFT[CurrentMatrix,qB1];
q6=TTApplyLFT[CurrentMatrix,qB2];
q7=TTApplyLFT[CurrentMatrix,qa1];
q8=TTApplyLFT[CurrentMatrix,qa2];

r1=ChordLength[q1,q2];
r2=ChordLength[q7,q8];
r3=ChordLength[q5,q6];
r4=ChordLength[q3,q4];

GetQuadrilateralFlared[q3,q4,q7,q8,q1,q2,q5,q6,lastDrawn];

If[r1>Threshold&&lastDrawn!=2,Output2=Append[Output2,-2]];
If[r2>Threshold&&lastDrawn!=1,Output2=Append[Output2,-1]];
If[r3>Threshold&&lastDrawn!=-2,Output2=Append[Output2,2]];
If[r4>Threshold&&lastDrawn!=-1,Output2=Append[Output2,1]];

Return[Output2];

];

GetQuadrilateralFlared[qa1_, qa2_,qA1_, qA2_,qb1_, qb2_,qB1_, qB2_,lastDrawn_]:=

Module[{geoA={},geoB={},geoa={},geob={}},
If[lastDrawn!=-1,
geoA={TTSegment[{qa1,qa2},DefaultTTSegmentThreshold]}];
If[lastDrawn!=1,
geoa={TTSegment[{qA1,qA2},DefaultTTSegmentThreshold]}];
If[lastDrawn!=2,
geob={TTSegment[{qb1,qb2},DefaultTTSegmentThreshold]}];
If[lastDrawn!=-2,
geoB={TTSegment[{qB1,qB2},DefaultTTSegmentThreshold]}];

AGeodesics=Join[AGeodesics,geoA,geoa];
BGeodesics=Join[BGeodesics,geoB,geob];
];

(*Returns the eight endpoints corresponding to a fundamental domain.*)
ReturnEndPointsFlared[A_,B_,p_]:=Chop[{
TTApplyLFT[Inverse[A],Part[FindFourPointsFlared[A,B,p],3]],
TTApplyLFT[Inverse[B].Inverse[A],Part[FindFourPointsFlared[A,B,p],1]],
Part[FindFourPointsFlared[A,B,p],3],
TTApplyLFT[A.Inverse[B].Inverse[A],Part[FindFourPointsFlared[A,B,p],1]],
TTApplyLFT[Inverse[B],Part[FindFourPointsFlared[A,B,p],4]],
TTApplyLFT[Inverse[B].Inverse[A],Part[FindFourPointsFlared[A,B,p],2]],
Part[FindFourPointsFlared[A,B,p],4],
TTApplyLFT[Inverse[A],Part[FindFourPointsFlared[A,B,p],2]]
}]

FindFourPointsFlared[A_,B_,x_] := 
(*FourPoints returns the first four points needed to construct a fundamental domain for a pair of generators A,B.*)
Module[
{FirstPoint,ImageFirstPoint,ImageAngle,PointonSide,HalfwayPointComm,HalfwayAngleComm,FixedPoint1 = Part[IdealFPs[B.A.Inverse[B].Inverse[A]],1],
FixedPoint2 = Part[IdealFPs[B.A.Inverse[B].Inverse[A]],2]},

(*Find half the angle between the two fixed points of the commutator.*)
HalfwayAngleComm = .5*Arg[ TTtoUnitDisc[FixedPoint2]/TTtoUnitDisc[FixedPoint1]];

(*Find the point on the unit disc halfway between the two fixed points of the commutator*)
HalfwayPointComm = Exp[I*HalfwayAngleComm]*TTtoUnitDisc[FixedPoint1];

(*The commutator geodesic splits the hyperbolic plane into two half-planes. Check to see that the halfway point is 'inside' the half plane enclosed by the commutator geodesic. If not, return the point directly opposite from the halfway point on the unit disc (the new point will now be 'inside').*)
PointonSide = If[N[Abs[Arg[TTtoUnitDisc[Part[IdealFPs[A],1]]/HalfwayPointComm]]]<=HalfwayAngleComm,Exp[-I*(Pi-HalfwayAngleComm)]*HalfwayPointComm,HalfwayPointComm];

(*Find the first point of the four points by displacing the halfway point over by the parameter x.*)
FirstPoint = Exp[I*x*HalfwayAngleComm]*PointonSide;

(*The remaining points are found by taking the angle between the first point and its image under the commutator baBA of A,B and displacing the first point by 0, .25, .5, .75 times the angle. The function returns the final four points.*)
ImageFirstPoint = TTtoUnitDisc[TTApplyLFT[B.A.Inverse[B].Inverse[A],TTtoUHP[FirstPoint]]];
ImageAngle = Arg[ImageFirstPoint/FirstPoint];
N[Map[TTtoUHP,{FirstPoint,Exp[I*.25*ImageAngle]*FirstPoint,Exp[I*.5*ImageAngle]*FirstPoint,Exp[I*.75*ImageAngle]*FirstPoint}]]
]


(* ::Item:: *)
(*Shared Functions with HyperbolicToolbox*)

TTApplyLFT::usage="TTApplyLFT[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 TTApplyLFT applies the linear fraction across the list.";
TTApplyLFT[{{a_,b_},{c_,d_}},z_]:=N[(a*z+b)/(c*z+d)]
TTApplyLFT[{{a_,b_},{c_,d_}},ComplexInfinity]:=N[a/c]
TTApplyLFT[M_,l_List]:=Map[TTApplyLFT[M,#]&,l]

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

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

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

TTSegment::usage="TTSegment[{z1,z2}] (\[DoubleStruckCapitalH]\[RightArrow]\[DoubleStruckCapitalD], Graphics) Given two points from the Upper Half Plane, TTSegment draws hyperbolic lines on the disc model.  An optional third argument specifies the minimum draw size.";
TTSegment[{z1_,z2_}]:=TTSegment[{z1,z2},.003]
TTSegment[{z1_,z2_},Threshold_]:=Module[{Disc1,Disc2,x1,y1,x2,y2,Output},
(*Move the coordinates to the disc model.*)
Disc1=TTtoUnitDisc[z1];
Disc2=TTtoUnitDisc[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=TTDrawArc[x1,y1,x2,y2,Threshold]];
Return[Output]
]
TTDrawArc[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]
]

TTFixedPointSet::usage="TTFixedPointSet[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.";
TTFixedPointSet[M_]:=N[Map[TTAffine[#]&,Eigenvectors[M]]] /; Abs[Tr[M]]!=2;
TTFixedPointSet[M_]:=N[TTAffine[Eigenvectors[M][[1]]]] /;Abs[Tr[M]]==2;

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