(*
    (c) 2010 Alex Kontorovich
*)

(* Parameters *)
sizemax=50; (* max size of spheres *)
MaxNum=30; (* max number of spheres *)
levelmax=30; (* max wordlength of spheres *)
b0={-11,21,28,27,25}; (* Root quintuple *)
labels=True; (* Whether to label circle centers by their curvature *)
fscale=1000; (* toggles size of text *)


(* Code *)
scale=1;
dim=3;
Smat={
{{-1,1,1,1,1},{0,1,0,0,0},{0,0,1,0,0},{0,0,0,1,0},{0,0,0,0,1}},
{{1,0,0,0,0},{1,-1,1,1,1},{0,0,1,0,0},{0,0,0,1,0},{0,0,0,0,1}},
{{1,0,0,0,0},{0,1,0,0,0},{1,1,-1,1,1},{0,0,0,1,0},{0,0,0,0,1}},
{{1,0,0,0,0},{0,1,0,0,0},{0,0,1,0,0},{1,1,1,-1,1},{0,0,0,0,1}},
{{1,0,0,0,0},{0,1,0,0,0},{0,0,1,0,0},{0,0,0,1,0},{1,1,1,1,-1}}
};


des[b1_,b2_,b3_,b4_]:=Solve[ (b1 +b2  +b3 +b4+b5 )^2== 3(b1^2 +b2^2 +b3^2 +b4^2 +b5^2), b5];

touch2[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}];

touch3[a1_,b1_,r1_,a2_,b2_,r2_,a3_,b3_,r3_]:=Solve[{
(x-a1)^2+(y-b1)^2+z^2==(r1)^2,
(x-a2)^2+(y-b2)^2+z^2==(r2)^2,
(x-a3)^2+(y-b3)^2+z^2==(r3)^2
},{x,y,z}];

circPlot[a_,b_,c_,r_]:=Graphics3D[{Opacity[1/r/sizemax],Sphere[{a,b,c},r]}];



morecircs[level_]:=
If[level<levelmax,
Module[{j},
For[j=1,j≤dim+2,j++,

b=Smat[[j]].b;
w=Smat[[j]].w;


If[
(* Not too many spheres*)
(num≤MaxNum)
&&
(* bj not too big *)
(b[[j]]<sizemax)
&&
(* no AA *)
(
from[[level-1]]≠j
)
&&
(* not ABA unless A<B *)
(
from[[level-2]]≠j
||
(
from[[level-2]]==j 
&&
from[[level-1]]>j
)
)
,

num++;

AppendTo[from,j];

AppendTo[quadrups,b];
cNew=circPlot[scale ((w[[j,1]])/(b[[j]])),scale ((w[[j,2]])/(b[[j]])),scale ((w[[j,3]])/(b[[j]])),(scale/(Abs[b[[j]]]))];
p=Graphics3D[{Opacity[Abs[b[[j]]]/sizemax],Text[Style[b[[j]],Abs[fscale/b[[j]]]],{scale ((w[[j,1]])/(b[[j]])),scale ((w[[j,2]])/(b[[j]])),scale ((w[[j,3]])/(b[[j]]))}]}];
AppendTo[pic,cNew];
If[labels,AppendTo[pic,p]];
morecircs[level+1];

from=Drop[from,-1];


];

b=Smat[[j]].b;
w=Smat[[j]].w;
]
];
];

(* Initializations *)
w0={{0,0,0},
{b0[[2]]/Abs[b0[[1]]]- 1,0,0},
blank,blank,blank
};

w0[[3]]=
b0[[3]]{x,y,0}/.touch2[0,0,1/b0[[1]]+1/b0[[3]],w0[[2,1]]/b0[[2]],0,1/b0[[2]]+1/b0[[3]]][[1]];

w0[[4]]=b0[[4]]({x,y,z}/.
(touch3[
0,0,1/b0[[1]]+1/b0[[4]],
w0[[2,1]]/b0[[2]],0,1/b0[[2]]+1/b0[[4]],
w0[[3,1]]/b0[[3]],w0[[3,2]]/b0[[3]],1/b0[[3]]+1/b0[[4]]
])[[2]])/.z->0;

w0[[5]]=b0[[5]]({x,y,z}/.
touch3[
0,0,1/b0[[1]]+1/b0[[5]],
w0[[3,1]]/b0[[3]],w0[[3,2]]/b0[[3]],1/b0[[3]]+1/b0[[5]]
,
w0[[2,1]]/b0[[2]],0,1/b0[[2]]+1/b0[[5]]
][[1]]/.z->0);

cs=ps={blank,blank,blank,blank,blank};
For[j=1,j≤dim+2,j++,
cs[[j]]=circPlot[scale/b0[[j]]w0[[j,1]],scale/b0[[j]]w0[[j,2]],scale/b0[[j]]w0[[j,3]],scale/Abs[b0[[j]]]];
ps[[j]]=Graphics3D[{Opacity[Abs[b0[[j]]]/sizemax],Text[Style[Abs[b0[[j]]],Floor[fscale/Abs[b0[[j]]]]],{scale/b0[[j]]w0[[j,1]],scale/b0[[j]]w0[[j,2]],scale/b0[[j]]w0[[j,3]]}]}];
];
cs[[1]]=Graphics3D[{Opacity[.1+Abs[b0[[1]]]/sizemax],Sphere[w0[[1]]/Abs[b0[[1]]],1/b0[[1]]]}];

pic0={cs[[1]],cs[[2]],cs[[3]],cs[[4]],cs[[5]]};
If[labels,AppendTo[pic0,{ps[[2]],ps[[3]],ps[[4]],ps[[5]]}]];
pic=pic0;

b=b0;
w=N[w0];

num=1;

quadrups={};

from={0,0};

(* Run the code *)

morecircs[3];

pic1=Show[pic,Graphics3D[Text[Style[-11,30],{-.17,0,0}]],Boxed->False]

(* This is to make a movie of the sphere spinning around; it takes a WHILE to rasterize *)
tab=Table[
Rasterize[
Show[pic,ViewPoint->RotationTransform[θ,{0,0,1}][{3,0,3}] , SphericalRegion->True]
,
ImageResolution->70
],{θ,0,2Pi,.05}];
Beep[] (* Beep when done *)

(* Check that the pictures are coming out ok *)
tab[[1]]
(* Check the length; more than 150 slides, and you'll have a 10Mb file... *)
Length[tab]

(* This might take as much as a half hour if not more; it beeps when done *)
Export["movSoddy.swf",tab,"CompressionMethod"->"GZ"];
Beep[]