(*
	(c) 2010 Alex Kontorovich
*)

(* Initial parameters *)

(* Keep only geodesics with endpoints in (windowMin,windowMax) which reach at least heightMin into the upper half plane *)
windowMin=-3;
windowMax=3;
heightMin=0.1;

shift=1;(* shift=1 for SL(2,Z) and 2 for Gamma[2] *)

Tmax=8; (* Will make a ball in SL(2,Z) of size Tmax *)

(* Code *)

xmin=-shift/2;
xmax=shift/2;
sph={{0,0,1}};
Dom=RegionPlot[(((x-sph[[1,1]]))^2)+(((y-sph[[1,2]]))^2)>(sph[[1,3]])^2&&y>0&&-shift/2<x<shift/2,{x,-shift,shift},{y,-.3,100},AspectRatio->Automatic,FrameTicks->None];
toFundDom[v_]:=
Module[{vFD=v,checkMore=True,j},
While[checkMore,
checkMore=False;
If[Abs[vFD[[1]]]>shift/2,
vFD[[1]]=vFD[[1]]-Round[vFD[[1]]/shift]shift;
checkMore=True;
,
If[(vFD[[1]]-sph[[1,1]])^2+(vFD[[2]]-sph[[1,2]])^2<sph[[1,3]]^2,
vFD={sph[[1,1]],sph[[1,2]]}+sph[[1,3]]^2(vFD-{sph[[1,1]],sph[[1,2]]})/Norm[(vFD-{sph[[1,1]],sph[[1,2]]})]^2;
vFD={-vFD[[1]],vFD[[2]]};

checkMore=True;
];
];
];
vFD
];

(* Some elements of G=SL(2,R), G=NAK, etc *)
a[s_]:=MatrixExp[s/2{{1,0},{0,-1}}];
n[x_]:=MatrixExp[x{{0,1},{0,0}}];
k[gt_]:=MatrixExp[gt{{0,1},{-1,0}}];

GAct[g_,z_]:=(g[[1,1]]z +g[[1,2]])/(g[[2,1]]z +g[[2,2]]);
CtoR2[z_]:={Re[z],Im[z]};

geod[x0_,s_]:=GAct[x0.a[s],I];

geodOrbit[mat_]:=(* makes pics of gammas orbit of geodesic; input: Vis+ mat is one endpoint, Vis- mat is another ; output into "plots" *)
Module[{j,p1,p2,pic,m},
For[j=1,j≤Length[gammas],j++,
m=gammas[[j]]. mat;
If[m[[2,1]]≠0,
p1=m[[1,1]]/(m[[2,1]]);
,
p1=I;
(*Limit[GAct[m. a[s],I],s->Infinity];*)
];
If[m[[2,2]]≠0,
p2=m[[1,2]]/(m[[2,2]]);
,
p2=I; (*Limit[GAct[m. a[s],I ],s->-Infinity,Assumptions->s<0];*)
];

(* Are we within visible range? *)
If[ 
((windowMin<p1<windowMax)
||
(windowMin<p2<windowMax))
&&
(Abs[p1-p2]>heightMin)
,

If[
(p1∈Reals) && (p2∈Reals),
(* Draw circles? *)
pic=
Circle[{(p1+p2)/2,0},Abs[(p1-p2)/2],{0,Pi}];
,
(* Draw lines *)
If[
p1∈Reals,
pic=
Line[{{ p1  ,0},{  p1 ,10}}]
;
,
pic=
Line[{{ p2 ,0},{  p2,10}}]
;
];
]; (* end if *)
AppendTo[plots,pic];
];
];
];

MakeGammas[Tmax_]:=(* Make SL(2,Z) (or Gamma[2] if shift=2) ball of radius T, output into "gammas" *)
Module[{aa,bb,cc,dd},
For[aa=1,aa≤Tmax,aa=aa+shift,
For[dd=-Tmax shift-1,dd≤shift Tmax+1,dd=dd+shift,
If[aa dd-1≠0,
bcs=Divisors[aa dd -1];
For[bb=1,bb≤Length[bcs],bb++,
cc=(aa dd -1 )/bcs[[bb]];
If[bcs[[bb]]<=Tmax && cc<=Tmax && Mod[bcs[[bb]],shift]==0&& Mod[cc,shift]==0,
AppendTo[gammas,
{{aa,bcs[[bb]]},{cc,dd}}
];
AppendTo[gammas,
{{aa,cc},{bcs[[bb]],dd}}
];
AppendTo[gammas,
{{aa,-bcs[[bb]]},{-cc,dd}}
];
AppendTo[gammas,
{{aa,-cc},{-bcs[[bb]],dd}}
];
];
];
,
(* else unipotent *)
For[bb=shift,bb≤5Tmax shift,bb=bb+shift,
AppendTo[gammas,
{{aa,bb},{0,dd}}
];
AppendTo[gammas,
{{aa,0},{bb,dd}}
];
AppendTo[gammas,
{{aa,-bb},{0,dd}}
];
AppendTo[gammas,
{{aa,0},{-bb,dd}}
];
];
];
];
];
(* aa=0 *)
If[shift==1, (* Not elements of Gamma[2] *)
For[dd=shift,dd≤5 shift Tmax,dd=dd+shift,
AppendTo[gammas,
{{0,1},{-1,dd}}
];
AppendTo[gammas,
{{0,-1},{1,dd}}
];
];
];
];


(* Now run some things: *)


(* Make a ball in SL(2,Z) of size Tmax *)
gammas={{{1,0},{0,1}}};
MakeGammas[Tmax];
Print["Length of gammas is ",Length[gammas]]; (* Make sure not trying to tile too make guys! *)


(* Make a tiling of Gamma[1] *)
plots={};
Quiet[geodOrbit[{{1,1},{2,0}}]];
plotsSL2Z=plots;

(* Tiling of golden mean geodesic *)
plots={};
Quiet[geodOrbit[{{1+Sqrt[5],1-Sqrt[5]},{2,2}}]];
plotsGoldenMean=plots;

(* Shows golden mean geodesic against SL(2,Z) tiling *)
Show[Dom,Graphics[{Thin,plotsSL2Z}],
Graphics[{Red,Thick,plotsGoldenMean}],
PlotRange->{{-.6,1.2},{-.1,2.3}},AxesOrigin->{0,0}]

(* Tiling of Low-Lying geod *)

(* Parameters *)

sFrom=-.52; (* start and end of geodesic flow a[s] *)
sTo=6.2;
sStep=.006;

plotRan={{-1.1,1.1},{-.1,2}}; (* plot Range *)

ga=FromContinuedFraction[{0,2,1,1,2}]; (* Use small numbers here to get low-lying geodesic *)

Print["alpha is ",ga];

(* the code *)

gg={{Numerator[ga],-Denominator[ga]},{Denominator[ga],Numerator[ga]}}; (* this is an (un-normalized) element of the tangent bundle, fixing I *)

plots={};
Quiet[geodOrbit[gg]];
plotsLowLying=plots;

FixedPic=
Show[Dom,
Graphics[plotsSL2Z],
(* Graphics[{Red,Thick,plotsGoldenMean}], *)
Graphics[{Blue,plotsLowLying}],

(* geodesic flow in H at gg *)
ParametricPlot[
CtoR2[GAct[
gg.
a[s],I]]
,
{s,sFrom,sTo},AxesOrigin->{0,0},PlotStyle->Thick]
,

(* geodesic flow in Gamma \ H at gg *)
ListPlot[
Table[
toFundDom[
CtoR2[GAct[
gg.
a[s],I]]],{s,sFrom,sTo,sStep}],AxesOrigin->{0,0},PlotStyle->Thick]
,PlotRange->plotRan,AxesOrigin->{0,0}
]

(* Now make a movie of the geodesic flow, showing the point in the tangent bundle *)

(* Parameters *)
dt=0.0001; (* tiny step for making an arrow *)
sStepT=0.05; (* frequency of sampling geodesic flow *)
lead=-.4; (* controls distance of arrowhead from point *)

tab=Table[
Show[
FixedPic
,

(* make a point in H *)
Graphics[{Blue,PointSize[ Large],Point[CtoR2[GAct[
gg.a[t],I]]]
}]
,
(* make its arrow *)
Graphics[{Blue,Arrow[
{CtoR2[GAct[
gg.a[t],I]],CtoR2[GAct[
gg.a[t+dt],I]]}
,{0,lead}
]
}]
,

(* make a point in Gamma \ H and its arrow *)
Graphics[{Blue,PointSize[ Large],Point[
toFundDom[CtoR2[GAct[
gg.a[t],I]]]]
}]
,
Graphics[{Blue,Arrow[
{toFundDom[CtoR2[GAct[
gg.a[t],I]]],
toFundDom[CtoR2[GAct[
gg.a[t+dt],I]]]
},{0,lead}
]
}]
,
PlotRange->plotRan,AxesOrigin->{0,0}]
,{t,sFrom,sTo,sStepT}];

(* Check that it looks good *)
tab[[15]]

(* Check how many frames there are *)
Length[tab ]

(* Make a movie out of tab; This will take a while *)
Export["movLowLying.swf",tab];
Beep[];