def kfun(k1,k2,k3,k4):

"""

The Descartes formula for the curvature of an inverted tangent circle.

"""

return 2*k1+2*k2+2*k3-k4

def circfun(c1,c2,c3,c4):

"""

Computes the inversion of circle 4 in the first three circles.

"""

newk = kfun(c1[3],c2[3],c3[3],c4[3])

newx = (2*c1[0]*c1[3]+2*c2[0]*c2[3]+2*c3[0]*c3[3]-c4[0]*c4[3])/newk

newy = (2*c1[1]*c1[3]+2*c2[1]*c2[3]+2*c3[1]*c3[3]-c4[1]*c4[3])/newk

if newk > 0:

newr = 1/newk

elif newk < 0:

newr = -1/newk

else:

newr = Infinity

return [newx, newy, newr, newk]

def mcircle(circdata, label = False, thick = 1/10):

"""

Draws a circle from the data. label = True

"""

if label==True and circdata[3] > 0 and circdata[2] > 1/2000:

lab = text(str(circdata[3]),(circdata[0],circdata[1]), fontsize = \

500*(circdata[2])^(.95), vertical_alignment = 'center', horizontal_alignment \

= 'center', rgbcolor = (0,0,0))

else:

lab = Graphics()

circ = circle((circdata[0],circdata[1]), circdata[2], rgbcolor = (0,0,0), \

thickness = thick)

return circ+lab

def add_circs(c1, c2, c3, c4, cutoff = 300):

"""

Find the inversion of c4 through c1,c2,c3. Add the result to circlist,

then (if the result is big enough) recurse.

"""

newcirc = circfun(c1, c2, c3, c4)

if newcirc[3] < cutoff:

circlist.append(newcirc)

add_circs(newcirc, c1, c2, c3)

add_circs(newcirc, c2, c3, c1)

add_circs(newcirc, c3, c1, c2)

zst1 = [0,0,1/2,-2]

zst2 = [1/6,0,1/3,3]

zst3 = [-1/3,0,1/6,6]

zst4 = [-3/14,2/7,1/7,7]

circlist = [zst1,zst2,zst3,zst4]

add_circs(zst1,zst2,zst3,zst4)

add_circs(zst2,zst3,zst4,zst1)

add_circs(zst3,zst4,zst1,zst2)

add_circs(zst4,zst1,zst2,zst3)

sum([mcircle(q, label = True, thick = 1/2) for q in \

circlist]).save('./Apollonian3.png',axes = False, figsize = [12,12], xmin = \

-1/2, xmax = 1/2, ymin = -1/2, ymax = 1/2)

## Thursday, January 22, 2009

### Integral Apollonian Packings with Sage

At the national math meetings this year I heard about some really interesting and fun work on integral Apollonian circle packings. The AMS has a nice introductory article about them. I couldn't resist trying to compute and draw some in Sage. Carl Witty greatly improved the speed of my first attempt, so what follows can be considered joint work of ours. (Sage can compute nicer PDF output of these, but blogger doesn't embed PDFs.)

Subscribe to:
Post Comments (Atom)

## 2 comments:

Nice article!

I have a problem. I need to put 15 balls inside one ball so that the balls take the maximum volume.

I guess that the procedure in your article may work.

Thanks, I wouldn't have known where to start!

Post a Comment