Fractals/Apollonian fractals

Theory

edit
  • wikipedia[1]
  • mathforum : Problem of Apollonius [2]

How to compute and draw Apollonian gasket ?

edit

Notation

edit
  • Apollonian gasket = curvilinear Sierpinski gasket = Leibniz packing = Apollonian packing
  • stage = level = step
  • incircle = inscribed circle = circle inside 3 circles = inner circle
  • gap = curvilinear triangle = ideal triangle ( because the sides are tangent at each vertix and the angle between them is zero) = region between 3 tangent circles

Algorithms

edit
 
Ancillary figure for van Roomen's solution to the problem of Apollonius.

Apollonian gasket can be made using these algorithms:[3]

  • Mandelbrot algorithm [4]( using circle inversion)[5][6]
  • Soddy's algorithm ( using circle Descartes theorem )[7][8]
  • IFS algorithm by Alexei Kravchenko and Dmitriy Mekhontsev [9]

All algorithms give state after n stages. It is an approximation of limit set which is made of infinite number of stages (cicles).

Mandelbrot algorithm

edit
 
Gasket made from 4 equal inner circles made with circle inversion
 
Animation of the solution to Apollonius' problem using the method of inversion.

(to do )

Soddy's algorithm

edit
 
Integer Apollonian gasket defined by circle curvatures of (32,32,33). This image is probably made of 19 stages

It will be explained by example gasket with initial configuration : 3 inner circles with integer curvatures

 

This algorithm uses:

  • previous circles defined by centers   and curvatures  
  • circle Descartes theorem to compute curvature of next circle  
 
  • complex circle Descartes theorem to compute center of next circle  
 

When center of outer circle is at origin one has to chose solutions in case of centers. It is easier when the all gasket is in the first quadrant of Cartesian plane. Then all centers have positive parts ( real and imaginary ).

Functions

edit
Functions for computing curvature and center
edit
 
 

Note that variable   can be curvature   or   so   can be used for computing both curvatures and centers

Functions for computing ck list
edit

It is easier to use a list   which defines circle. It consist of two element center   and curvature   :

 

 

 

Now one can define new functions ( Maxima CAS code) for computing 4-th circle mutually tangent to 3 given circles:

f_pp(ck1,ck2,ck3):=block
([c4,k4,ck4],
k4:f_p(ck1[2],ck2[2],ck3[2]),
c4:f_p(ck1[1]*ck1[2],ck2[1]*ck2[2],ck3[1]*ck3[2])/k4,
ck4:[c4,k4])$
f_pm(ck1,ck2,ck3):=block
([c4,k4,ck4],
 k4:f_p(ck1[2],ck2[2],ck3[2]),
 c4:f_m(ck1[1]*ck1[2],ck2[1]*ck2[2],ck3[1]*ck3[2])/k4,
 ck4:[c4,k4]);
f_mm(ck1,ck2,ck3):=block
([c4,k4,ck4],
 k4:f_m(ck1[2],ck2[2],ck3[2]),
 c4:f_m(ck1[1]*ck1[2],ck2[1]*ck2[2],ck3[1]*ck3[2])/k4,
 ck4:[c4,k4]);
 f_mp_ck4(ck1,ck2,ck3):=block
([c4,k4],
k4:solve_m_eq(ck1[2],ck2[2],ck3[2]),
c4:solve_p_eq(ck1[1]*ck1[2],ck2[1]*ck2[2],ck3[1]*ck3[2])/k4,
return([c4,k4])
);

Function   uses   function for computing both center and curvature. It is used for computing almost all circles without outer circle and hard circles

Function   uses   function for computing curvature and   function for computing center. It is used for computing hard circles.

Function   uses   function for computing both center and curvature. It is used for computing outer circle.

Functions for filling gaps
edit
fill_easy(ckla,cklb,cklc,max_stage):=block
([ckm],
 if max_step>0 then (
  /* circle in the middle */
  ckm:give_pp_ck4(ckla,cklb,cklc),
  /* 3 circles around */
  fill_easy(ckla,cklb,ckm,max_stage-1),
  fill_easy(cklb,cklc,ckm,max_stage-1),
  fill_easy(ckla,cklc,ckm,max_stage-1)))$
fill_hard(ck_ic,ck_ia,ck_0out,max_stage):=block
([ck_1c,ck_2cc,ck_3ccc,ck_4ccca,ck_5cccac], /* hard circles */
 if 	max_stage>0 then (
  /* step 1 = circle in the middle */
  ck_1c:f_pm(ck_ia,ck_ic,ck_0out),      /* 1c */
  /* step 2 = 3 circles around */
  fill_easy(ck_ia,ck_0out,ck_1c,max_stage-1), /* 2ca */
  fill_easy(ck_ia,ck_ic,ck_1c,max_stage-1), /* 2cb */
  /* hard subgap 2cc */
  if max_stage>1 then ck_2cc:f_pm(ck_ic,ck_1c,ck_0out), /* 2cc */
  /* step 3 for subgap 2cc */
  fill_easy(ck_0out,ck_1c,ck_2cc,max_stage-2), /* 3cca */
  fill_easy(ck_ic,ck_1c,ck_2cc,max_stage-2),  /* 3ccb */
  if max_stage>2 then ck_3ccc:f_pm(ck_ic,ck_0out,ck_2cc),	/* 3ccc */
  /* step 4 for subgap 3ccc */
  if max_stage>3 then ck_4ccca:f_pm(ck_0out,ck_2cc,ck_3ccc),
  fill_easy(ck_ic,ck_2cc,ck_3ccc,max_stage-3),  /* 4cccb */
  fill_easy(ck_ic,ck_0out,ck_3ccc,max_stage-3),  /* 4cccc */
  /* step 5 for subgap 4ccca */
  fill_easy(ck_0out,ck_2cc,ck_4ccca,max_stage-4),  /*  5cccaa */
  fill_easy(ck_2cc,ck_3ccc,ck_4ccca,max_stage-4),  /* 5cccab */
  if max_stage>4 then ck_5cccac:f_pm(ck_0out,ck_3ccc,ck_4ccca)))$

Construction

edit
 
Apollonian gasket with marked hard circles
  • initial stage : 3 mutually tangent circles =  
  • zero stage : 2 circles tangent to previous : outer circle   and inner circle  , where :
 
 
  • first stage : there are 6 gaps ( = curvilinear triangles) between previous circles. Three gaps are peripheral ( near outer circle) and three are medial ( near inner circle). Put one circle in each gap :
 
 
 
 
 
 

Each new circle placed in the gap makes three new gaps (trifurcation).

  • second stage : there are 18 (6*3) new gaps. Put one circle in each gap using   function except one hard circle :
 
  • next stages : Put one circle in each gap at each stage. At each stage there is only one hard circle and it is in   gap.

Algorithm for programmers

edit

Because all circles are inscribed in outer circle so it is easier to start from it.

  • choose size of quadratic image  
  • compute radius of outer circle :
 
  • place outer circle in first quadrant of Cartesian plane :
 
  • compute radius of 3 equal initial circles [10]
 
 
  • find centers of initial circles. For example place one circle   in upper row and 2 circles (  and  ) in lower row.

Then centers of these circles are vertices of equilateral triangle, with side of the triangle   so :

 

 

 

 

  • compute inner circle  
 
  • "Now we have a configuration to start" [11] There is 6 gaps to fill. 5 gaps are easy and one gap   is hard to fill.

Relation between circles

edit
 
 
 
 
 

Number of circles

edit
  • number of new circles at stage n =  
  • total number of circles after n stages =  

where  

For example :

  • stage -1 = initial configuration ( 3 inner circles )
  • stage 0 gives 2 new circles ( one outer and one most inner ). All circles = 5
  • stage 1 gives 6 new circles ( all circles = 11 )
  • stage 2 gives 18 new circles ( all circles = 29 )
  • stage 3 gives 54 new circles ( all circles = 83)
  • stage 4 gives 162 new circles ( all circles = 245)
  • stage 5 gives 486 new circles ( all circles = 731)
  • stage 6 gives 1 458 new circles ( all circles = 2 189 )
  • stage 7 gives 4 374 new circles ( all circles = 6 563 )
  • stage 8 gives 13 122 new circles ( all circles = 19 685 )
  • stage 9 gives 39 366 new circles ( all circles = 59 051 )
  • stage 10 gives 118 098 new circles ( all circles = 177 149 )
  • stage 11 gives 354 294 new circles ( all circles = 531 443 )
  • stage 12 gives 1 062 882 new circles ( all circles = 1 594 325 )
  • stage 13 gives 3 188 646 new circles ( all circles = 4 782 971 )
  • stage 14 gives 9 565 938 new circles ( all circles = 14 348 909 )
  • stage 15 gives 28 697 814 new circles ( all circles = 43 046 723 )
  • stage 16 gives 86 093 442 new circles ( all circles = 129 140 165 )
  • stage 17 gives 258 280 326 new circles ( all circles = 387 420 491 )
  • stage 18 gives 774 840 978 new circles ( all circles = 1 162 261 469 )
  • stage 19 gives 2 324 522 934 new circles ( all circles = 3 486 784 403 )
  • stage 20 gives 6 973 568 802 new circles ( all circles = 10 460 353 205 )

It can be computed using this Maxima CAS code :

 NumberOfStageCircles(stage):= 
  if stage =-1 then 3
  elseif stage=0 then 2
  else 2*3^stage;

 NumberOfAllCircles(stage_max):=sum(NumberOfStageCircles(i),i,-1,stage_max);

for i:-1 step 1 thru 20 do print("stage ",i," gives ",NumberOfStageCircles(i)," new circles ( all circles =  ",NumberOfAllCircles(i), " )");

Number of circles grows rapidly and makes problems caused by big file size. One can :

  • use only stages up to 5-7. It seems to be sufficient to give good image in a short time.[12]
  • do not draw circles which are small ( have a small radius < pixel size) [13]
  • convert svg file to png ( for example using Image Magic : "convert a.svg a.png" ). Then for stage 12 svg image has 128 MB and png file only 57 KB.
  • be happy that you have so detailed picture and wait a long time to load file or see how system hangs up (:-))

Curvatures

edit

Note that only one curvature has been computed with   function.

It is  . It is also only one negative curvature.

All the rest have been comuted with   function.

stage   curvatures
-1  
0  
1  
2  

Example programs

edit
Haskell and EDSL Diagrams
edit
  • Image and code by Brent Yorgey[14]
  • diagrams code [15]
Java
edit

CirclePack program by Ken Stephenson

Maxima CAS
edit

Program is explained above, and src is in description of Apollonian_rc_6.svg

Common Lisp
edit

To run put apollonian.lisp file in your home directory and in Emacs/slime. To get another step change *max-step*

Load file to run code ( batch mode ) :

(load  "apollonian.lisp")
; first run
;(require :asdf)
;(require :asdf-install)
;(asdf-install:install :zpng)
;(asdf-install:install :Vecto)
; you must press 2 and 0 when the program asks

; next times load packages from disk
(asdf:operate 'asdf:load-op 'vecto)

(in-package vecto)

;-----------------------------------------------------------------------
;---------------------------- definitions of functions ---------------
;-----------------------------------------------------------------------

; kc is a list describing circle (list k c)  where k is curvature and c is a center ( complex number)
; k can be positive ( inner circle) , negative ( outer circle) or zero ( =  line )

(defun draw-vecto-circle (kc)
  (let* ((c (second kc))
	 (r (abs (/ 1 (first kc)))))	; r = abs(1/k)
    (centered-circle-path (realpart c) (imagpart c) r))) ; vecto

; Desfirsttes' circle theorem 
(defun solve-equation-m (k1 k2 k3)
  (- (+ k1 k2 k3 ) (* 2 (sqrt (+ (* k1 k2) (* k2 k3) (* k1 k3))))))

(defun solve-equation-p (k1 k2 k3)
  (+ (+ k1 k2 k3 ) (* 2 (sqrt (+ (* k1 k2) (* k2 k3) (* k1 k3))))))

; compute, draws and gives back 4-th circle usind solve-equation-p and solve-equation-p
(defun draw-forth-circle-pp (kc1 kc2 kc3)
  (let* ((k1 (first kc1))
	 (k2 (first kc2))
	 (k3 (first kc3))
	 (p1 (* k1 (second kc1))) ; pn = kn * cn 
	 (p2 (* k2 (second kc2)))
	 (p3 (* k3 (second kc3)))
	 (k4 (solve-equation-p k1 k2 k3)) ; find k
	 (c4 (/ (solve-equation-p p1 p2 p3) k4)) ; find c
	 (kc4 (list k4 c4)))
    (draw-vecto-circle kc4) ; draw 
    kc4)) ; return 4-th circle

; compute, draws and gives back 4-th circle usind solve-equation-p and solve-equation-m
(defun draw-forth-circle-pm (kc1 kc2 kc3)
  (let* ((k1 (first kc1))
	 (k2 (first kc2)) 
	 (k3 (first kc3))
	 (p1 (* k1 (second kc1)))
	 (p2 (* k2 (second kc2)))
	 (p3 (* k3 (second kc3)))
	 (k4   (solve-equation-p k1 k2 k3))
	 (c4 (/ (solve-equation-m p1 p2 p3) k4))
	 (kc4 (list k4 c4)))
    (draw-vecto-circle kc4)
    kc4))

; easy = pp
(defun fill-easy-gap (kc1 kc2 kc3 steps)
  (when (> steps 0)
      (let* ((kc4 (draw-forth-circle-pp kc1 kc2 kc3))) ; 4-th circle
	(fill-easy-gap kc1 kc2 kc4 (- steps 1)) ; 3 subgaps
	(fill-easy-gap kc2 kc3 kc4 (- steps 1))  
	(fill-easy-gap kc3 kc1 kc4 (- steps 1)))))

; only for gap 1c = hard gap
(defun fill-hard-gap (kc-ic kc-ia kc-0out steps)
  (let* (kc-1c kc-2cc  kc-3ccc  kc-4ccca kc-5cccac) ; hard circles , also kc-6cccacc
    ;------------ gap 1c --------------------------------------------   
    (when (> steps 0) (setq kc-1c (draw-forth-circle-pm kc-ic kc-ia kc-0out)) 						; kc-1c hard
	  (fill-easy-gap kc-ia kc-0out kc-1c (- steps 1))								; kc-2ca
	  (fill-easy-gap kc-ia kc-ic kc-1c (- steps 1))									; kc-2cb
	  ; ----------- hard subgap 2cc ------------------------------------
	  (when (> steps 1)  (setq kc-2cc (draw-forth-circle-pm kc-ic kc-0out kc-1c))					; kc-2cc hard
		(fill-easy-gap kc-0out kc-1c kc-2cc (- steps 2))							; kc-3cca
		(fill-easy-gap kc-ic kc-1c kc-2cc (- steps 2))	        						; kc-3ccb
		; ----------- hard subgap 3ccc ------------------------------------
		(when (> steps 2)   (setq kc-3ccc (draw-forth-circle-pm kc-ic kc-0out kc-2cc))				; kc-3ccc hard
		      (fill-easy-gap kc-ic kc-2cc kc-3ccc (- steps 3))							; kc-4cccb
		      (fill-easy-gap kc-ic kc-0out kc-3ccc (- steps 3))							; kc-4cccc
		      ; ----------- hard subgap 4ccca ------------------------------------
		      (when (> steps 3)  (setq kc-4ccca (draw-forth-circle-pm kc-0out kc-2cc kc-3ccc))			; kc-4ccca hard
			    (fill-easy-gap kc-0out kc-2cc kc-4ccca (- steps 4))						; kc-5cccaa
			    (fill-easy-gap kc-2cc kc-3ccc kc-4ccca (- steps 4))						; kc-5cccab
			    (when (> steps 4) (setq kc-5cccac (draw-forth-circle-pm kc-0out kc-3ccc kc-4ccca)))))))))   ; kc-5cccac hard

; example of use : (draw-apollonian-gasket "a.png" 800 2)
; classical Desfirsttes configuration : 3 mutually tangent circles and 4-th inside ( and fith outside)
(defun draw-apollonian-gasket-n3 (file side step)
  (with-canvas (:width side :height side) ; vecto
	       (set-rgb-stroke 0 0 0) ; vecto
	       (set-line-width 1) ; vecto

	       
	      (let* (	(r0out (/ side 2))
			(k0out (- (/ 1 r0out)))
			(c0out (complex r0out r0out))
			(kc0out (list k0out c0out)) ; list defining circle 
					; a1:float(1 + 2 / sqrt(3)),  http://www2.stetson.edu/~efriedma/cirincir/ 
			(a (+ 1 (/ 2 (sqrt 3))))
			(ri   (/ r0out a)) ; three equal inner circles = initial circles
			(ki (/ 1 ri))	
			(cia (complex r0out (- side ri))) ; one circle in upper row
			(kcia (list ki cia))
			(h (* ri (sqrt 3))) ; height of equilateral triangle with side = 2*ri
					; 2 circles in lower row
			(cib (complex  (+ r0out ri) (- side (+ h ri))))
			(cic (complex  (- r0out ri) (- side (+ h ri))))
			(kcib (list ki cib))
			(kcic (list ki cic))
			kc0in ) 
			; drawing code
			; step -1 = three equal inner circles
		 (draw-vecto-circle kcia ) 		; ia
		 (draw-vecto-circle kcib ) 		; ib
		 (draw-vecto-circle kcic ) 		; ic
			; step 0 = outer and most inner circle
		 (draw-vecto-circle kc0out )		; 0out
		 (setq kc0in ( draw-forth-circle-pp kcia kcib kcic )) ; 0in
			; this is starting configuration. 
			; One can go to the next steps now
			; Fill 6 gaps :   
			; 3 peripheral gaps
		 (fill-easy-gap kcia kcib kc0out step) ; 1a
		 (fill-easy-gap kcib kcic kc0out step) ; 1b	
		 (fill-hard-gap kcic kcia kc0out step) ; 1c
			; 3 medial gaps
		 (fill-easy-gap kcia kcib kc0in step)  ; 1d
		 (fill-easy-gap kcib kcic kc0in step)  ; 1e
		 (fill-easy-gap kcic kcia kc0in step)  ; 1f
			; rest of drawing code
		 (stroke)) ; before save ( vecto procedure)
	       (print (save-png file)))) ; info text and vecto procedure

;---------------compile ---------------------------------
(compile 'draw-vecto-circle )
(compile 'solve-equation-m)
(compile 'solve-equation-p)
(compile 'draw-forth-circle-pp)
(compile 'draw-forth-circle-pm)
(compile 'fill-easy-gap)
(compile 'fill-hard-gap)
(compile 'draw-apollonian-gasket-n3)

;----------global var ----------------------

(defparameter *max-step* 0 " maximal step of drawing apollonian gasket from. It is an integer >= 0 ")

(defparameter *file-name*
  (make-pathname 
   :name (concatenate 'string "apollonian-n3-s" (write-to-string *max-step*))
   :type "png")
  "name (or pathname) of png file ")

;====================================== run ==================================
(draw-apollonian-gasket-n3 *file-name* 800 *max-step*)
JavaScript
edit
 
Random Apollonian circle fractal

This is an incremental algorithm that maintains a queue of triples of touching circles. The main loop starts at the front of the queue, removing a triple and filling it with the required circle, adding the resulting circle to a list of circles to draw as SVG and adding the three triples that result from adding the circle to the back of the computation queue. As a bonus feature, the circles are coloured so that no two neighbours have the same colour.

Usage: concatenate the two SVG boilerplate parts and the middle Javascript part into one SVG file and load into a web browser (or another SVG viewer that supports Javascript or ECMAScript).

<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN" "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
<svg
 xmlns="http://www.w3.org/2000/svg"
 xmlns:xlink="http://www.w3.org/1999/xlink"
 width="100%" height="100%" viewBox="0 0 768 768"
 onload="init(evt)"
><title>Apollonian Fractal</title
><desc>Randomized Apollonian fractal circle packing</desc
><script type="text/ecmascript"><![CDATA[
// SVG properties
var ns = "http://www.w3.org/2000/svg";
var xns = "http://www.w3.org/1999/xlink";
var SVGDocument = null;
var size = 768; // should match 'viewBox'
var minsize = 0.001;

// convert internal coordinates and lengths to SVG
function coord(x) { return x*size/2 + size/2; }
function dist(d) { return d*size/2; }

// complex numbers
function C(r,i)    { this.r = r; this.i = i; }
function Cconj(c)  { return new C(c.r, -c.i); }
function Cadd(c,d) { return new C(c.r+d.r,c.i+d.i); }
function Csub(c,d) { return new C(c.r-d.r,c.i-d.i); }
function Cmul(c,d) { return new C(c.r*d.r-c.i*d.i,c.r*d.i+c.i*d.r); }
function Csqrt(c)  { return new C(Math.sqrt(Math.sqrt(c.r*c.r+c.i*c.i)) * Math.cos(Math.atan2(c.i,c.r)/2),
                                  Math.sqrt(Math.sqrt(c.r*c.r+c.i*c.i)) * Math.sin(Math.atan2(c.i,c.r)/2)); }

// sort points into order
function clockwise(a, b, c) {
  var x1 = a.x-b.x; var y1 = a.y-b.y;
  var x2 = c.x-b.x; var y2 = c.y-b.y;
  if ((x1*y2-y1*x2) >= 0) return [a,b,c]; else return [a,c,b];
}

// compute the 4th circle touching 3 circles, each of which touch the other two
function Kiss(a, b, c, initial) {
  var k1 = 1 / a.r; var z1 = new C(a.x, a.y); var kz1 = Cmul(new C(k1,0),z1);
  var k2 = 1 / b.r; var z2 = new C(b.x, b.y); var kz2 = Cmul(new C(k2,0),z2);
  var k3 = 1 / c.r; var z3 = new C(c.x, c.y); var kz3 = Cmul(new C(k3,0),z3);
  var k4p = k1 + k2 + k3 + 2*Math.sqrt(k1*k2 + k2*k3 + k3*k1);
  var k4m = k1 + k2 + k3 - 2*Math.sqrt(k1*k2 + k2*k3 + k3*k1);
  var kz4p = Cadd(Cadd(Cadd(kz1,kz2),kz3),Cmul(new C(2,0),
    Csqrt(Cadd(Cadd(Cmul(kz1,kz2),Cmul(kz2,kz3)),Cmul(kz3,kz1)))));
  var kz4m = Csub(Cadd(Cadd(kz1,kz2),kz3),Cmul(new C(2,0),
    Csqrt(Cadd(Cadd(Cmul(kz1,kz2),Cmul(kz2,kz3)),Cmul(kz3,kz1)))));
  var k4;
  var kz4;
  var k4b;
  var kz4b;
  var cs = new Array();
  if (k4p > k4m) { k4 = k4p; kz4 = kz4p; k4b = k4m; kz4b = kz4m; }
  else           { k4 = k4m; kz4 = kz4m; k4b = k4p; kz4b = kz4p; }
  var cc = new Circle(kz4.r/k4,kz4.i/k4,Math.abs(1/k4), 6 - a.colour - b.colour - c.colour);
  var dx = a.x - cc.x
  var dy = a.y - cc.y
  var dr = a.r + cc.r
  if (Math.abs(dx * dx + dy *dy - dr * dr) > 0.0001) {
    cc = new Circle(kz4b.r/k4,kz4b.i/k4,Math.abs(1/k4), 6 - a.colour - b.colour - c.colour);
  }
  cs[cs.length] = cc;
  if (initial) {
    cc = new Circle(kz4b.r/k4b,kz4b.i/k4b,Math.abs(1/k4b), 6 - a.colour - b.colour - c.colour);
    cs[cs.length] = cc;
  }
  return cs;
}

// called once on load
function init(evt) {
  // get document
  SVGDocument = evt.target.ownerDocument;

  // initial bounding circle
  var b = new Circle(0, 0, -1, 0);

  // insert two randomly positioned touching circles
  var tr = 1-Math.random()/2;
  var pa = Math.PI/2 - Math.asin(Math.random()*(1-tr)/tr);
  var px = tr * Math.cos(pa);
  var py = tr * Math.sin(pa);
  var pr = 1 - tr;
  var qr = (1 - pr) * (1 - Math.cos(pa + Math.PI/2))
         / (1 + pr - (1 - pr) * Math.cos(pa + Math.PI/2));
  var qx = 0;
  var qy = qr - 1;
  var p = new Circle(px, py, pr, 1);
  var q = new Circle(qx, qy, qr, 2);

  // a queue to contain triples of touching circles
  var queue = new Array();
  var cs = Kiss(b,p,q, true);
  queue[queue.length] = [b,p,cs[0]];
  queue[queue.length] = [b,cs[0],q];
  queue[queue.length] = [cs[0],p,q];
  queue[queue.length] = [b,p,cs[1]];
  queue[queue.length] = [b,cs[1],q];
  queue[queue.length] = [cs[1],p,q];

  // a queue to contain circles to draw
  var draw = new Array();
  draw[draw.length] = b;
  draw[draw.length] = p;
  draw[draw.length] = q;
  draw[draw.length] = cs[0];
  draw[draw.length] = cs[1];

  // add 10000 more circles to the draw queue
  // adding new triples to the compute queue
  var c;
  for (c = 0; c < Math.min(queue.length, 10000); c = c + 1) {
    var c0 = queue[c][0];
    var c1 = queue[c][1];
    var c2 = queue[c][2];
    var ncs = Kiss(c0,c1,c2);
    var nc = ncs[0];
    if (nc.r > minsize) {
      queue[queue.length] = [nc,c1,c2];
      queue[queue.length] = [c0,nc,c2];
      queue[queue.length] = [c0,c1,nc];
      draw[draw.length] = nc;
    }
  }

  // add all generated circles to an SVG <g> element
  var g = SVGDocument.createElementNS(ns, "g");
  g.setAttributeNS(null, "stroke", "black");
  g.setAttributeNS(null, "stroke-width", "1");
  var i;
  for (i = 0; i < draw.length; i = i + 1) {
    g.appendChild(draw[i].draw());
  }

  // link the <g>s into the DOM so they are displayed
  var svg = SVGDocument.documentElement;
  svg.appendChild(g);

} // init()

// circle class
function Circle(x, y, r, colour) {
  // properties of a circle
  this.x = x;
  this.y = y;
  this.r = r;
  this.colour = colour;
  // convert to svg
  this.draw = function() {
    var colours = ["white", "red", "yellow", "cyan"];
    var c = SVGDocument.createElementNS(ns, "circle");
    c.setAttributeNS(null, "fill", colours[this.colour]);
    c.setAttributeNS(null, "cx", coord(this.x));
    c.setAttributeNS(null, "cy", coord(this.y));
    c.setAttributeNS(null, "r",  dist(this.r));
    return c;
  };
} // Circle()
]]></script></svg>

See also

edit

References

edit
  1. Apollonian gasket
  2. mathforum : Problem of Apollonius
  3. Apollonius' Tangency Problem at MathPages
  4. Mandelbrot's Algorithm at Fractal Geometry from Yale University by Michael Frame, Benoit Mandelbrot, and Nial Neger
  5. Drawing 2D Apollonian gasket with n identical circles in Matlab by Guillaume JACQUENOT
  6. Apollonian Gasket in Mathematica and Povray by Paul Nylander
  7. Drawing the Apollonian Gasket with Common Lisp and Vecto by Luis Diego Fallas
  8. Leibnitz packing by Takaya Iwamoto with program in AutoLisp
  9. Apollonian Gasket - by Paul Bourke in BAsic and C
  10. How to pack n circles inside unit circle by Erich Friedman
  11. SVG Math Animation Example: Apollonian Gasket (-15,32,32,33) by KiHyuck Hong. See the src code.
  12. SVG Math Animation Example: Apollonian Gasket (-15,32,32,33) by KiHyuck Hong
  13. File Apolonian gasket containes about 8140 circles, but it looks like made of about 20 steps. So it should have about 10 460 353 205 circles
  14. Apollonian gasketin Haskell and Diagrams by Brent Yorgey
  15. Apollonian.hs - haskell code