(*
	(c) 2008 Alex Kontorovich
*)

(* First make the Underlying Apollonian circle packing *)
b0={-10,18,23,27}; (* root quadruple *)
fscale=2000; (* text size *)

scale=1;
CirThru[a1_,b1_,a2_,b2_,a3_,b3_]:=Solve[{
(a1-x0)^2+(b1-y0)^2==(r)^2,(a2-x0)^2+(b2-y0)^2==(r)^2,(a3-x0)^2+(b3-y0)^2==(r)^2},{x0,y0,r}];des[b1_,b2_,b3_]:=Solve[ (b1 +b2  +b3 +b4 )^2== 2(b1^2 +b2^2 +b3^2 +b4^2 ), b4];touch[a1_,b1_,r1_,a2_,b2_,r2_]:=Solve[{(x-a1)^2+(y-b1)^2==(r1)^2,(x-a2)^2+(y-b2)^2==(r2)^2},{x,y}];
circPlot[a_,b_,r_]:=Graphics[ Circle[{a,b},r]];

b0[[1]]=b4/.des[b0[[2]],b0[[3]],b0[[4]]][[1]];
w0={blank,0,b0[[3]]I (1/b0[[2]]+1/b0[[3]]),blank};
w0[[4]]=b0[[4]](x+I y) /.touch[0,0,1/b0[[2]]+1/b0[[4]],0,1/b0[[2]]+1/b0[[3]],1/b0[[3]]+1/b0[[4]]][[2]];
w0[[1]]=b4/.des[w0[[2]],w0[[3]],w0[[4]]][[1]];
cs=ps={blank,blank,blank,blank};
For[j=1,j≤4,j++,
cs[[j]]=circPlot[scale/b0[[j]]Re[w0[[j]]],scale/b0[[j]]Im[w0[[j]]],scale/Abs[b0[[j]]]];
ps[[j]]=Graphics[Text[Style[b0[[j]],fscale/b0[[j]]],{scale Re[((w0[[j]])/(b0[[j]]))],scale Im[((w0[[j]])/(b0[[j]]))]}]];
];
b=b0;
w=w0;
pic0={cs[[1]],cs[[2]],cs[[3]],cs[[4]],ps[[2]],ps[[3]],ps[[4]]};
pic=pic0;
Show[pic]

(* Now make the dual quadruple *)
touchs={
{blank,blank,blank,blank},
{blank,blank,blank,blank},
{blank,blank,blank,blank},
{blank,blank,blank,blank}
};
For[i=1,i<4,i++,
For[j=i+1,j≤4,j++,
touchs[[i,j]]=x+I y/.touch[
Re[w[[j]]/b[[j]]],
Im[w[[j]]/b[[j]]],
1/b[[j]],
Re[w[[i]]/b[[i]]],
Im[w[[i]]/b[[i]]],
1/b[[i]]
][[1]];
];
];
dualCs={blank,blank,blank,blank};
s={blank,blank,blank,blank};
sph={blank,blank,blank,blank};
For[j=1,j≤4,j++,
oth=Select[{1,2,3,4},#!=j&];
{rad,xval,yval}=({r,x0,y0}/.
CirThru[
Re[touchs[[oth[[1]],oth[[2]]]]],
Im[touchs[[oth[[1]],oth[[2]]]]],
Re[touchs[[oth[[1]],oth[[3]]]]],
Im[touchs[[oth[[1]],oth[[3]]]]],
Re[touchs[[oth[[2]],oth[[3]]]]],
Im[touchs[[oth[[2]],oth[[3]]]]]
][[2]]);
dualCs[[j]]=Graphics[{Red,Thick, Dashed,Circle[{scale xval,scale yval},scale rad]}];
s[[j]]=Graphics3D[Sphere[{xval,yval,0},rad]];
sph[[j]]={xval,yval,rad};
];

(* Show the dual on top of the root *)
pic1=Show[
cs[[1]],cs[[2]],cs[[3]],cs[[4]],
dualCs[[1]],dualCs[[2]],dualCs[[3]],dualCs[[4]]
]

(* Now show a fundamental domain for this action in H^3 *)
bot=Plot3D[0,{x,-20,20},{y,-20,20}]; (* bottom is the boundary of H^3 *)
fund=Show[s[[1]],s[[2]],s[[3]],s[[4]],bot,PlotRange->{{-1.2,.4},{-.5,.6},{0,.5}},Boxed->True,Axes->True]

(* Now for low-lying horocycles*)

(* send point to exterior of four spheres *)
toFundDom[v_]:=
Module[{vFD=v,checkMore=True,j},
While[checkMore,
checkMore=False;
For[j=1,j≤4,j++,
If[(vFD[[1]]-sph[[j,1]])^2+(vFD[[2]]-sph[[j,2]])^2+(vFD[[3]])^2<sph[[j,3]]^2,
vFD={sph[[j,1]],sph[[j,2]],0}+sph[[j,3]]^2(vFD-{sph[[j,1]],sph[[j,2]],0})/Norm[(vFD-{sph[[j,1]],sph[[j,2]],0})]^2;
checkMore=True;
];
];
];
vFD
]

(* Here is an individual screen shot *)
height=.115;
list0=ListPointPlot3D[Flatten[Table[toFundDom[{x,y,height}],{x,-1.5,.5,.01},{y,-1,1,.01}],1],PlotStyle->Opacity[.4]];
pic0=Show[fund,list0,PlotRange->{{-1.5,0.5},{-1,1},{0,2}} (*,ViewPoint->{2,10,1.125}*)]

(* Now for a movie *)

(* make the frames; this will take a while *)
tab=Table[
Show[fund,
ListPointPlot3D[Flatten[Table[toFundDom[{x,y,height^2}],{x,-1.5,.5,.01},{y,-1,1,.01}],1],PlotStyle->Opacity[.4]]
,PlotRange->{{-1.5,0.5},{-1,1},{0,2}},ViewPoint->{2,2,.125}
]
,{height,.8,.01,-.01}
];
Beep[]

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