Apollonian Gasket

Author

Eric W. Weisstein
January 5, 2011

This notebook downloaded from http://mathworld.wolfram.com/notebooks/Fractals/ApollonianGasket.nb.

For more information, see Eric's MathWorld entry http://mathworld.wolfram.com/ApollonianGasket.html.

©2010 Wolfram Research, Inc. except for portions noted otherwise

Additional Authors

Nested version of code by
Michael Trott
The Mathematica Guidebook, Vol. 1: Programming in Mathematica.
New York: Springer-Verlag, 2004.

Code

Code from Mathematica Guidebooks
By Michael Trott

makeTouchingCircles[{p1_,p2_,p3_},iter_,howLong_:10^-3]:=Module[{newton,newCircleData,rMax},(*the derivative for the Newton method*)newton[{{{x1_,y1_},r1_},{{x2_,y2_},r2_},{{x3_,y3_},r3_}},{{xn_,yn_},rn_}]:={{#[[1]],#[[2]]},#[[3]]}&[{xn,yn,rn}-1/2*Inverse[{{(xn-x1),(yn-y1),-(rn+r1)},{(xn-x2),(yn-y2),-(rn+r2)},{(xn-x3),(yn-y3),-(rn+r3)}}].{(xn-x1)^2+(yn-y1)^2-(r1+rn)^2,(xn-x2)^2+(yn-y2)^2-(r2+rn)^2,(xn-x3)^2+(yn-y3)^2-(r3+rn)^2}];
(*the next smaller circle*)newCircleData[{pp1_,pp2_,pp3_}]:={Circle@@#,{{pp1,pp2,#[[1]]},{pp1,pp3,#[[1]]},{pp2,pp3,#[[1]]}}}&[Module[{r12,r23,r13,r1,r2,r3,startx,starty,startr,ε=10^-10},{r12,r23,r13}=Sqrt[#.#]&/@N[{pp1-pp2,pp2-pp3,pp1-pp3}];
r1=(r12+r13-r23)/2;r2=(r12-r13+r23)/2;
r3=(-r12+r13+r23)/2;
startr=Sqrt[#.#]&[N[pp1-({startx,starty}=N[(pp1+pp2+pp3)/3])]]-r1;
(*iterating the Newton method*)FixedPoint[newton[{{pp1,r1},{pp2,r2},{pp3,r3}},#]&,{{startx,starty},startr},SameTest->(#.#&[Flatten[#1-#2]]<ε&)]]];
Join[Module[{r1,r2,r3},(*the start circles*){r12,r23,r13}=Sqrt[#.#]&/@N[{p1-p2,p2-p3,p1-p3}];
r1=(r12+r13-r23)/2;
r2=(r12-r13+r23)/2;
r3=(-r12+r13+r23)/2;
rMax=Max[r1,r2,r3];
{Circle[p1,r1],Circle[p2,r2],Circle[p3,r3]}],(*iterating the calculation of new circles*)Map[First,df=NestList[newCircleData/@Flatten[Map[Last,Select[#,(#[[1,2]]/rMax>howLong)&],{1}],1]&,{newCircleData[{p1,p2,p3}]},iter],{2}]]]

CircleIteration[n_,colors_]:=
Graphics[{circles,
Table[{{colors[[i,k]],disks[[i,k]]/.Circle:>Disk},disks[[i,k]]},{i,n},{k,Length[disks[[i]]]}]
},PlotRange->{{-.4,.4},{-.25,.5}},
AspectRatio->Automatic
]

ApollonianGasket[n_]:=Module[{i},
Drop[
makeTouchingCircles[Table[{Cos[#],Sin[#]}&[2.Pi (i+.25)/3],{i,0,2}],n],
0]
]

Eric's Simple Code
Needs to be rewritten

ApollonianGasket::usage="ApollonianGasket[{v1,v2,v3},n] gives the nth iteration of the Apollonian gasket defined bu circles with centers vi.  There is is currently a limitation of n<=3."

TangentCirclesTriangle[tri_List]:=Module[{i},Table[Circle[TangentCircles[tri][[i,1]],TangentCircles[tri][[i,2]],Sort[{TriSideAngle[tri,{i,Add1[i]}],TriSideAngle[tri,{i,Add1[Add1[i]]}]}]],{i,3}]]

gasket[t_Triangle,1]:=SoddyCircles[t][[2]]

gasket[t_Triangle,2]:=Module[{i,tri=Vertices@t},Table[SoddyCircles[Triangle@{tri[[i]],tri[[Add1[i]]],gasket[tri,1][[1]]}][[2]],{i,3}]]

gasket[t_Triangle,3]:=Module[{i,tri=Vertices@t},{Table[SoddyCircles[Triangle@{tri[[i]],tri[[Add1[i]]],gasket[tri,2][[i,1]]}][[2]],{i,3}],Table[{SoddyCircles[Triangle@{tri[[i]],gasket[tri,2][[i,1]],gasket[tri,1][[1]]}][[2]],SoddyCircles[Triangle@{tri[[Add1[i]]],gasket[tri,2][[i,1]],gasket[tri,1][[1]]}][[2]]},{i,3}]}]

ApollonianGasket2[t_Triangle,n_Integer/;0≤n≤3]:=Module[{i},Graphics[{TangentCirclesTriangle[t],Table[gasket[t,i],{i,n}]}]]

Graphics

Schematic

tri={{0,0},{1,0},{.3,.8}};
Show[GraphicsArray[Partition[
Table[Show[ApollonianGasket[tri,i],AspectRatio->Automatic,DisplayFunction->Identity],
{i,0,3}],4]]]

g=Graphics[ApollonianGasket[2]]

Randomly Colored

V5

its=8;

Print[Timing[g=ApollonianGasket[its];]];

circles=Take[g,3];
disks=Drop[g,3];
colors=Table[Hue[Random[]],{Length[#]}]&/@disks;

Show[GraphicsArray[
Partition[Table[CircleIteration[n,colors],{n,0,5}],3]]]

V6

its=8;

Print[Timing[g=ApollonianGasket[its];]];

circles=Take[g,3];
disks=Drop[g,3];
colors=Table[Hue[Random[]],{Length[#]}]&/@disks;

GraphicsGrid[Partition[Table[CircleIteration[n,colors],{n,0,5}],3]]

Colored by Radius

its=5;

Print[Timing[g=ApollonianGasket[its];]];

circles=Take[g,3];
disks=Drop[g,3];
colors=Table[Hue[1-4disks[[i,k,2]]],{i,its},{k,Length[disks[[i]]]}];

Show[GraphicsArray[
Partition[Table[CircleIteration[n,colors],{n,0,5}],3]]]

Animation

movie=Show/@Table[CircleIteration[n,colors],{n,0,its}];

Export[gifs<>"/apollgas.gif",movie]//Timing

Apollonius Spheres

Code by
Michael Trott
October 5, 2007

(* position a sphere between four given spheres *)
centerSphere[{Sphere[{x1_, y1_, z1_}, r1_], Sphere[{x2_, y2_, z2_}, r2_], 
              Sphere[{x3_, y3_, z3_}, r3_], Sphere[{x4_, y4_, z4_}, r4_]}] := 
Module[{x, y, z, r}, Sphere[Most[#], Last[#]]&[
    Sort[Select[({x, y, z, r} /. #&) /@ 
        Solve[{(x1 - x)^2 + (y1 - y)^2 + (z1 - z)^2 == (r1 + r)^2, 
               (x2 - x)^2 + (y2 - y)^2 + (z2 - z)^2 == (r2 + r)^2, 
               (x3 - x)^2 + (y3 - y)^2 + (z3 - z)^2 == (r3 + r)^2, 
               (x4 - x)^2 + (y4 - y)^2 + (z4 - z)^2 == (r4 + r)^2}, {x, y, z, r}], 
               #[[4]] > 0&], Last[#2] > Last[#1]&][[1]]]]

(* form 4-tuples of spheres*)
newSphereSet[{s1_, s2_, s3_, s4_}] := 
With[{S = centerSphere[{s1, s2, s3, s4}]}, 
     {{S, s1, s2, s3}, {S, s1, s2, s4}, {S, s1, s3, s4}, {S, s2, s3, s4}}]

(* start with spheres at the vertices of a tetrahedron *)
startSpheres = Sphere[#, 1/2]& /@ PolyhedronData["Tetrahedron", "VertexCoordinates"];

gr = Graphics3D[{Directive[Yellow, Specularity[Blue, 15]],
DeleteCases[Union[Flatten[(* iterate sphere adding *)
 Nest[Flatten[newSphereSet[#]& /@ #, 1] &, {startSpheres // N}, 3]], 
                  SameTest -> (Norm[#1[[1]] - #2[[1]]] < 10^-6&)], 
           (* remove the four large original spheres*) Sphere[_, _?(# > 1/3&)]]},
           PlotRange -> All, Boxed -> False]

gr = Graphics3D[{Directive[Yellow, Specularity[Blue, 15]],
DeleteCases[Union[Flatten[(* iterate sphere adding *)
 Nest[Flatten[newSphereSet[#]& /@ #, 1] &, {startSpheres // N}, 6]], 
                  SameTest -> (Norm[#1[[1]] - #2[[1]]] < 10^-6&)], 
           (* remove the four large original spheres*) Sphere[_, _?(# > 1/3&)]]},
           PlotRange -> All, Boxed -> False]

GraphicsRow[
{Show[gr , ViewPoint ->3 Normalize[{0,0,-0.20412414523193154}]] ,
Show[gr , ViewPoint -> 3Normalize[{-0.2886751345948129,-0.5,-0.20412414523193154}]]/.
Yellow->Red,
Show[gr , ViewPoint ->3 Normalize[{-0.2886751345948129,0,-0.20412414523193154}]]/.
Yellow->Purple
}
]

Construction

?CirclePoints

Graphics[{Circle/@CirclePoints[2Sqrt[3]/3,3],Blue,Circle[{0,0},2.15]}]

https://arxiv.org/pdf/math/0101066v1.pdf

?LinearSolve

Q2=SparseArray[{i_,i_}:>1,{4,4},-1]

MatrixForm[%]