Fractals/Apollonian fractals

=Theory= =How to compute and draw Apollonian gasket ? =
 * wikipedia
 * mathforum : Problem of Apollonius

Notation

 * 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


Apollonian gasket can be made using these algorithms:
 * Mandelbrot algorithm ( using circle inversion)
 * Soddy's algorithm ( using circle Descartes theorem )
 * IFS algorithm by Alexei Kravchenko and Dmitriy Mekhontsev

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


(to do )

Soddy's algorithm


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


 * $$( k_{ia}, k_{ib}, k_{ic} ) = (32, 32, 33)$$

This algorithm uses:
 * previous circles defined by centers $$c_n$$ and curvatures $$k_n$$
 * circle Descartes theorem to compute curvature of next circle $$k_4$$
 * $$ k_4 = k_1 + k_2 + k_3 \pm2 \sqrt{k_1 k_2 + k_2 k_3 + k_3 k_1}$$


 * complex circle Descartes theorem to compute center of next circle $$c_4$$
 * $$c_4 = \frac{c_1 k_1 + c_2 k_2 + c_3 k_3 \pm 2 \sqrt{k_1 k_2 c_1 c_2 + k_2 k_3 c_2 c_3 + k_1 k_3 c_1c_3} }{k_4}$$

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 for computing curvature and center

 * $$ f_p(a_1,a_2,a_3) \overset{\underset{\mathrm{def}}{}}{=} a_1 + a_2 + a_3 + 2\cdot\sqrt{a_1 a_2 + a_2 a_3 + a_3 a_1}$$


 * $$ f_{\color{Red}m}(a_1,a_2,a_3) \overset{\underset{\mathrm{def}}{}}{=} a_1 + a_2 + a_3 {\color{Red}-} 2\cdot\sqrt{a_1 a_2 + a_2 a_3 + a_3 a_1}$$

Note that variable $$a_i$$ can be curvature $$k_i$$ or $$c_i*k_i$$ so $$f_p$$ can be used for computing both curvatures and centers

Functions for computing ck list
It is easier to use a list $$ck$$ which defines circle. It consist of two element center $$c$$ and curvature $$k$$ :

$$ck = [c,k]$$

$$c = ck[1]$$

$$k = ck[2]$$

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_p m (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_ m p_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 $$f_{pp}$$ uses $$f_p$$ function for computing both center and curvature. It is used for computing almost all circles without outer circle and hard circles

Function $$f_{pm}$$ uses $$f_p$$ function for computing curvature and $$f_m$$ function for computing center. It is used for computing hard circles.

Function $$f_{mm}$$ uses $$f_m$$ function for computing both center and curvature. It is used for computing outer circle.

Functions for filling gaps
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

 * initial stage : 3 mutually tangent circles = $$( k_{ia}, k_{ib}, k_{ic} ) $$
 * zero stage : 2 circles tangent to previous : outer circle $$ck_{0out}$$ and inner circle $$ck_{0in}$$, where :
 * $$ {\color{Red}ck}_{0out} = f_(ck_{ia}, ck_{ib}, ck_{ic}) $$


 * $$ ck_{0in} = f_{pp}(ck_{ia}, ck_{ib}, ck_{ic})$$


 * 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 :


 * $$ ck_{1a} = f_{pp}(ck_{ia},ck_{ib},ck_{0out})$$
 * $$ ck_{1b} = f_{pp}(ck_{ib},ck_{ic},ck_{0out})$$
 * $$ {\color{Red}c}k_{1c} = f_{p{\color{Red}m}}(ck_{ic},ck_{ia},ck_{0out}) $$
 * $$ ck_{1d} = f_{pp}(ck_{ia},ck_{ib},ck_{0in})$$
 * $$ ck_{1e} = f_{pp}(ck_{ib},ck_{ic},ck_{0in})$$
 * $$ ck_{1f} = f_{pp}(ck_{ic},ck_{ia},ck_{0in})$$

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 $$ f_{pp} $$ function except one hard circle :


 * $$ {\color{Red}c}k_{2cc}:f_{p {\color{Red}m}}(ck_{ic},ck_{1c},ck_{0out}) $$


 * next stages : Put one circle in each gap at each stage. At each stage there is only one hard circle and it is in $$1c$$ gap.

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


 * choose size of quadratic image $$side$$
 * compute radius of outer circle :
 * $$r_{0out} = \frac{side}{2}$$


 * place outer circle in first quadrant of Cartesian plane :
 * $$c_{0out} = r_{0out} + r_{0out}*i$$


 * compute radius of 3 equal initial circles
 * $$a = 1 + \frac{2}{\sqrt{3}}$$


 * $$r_i = \frac{r_{0out}}{a}$$

Then centers of these circles are vertices of equilateral triangle, with side of the triangle $$= 2 * r_i$$ so :
 * find centers of initial circles. For example place one circle $$ia$$ in upper row and 2 circles ($$ib$$ and $$ic$$) in lower row.

$$height = r_i * \sqrt{3}$$

$$c_{ia} = r_{0out} + (2*r_0out - r_i)*i$$

$$c_{ib} = (r_{0out} + r_i) + (2*r_0out - r_i - height)*i$$

$$c_{ic} = (r_{0out} - r_i) + (2*r_0out - r_i - height)*i$$


 * compute inner circle $$0in$$
 * $$ ck_{0in} = f_{pp}(ck_{ia}, ck_{ib}, ck_{ic})$$


 * "Now we have a configuration to start" There is 6 gaps to fill. 5 gaps are easy and one gap $$1c$$ is hard to fill.

Relation between circles

 * $$ {\color{Red}ck}_{0out} = f_(ck_{ia}, ck_{ib}, ck_{ic}) $$


 * $$ ck_{0in} = f_{pp}(ck_{ia}, ck_{ib}, ck_{ic})$$


 * $${\color{Red}ck}_{ia}=f_(ck_{ib},ck_{ib},ck_{0out})$$


 * $$c{\color{Red}k}_{ib}=f_{{\color{Red}m}p}(ck_{ia},ck_{ic},ck_{0out})$$


 * $${\color{Red}ck}_{ic}=f_(ck_{ia},ck_{ib},ck_{0out})$$

Number of circles

 * number of new circles at stage n =  $$2*3^n$$
 * total number of circles after n stages = $$2 + 3^{n+1} $$

where $$n>=0$$

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.
 * do not draw circles which are small ( have a small radius < pixel size)
 * 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
Note that only one curvature has been computed with $$ f_{\color{Red}m}$$ function.

It is $$k_{0out}$$. It is also only one negative curvature.

All the rest have been comuted with $$f_p$$ function.

Haskell and EDSL Diagrams

 * Image and code by Brent Yorgey
 * diagrams code

Java
CirclePack program by Ken Stephenson

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

Common Lisp
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")

JavaScript


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).

=See also=
 * mathcurve : Apollonian gasket and apollonian packing
 * http://packomania.com/
 * My God, It’s Full of Dots! Posted on 27 September 2019 by Brian Hayes

=References=
 * sage code
 * Fractals and Steiner chains by Anders Sandberg with matlab code
 * Looking through the Apollonian Window by Jerzy Kocik