(*
	(c) 2008 Alex Kontorovich
*)

(* Parameters *)

xmin=-5; (* Ranges of plots *)
xmax=8.5;
actmin=-5;
actmax=8.5;
xstep=.005;

(* These spheres {x,y,rad} reflections will be the group (really GL(2,R)... Det=-1...) *)
sph={
{0,0,1/2},
{59/20,0,47/20},
{13/2,0,11/10}
};


(* The Code *)

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

(* Running the program *)
height=Table[x^2,{x,3,.1,-.1}]; (* these are the heights of the horocycle *)
pic=height; (* Stupid init *)

(* Make frames *)
For[j=1,j≤Length[height],j++,
list=ListPlot[Table[toFundDom[{x,height[[j]]}],{x,actmin,actmax,xstep}]];
pic[[j]]=Show[Dom,list(*,Graphics[{Darker[Blue],Thickness[.01],Line[{{xmin,height[[j]]},{actmin,height[[j]]}}]}]
,Graphics[{Darker[Blue],Thickness[.01],Line[{{actmax,height[[j]]},{xmax,height[[j]]}}]}]*)];
];

(* Check a frame *)
pic[[30]]

(* This will take a while *)
Export["mov.swf",pic];
Beep[]

(* Now do the same thing, just in the Disk model *)

(* The parameters *)
actmin=-35;
actmax=35;
xstep=.01;


(* The Code *)

HtoD[z_]:=(z-I)/(z+I)//Simplify;
CtoR2[w_]:={Re[w],Im[w]};
R2toC[p_]:=p[[1]]+I p[[2]];
p1p=CtoR2[HtoD[sph[[1,1]]+sph[[1,3]]]];
p1m=CtoR2[HtoD[sph[[1,1]]-sph[[1,3]]]];
p2p=CtoR2[HtoD[sph[[2,1]]+sph[[2,3]]]];
p2m=CtoR2[HtoD[sph[[2,1]]-sph[[2,3]]]];
p3p=CtoR2[HtoD[sph[[3,1]]+sph[[3,3]]]];
p3m=CtoR2[HtoD[sph[[3,1]]-sph[[3,3]]]];
p1c={x,y}/.Solve[
{
({x,y}-p1p).p1p==0,
({x,y}-p1m).p1m==0
}
,{x,y}
][[1]];
p2c={x,y}/.Solve[
{
({x,y}-p2p).p2p==0,
({x,y}-p2m).p2m==0
}
,{x,y}
][[1]];
p3c={x,y}/.Solve[
{
({x,y}-p3p).p3p==0,
({x,y}-p3m).p3m==0
}
,{x,y}
][[1]];
sphD={
{p1c[[1]],p1c[[2]],Norm[p1c-p1p]},
{p2c[[1]],p2c[[2]],Norm[p2c-p2p]},
{p3c[[1]],p3c[[2]],Norm[p3c-p3p]}
};
FDom=RegionPlot[x^2+y^2<1
&&(((x-sphD[[1,1]]))^2)+(((y-sphD[[1,2]]))^2)>(((sphD[[1,3]]))^2)
&&(((x-sphD[[2,1]]))^2)+(((y-sphD[[2,2]]))^2)>(((sphD[[2,3]]))^2)
&&(((x-sphD[[3,1]]))^2)+(((y-sphD[[3,2]]))^2)>(((sphD[[3,3]]))^2)
,{x,-1,1},{y,-1,1},AspectRatio->Automatic,FrameTicks->None,PlotPoints->40];
Bnd=Graphics[Circle[{0,0},1]];
Dom=Show[FDom,Bnd];

toFundDom[v_]:=CtoR2[HtoD[R2toC[toFundDomH[v]]]];
toFundDomH[v_]:=
Module[{vFD=v,checkMore=True,j},
While[checkMore,
checkMore=False;
For[j=1,j≤numSph,j++,
If[(vFD[[1]]-sph[[j,1]])^2+(vFD[[2]]-sph[[j,2]])^2<sph[[j,3]]^2,
vFD={sph[[j,1]],sph[[j,2]]}+sph[[j,3]]^2(vFD-{sph[[j,1]],sph[[j,2]]})/Norm[(vFD-{sph[[j,1]],sph[[j,2]]})]^2;
vFD[[1]]=2sph[[j,1]]-vFD[[1]];
checkMore=True;
];
];
];
vFD
];

(* Running the code *)

height=Table[x^2,{x,3,.1,-.1}];
pic=height;

(* Make frames; takes a minute *)

For[j=1,j≤Length[height],j++,
list=ListPlot[Table[toFundDom[{x,height[[j]]}],{x,actmin,actmax,xstep}]];
pic[[j]]=Show[Dom,list(*,Graphics[{Darker[Blue],Thickness[.01],Line[{{xmin,height[[j]]},{actmin,height[[j]]}}]}]
,Graphics[{Darker[Blue],Thickness[.01],Line[{{actmax,height[[j]]},{xmax,height[[j]]}}]}]*)];
];

Beep[]
(* A sample frame*)

pic[[25]]

(* This will take a while *)
Export["mov.swf",pic];
Beep[]