#
# Load the following
#
precision:=60: # sets the precision, the higher the better (but also slower)
basis:=8: # sets the maximum degree of an algebraic extension in which to try and find numeric values
#
StandardPack:=proc(a) # find radius and centers in a standard packing
global r,x,y,rAlg,xAlg,yAlg;
local n,d,i,j,currentSum,notDone,aa,bb,cc,tt,plotted,nplot; # variables used to keep track
n:=nops(a); # load number of vertices
for i from 1 to n do d[i]:=nops(a[i]); end do; # load degrees
#
# we implement algorithm of Collins and Stephenson to find radii in SOME configuration (namely first three touching and with radii=100) to high precision
#
for i from 1 to 3 do r[i]:=100.0 end do;
for i from 4 to n do r[i]:=8.0 end do;
Digits:=precision+10;
notDone:=1;
while (notDone>0) do
notDone:=0;
for i from 4 to n do
currentSum:=angleSum(i,d,a,r);
if (abs(evalf(currentSum-2*Pi))>1/10^(precision+5)) then notDone:=1; end if;
aa:=sin(currentSum/(2*d[i]));
bb:=sin(Pi/d[i]);
r[i]:=evalf(aa*(1-bb)*r[i]/((1-aa)*bb));
end do;
end do;
#
# with the radii in place we now pack the circles in the plane (done to find the centers of the circles, important for the next step which is inversion)
#
x[1]:=100.0; x[2]:=0.0; x[3]:=200.0;
y[1]:=100.0*sqrt(3.0); y[2]:=0.0; y[3]:=0.0;
plotted[1]:=1; plotted[2]:=1; plotted[3]:=1; # we have placed center of first three circles
for i from 4 to n do plotted[i]:=0; end do;
for i from 1 to n do # we now place center of remaining circles (in order)
for j from 2 to d[i] do
if plotted[a[i][j]]=0 then # we have a circle that hasn't been placed tangent to two which have, now we place it
plotted[a[i][j]]:=1;
aa:=r[i];
bb:=r[a[i][j-1]];
cc:=r[a[i][j]];
tt:=arccos(((aa+bb)^2+(aa+cc)^2-(bb+cc)^2)/(2*(aa+bb)*(aa+cc)));
x[a[i][j]]:=((aa+cc)/(aa+bb))*(cos(tt)*(x[a[i][j-1]]-x[i])+sin(tt)*(y[a[i][j-1]]-y[i]))+x[i];
y[a[i][j]]:=((aa+cc)/(aa+bb))*(-sin(tt)*(x[a[i][j-1]]-x[i])+cos(tt)*(y[a[i][j-1]]-y[i]))+y[i];
end if;
end do;
end do;
#
# with the radii and the centers in place, we now invert at point of tangency between circles 1 and 2 finding the new radii then scale so that largest radius is 1
#
r[1]:=0; r[2]:=0;
for i from 3 to n do
r[i]:=r[i]/((50.0-x[i])^2+(50.0*sqrt(3.0)-y[i])^2-(r[i])^2);
if r[i]>aa then aa:=r[i]; end if;
end do;
for i from 4 to n do
r[i]:=r[i]/r[3];
end do;
r[3]:=1;
#
# we now have the radii of all the circles in the standard packing, so we now pack them (a little subtler than previously since we now have two zero bend circles to contend with; we have not tried to optimize the code here!)
#
for i from 4 to n do plotted[i]:=0; end do; plotted[1]:=-1; plotted[2]:=-1; plotted[3]:=1;
x[0]:=0; y[0]:=0; x[1]:=0; y[1]:=0;
x[3]:=0; y[3]:=0; # anchor circle three at (0,0)
aa:=0.0;
for i from 3 to d[1] do # position circles along the bottom
plotted[i+1]:=1;
y[i+1]:=-1.0+r[i+1];
x[i+1]:=aa+2.0*sqrt(r[i]*r[i+1]);
aa:=x[i+1];
end do;
nplot:=n-1-d[1]; #used to keep track of how many circles have already been placed.
while (nplot>0) do
for i from 1 to n do
if (plotted[i]=1) then # circle is already placed so see if we can place any of its neighbors
for j from 2 to d[i] do
if ((plotted[a[i][j-1]]=1) and (plotted[a[i][j]]=0)) then
plotted[a[i][j]]:=1; nplot:=nplot-1;
aa:=r[i];
bb:=r[a[i][j-1]];
cc:=r[a[i][j]];
tt:=arccos(((aa+bb)^2+(aa+cc)^2-(bb+cc)^2)/(2*(aa+bb)*(aa+cc)));
x[a[i][j]]:=((aa+cc)/(aa+bb))*(cos(tt)*(x[a[i][j-1]]-x[i])+sin(tt)*(y[a[i][j-1]]-y[i]))+x[i];
y[a[i][j]]:=((aa+cc)/(aa+bb))*(-sin(tt)*(x[a[i][j-1]]-x[i])+cos(tt)*(y[a[i][j-1]]-y[i]))+y[i];
end if;
end do;
if ((plotted[a[i][d[i]]]=1) and (plotted[a[i][1]]=0)) then
plotted[a[i][1]]:=1; nplot:=nplot-1;
aa:=r[i];
bb:=r[a[i][d[i]]];
cc:=r[a[i][1]];
tt:=arccos(((aa+bb)^2+(aa+cc)^2-(bb+cc)^2)/(2*(aa+bb)*(aa+cc)));
x[a[i][1]]:=((aa+cc)/(aa+bb))*(cos(tt)*(x[a[i][j-1]]-x[i])+sin(tt)*(y[a[i][j-1]]-y[i]))+x[i];
y[a[i][1]]:=((aa+cc)/(aa+bb))*(-sin(tt)*(x[a[i][j-1]]-x[i])+cos(tt)*(y[a[i][j-1]]-y[i]))+y[i];
end if;
end if;
end do;
end do;
#
# have MAPLE attempt to identify numbers algebraically (but first round)
#
Digits:=precision;
for i from 1 to n do
r[i]:=evalf(10^(-precision)*trunc(10^(precision)*r[i]));
x[i]:=evalf(10^(-precision)*trunc(10^(precision)*x[i]));
y[i]:=evalf(10^(-precision)*trunc(10^(precision)*y[i]));
rAlg[i]:=identify(evalf(r[i]),all=false,Rational=true,Algebraic=true,BasisSizeAlg=basis);
xAlg[i]:=identify(evalf(x[i]),all=false,Rational=true,Algebraic=true,BasisSizeAlg=basis);
yAlg[i]:=identify(evalf(y[i]),all=false,Rational=true,Algebraic=true,BasisSizeAlg=basis);
end do;
#
# test to see that the algebraic answers satisfy the tangency conditions
#
Digits:=precision+5;
aa:=1;
for i from 2 to d[1] do #test along bottom
if (simplify(yAlg[a[1][i]]+1-rAlg[a[1][i]])<>0) then aa:=0; end if;
end do;
for i from 2 to d[2] do #test along top
if (simplify(yAlg[a[2][i]]-1+rAlg[a[2][i]])<>0) then aa:=0; end if;
end do;
for i from 3 to n do #test in the interior
for j from 1 to d[i] do
if (i<a[i][j]) then # test adjacency relationship
if (simplify((xAlg[i]-xAlg[a[i][j]])^2+(yAlg[i]-yAlg[a[i][j]])^2-(rAlg[i]+rAlg[a[i][j]])^2)<>0) then aa:=0; end if;
end if;
end do;
end do;
#
# give output
#
if (aa=1) then
printf("MAPLE found the algebraic answers, in addition to the horizontal lines y=1 and y=-1 the rest of the circles are listed (algebraically) below as [x,y,r]");
for i from 3 to n do print([xAlg[i],yAlg[i],rAlg[i]]); end do;
else
printf("MAPLE did not find the algebraic answers, in addition to the horizontal lines y=1 and y=-1 the rest of the circles are listed (numerically) below as [x,y,r]");
for i from 3 to n do print([x[i],y[i],r[i]]); end do;
end if;
end proc:
#
# A function to find angle sum around a circle
#
angleSum:=proc(v,d,a,r)
local i,rSum;
rSum:=0;
for i from 1 to d[v]-1 do
rSum:=rSum+arccos(((r[v]+r[a[v][i]])^2+(r[v]+r[a[v][i+1]])^2-(r[a[v][i]]+r[a[v][i+1]])^2)/(2*(r[v]+r[a[v][i]])*(r[v]+r[a[v][i+1]])));
end do;
rSum:=rSum+arccos(((r[v]+r[a[v][1]])^2+(r[v]+r[a[v][d[v]]])^2-(r[a[v][1]]+r[a[v][d[v]]])^2)/(2*(r[v]+r[a[v][1]])*(r[v]+r[a[v][d[v]]])));
evalf(rSum);
end proc:#
# StandardPack.mw - a user's guide:
#
# This is a supplement to "Irreducible Apollonian configurations and packings"
# by Steve Butler, Ron Graham, Gerhard Guettler and Colin Mallows
# (code written by Steve Butler. Last updated 4 May 2009)
#
# This MAPLE worksheet is designed to take a tangency graph of an Apollonian
# configuration and produce the centers and radius of a standard packing. The
# approach is straightforward and can be found by examining the above code for
# the procedure StandardPack.
#
# StandardPack is the main function and it works as follows:
# > StandardPack(A);
# where A is an adjacency list of the tangency graph. As an output MAPLE will
# first find the radius and centers to high precision (which can be modified
# above) and then will attempt to identify the numbers algebraically. If it
# can succesfully identify all of the numbers algebraically and these satsify
# the tangency conditions it will output the algebraic numbers, otherwise
# it will output the numbers to some number of digits (as set above).
#
# In addition after calling the function it will store all of the information
# in the variables x,y,r,xAlg,yAlg,rAlg where x,y,r are the high precision
# numeric approximations for the centers and radii while xAlg,yAlg,rAlg are
# MAPLE's attempt to identify them algebraically. The users can then play with
# these numbers. Note that sometimes MAPLE will be able to identify most but
# not all of the numbers algebraically and will thus ouput the high precision
# numbers. So it is good to look at the outputs xAlg,yAlg,rAlg. This can be
# done for example using the following command.
# > print(rAlg);
# We might then be able to fill in any of the missing pieces.
#
# The convention for creating an adjacency list from a triangulated graph is
# to first pick an orientation in the plane (i.e., clockwise or counter-
# clockwise), then pick a vertex and label it 1, pick a vertex adjacent to 1
# and label it 2. The rest is mechanical. Starting with vertex 1 we make a
# list of all the vertices adjacent to it by starting with the lowest label
# (in this case 2) and going in the direction of the orientation we list the
# vertices we meet, where any unlabeled vertex we enounter receives the
# lowest label not yet used. We then go to the next vertex, find the lowest
# label it is adjacent to and repeat the procedure. We keep this up until we
# have made a list for every vertex. So for example K_4 embedded in the plane
# would have A=[[2,3,4],[1,4,3],[1,2,4],[1,3,2]], i.e., A is the set of lists
# of adjacent vertices. Order is IMPORTANT in these lists!
#
# (This convention follows that given in the plantri documentation of
# Gunnar Brinkmann and Brendan McKay.)
#
# The standard packing that is produced corresponds to where the vertex
# labelled 1 is the line y=-1, the vertex labelled 2 is the line y=1 and the
# vertex labelled 3 is the unit circle at the origin. So if there is a
# specific standard packing that you are after choose the vertices accorgingly.
#
# Examples of a few packings from the paper are given below. They are a good
# place to start looking at how examples work.
##
# The tetrahedron graph - gives Apollonian packing
#
StandardPack([[2,3,4],[1,4,3],[1,2,4],[1,3,2]]);#
# The octahedron graph - gives Guettler-Mallow packing
#
StandardPack([[2,3,4,5],[1,5,6,3],[1,2,6,4],[1,3,6,5],[1,4,6,2],[2,5,4,3]]);#
# The icosahedron graph - an uninvestigated packing
#
StandardPack([[2,3,4,5,6],[1,6,7,8,3],[1,2,8,9,4],[1,3,9,10,5],[1,4,10,11,6],[1,5,11,7,2],[2,6,11,12,8],[2,7,12,9,3],[3,8,12,10,4],[4,9,12,11,5],[5,10,12,7,6],[7,11,10,9,8]]);#
# The subdivision rule shown in Figure 10
#
StandardPack([[2,3,4,5,6,7],[1,7,8,9,10,3],[1,2,10,4],[1,3,10,11,12,5],[1,4,12,13,14,6],[1,5,14,7],[1,6,14,15,8,2],[2,7,15,9],[2,8,15,16,17,10],[2,9,17,11,4,3],[4,10,17,16,13,12],[4,11,13,5],[5,12,11,16,18,14],[5,13,18,15,7,6],[7,14,18,16,9,8],[9,15,18,13,11,17],[9,16,11,10],[13,16,15,14]]);#
# The subdivision rule shown in Figure 7
#
StandardPack([[2,3,4,5,6,7],[1,7,8,9,10,3],[1,2,10,11,12,4],[1,3,12,13,14,5],[1,4,14,6],[1,5,14,7],[1,6,14,15,8,2],[2,7,15,9],[2,8,15,10],[2,9,15,13,11,3],[3,10,13,12],[3,11,13,4],[4,12,11,10,15,14],[4,13,15,7,6,5],[7,14,13,10,9,8]]);#
# The packing shown in Figure 13
#
StandardPack([[2,3,4,5,6,7],[1,7,8,9,10,11,12,13,14,3],[1,2,14,15,16,4],[1,3,16,5],[1,4,16,15,17,6],[1,5,17,7],[1,6,17,15,8,2],[2,7,15,9],[2,8,15,10],[2,9,15,11],[2,10,15,12],[2,11,15,13],[2,12,15,14],[2,13,15,3],[3,14,13,12,11,10,9,8,7,17,5,16],[3,15,5,4],[5,15,7,6]]);#
# The packing shown in Figure 17
#
StandardPack([[2,3,4,5,6,7,8,9],[1,9,10,11,12,13,14,3],[1,2,14,4],[1,3,14,5],[1,4,14,6],[1,5,14,15,16,7],[1,6,16,8],[1,7,16,17,18,19,12,11,10,9],[1,8,10,2],[2,9,8,11],[2,10,8,12],[2,11,8,19,20,13],[2,12,20,14],[2,13,20,21,22,15,6,5,4,3],[6,14,22,16],[6,15,22,21,20,17,8,7],[8,16,20,18],[8,17,20,19],[8,18,20,12],[12,19,18,17,16,21,14,13],[14,20,16,22],[14,21,16,15]]);