Internet regorge d'informations sur l'empilement de cercles dans un cercle. Je ne vais pas donner plus d'explications ici. Vous pouvez commencer par lire la page WIKIPEDIA sur ce sujet.

Entrons tout de suite dans le vif du sujet. Premier point, la caméra avec ces paramètres  :

camera {
  orthographic
  location <0, 0, -48>
  up y
  right x*image_width/image_height
  look_at <0, 0, 0>
  angle 36
  }

La caméra utilise des rayons parallèles pour créer l'image, ce qui signifie qu'il n'y a pas de distorsion des objets due à la perspective. La distance par rapport à la caméra est fixée à 48 et l'angle de vue horizontal à 36°. Avec ces paramètres, le plus grand cercle visible avec son centre en <0,0,0> a un rayon de 15 :Dans l'image ci-dessous, c'est le cercle de centre C4.



Prenons un rayon au hasard pour le cercle C1, disons 5,25 par exemple. Il est tangent intérieurement au cercle C4. Le cercle C2 a alors un rayon de 9,75. Il est aussi tangent intérieurement au C4 mais à l'opposé du diamètre. Il ne reste plus qu'à trouver le rayon du cercle C3 et la position de son centre. C'est le problèmes d'Apollonius (enfin, une partie). Lorsque l'on parle de cercles, d'intersections et de tangentes, on tombe très souvent sur la résolution d'une équation quadratique.

Nous utiliserons cette macro pour en trouver les solutions.

#macro ResolveQuadratic (q1,q2,q3,Off)
  #local Ct = 0.5;
  #local b = -(q1+q2+q3);
  #local c = q1*q1+q2*q2+q3*q3-0.5*(q1+q2+q3)*(q1+q2+q3)-Off;
  #local d = b*b-4*Ct*c;
  #local d = (abs(d)<0.000001?0.0:d);
  array {(-b+sqrt(d))/(2.0*Ct), (-b-sqrt(d))/(2.0*Ct)}
#end

Notez que la macro renvoie deux valeurs dans un tableau. Notez aussi que le test pour le déterminant n'est pas fait par rapport au zéro strict, mais avec une tolérance de 1e-6.

Nous devons introduire maintenant la notion de "courbure" qui, pour un cercle vaut simplement \( k = 1/r \). C'est une valeur signée. Il est possible de calculer le rayon et le centre du cercle C3 de cette façon. :

#declare R = 15;

#declare b4 = -1/R;
#declare c4 = <0,0,0>;
#declare r4 = abs(1/b4);

#declare c1 = <1.3*R/2, 0, 0>;
#declare r1 = R - c1.x;
#declare b1 = 1/r1;

#declare c2 = ;
#declare r2 = c1.x;
#declare b2 = 1/r2;

#declare b3 = ResolveQuadratic(b4,b1,b2,0);

#declare bx3 = ResolveQuadratic(b4*c4.x,b1*c1.x,b2*c2.x,2);
#declare c3X = bx3[0]/b3[0];

#declare by3 = ResolveQuadratic(b4*c4.y,b1*c1.y,b2*c2.y,2);
#declare c3Y = by3[0]/b3[0];

#declare c3 = <c3X, c3Y, 0>;
#declare r3 = c3R;
#declare b3 = 1/r3;

Pour ne pas garder un aspect top symétrique, l'ensemble des trois cercles a subi une rotation autour de l'axe z. Cette rotation est appliquée grace à la macro :

#macro RotZ(Pt,Angle)
  <Pt.x*cosd(Angle)-Pt.y*sind(Angle), Pt.x*sind(Angle)+Pt.y*cosd(Angle), 0>
#end

Les fonctions sind ( ) et cosd ( ) font parties de la librairie math.inc. J'ai choisi ici un angle de 38°. En demandant d'afficher les valeurs, on obtient :

#declare C1 = < 9.750000,  0.000000,  5.250000>;
#declare C2 = <-5.250000,  0.000000,  9.750000>;
#declare C3 = < 5.825243,  8.834951,  4.417476>;
#declare C4 = < 0.000000,  0.000000, 15.000000>;

Comme la composante z n'est pas utilisée, je m'en sers pour y stocker le rayon. On retrouvera ces valeurs de base dans le code source complet des images.



Ci dessous deux extraits des images en 2000x2000 pixels réalisées avec POVRay.

Le niveau de récursion est de 12 pour l'image en noir et blanc et 14 pour celle en couleurs. Comme cette figure est une fractale, il n'y a pas de limite théorique. Mais il faut en fixer une avec nos petits ordinateurs et c'est le rôle de ce paramètre qui permet de limiter le calcul et le dessin des cercles. Plus ce niveau est élevé, plus le rayon des cercles diminuent, plus le temps de calcul augmente et plus l'utilisation de la mémoire augmente aussi. Avec 16Go de RAM je n'ai pas pu dépasser 14.

Pour l'image en noir et blanc, les cercles sonr représentés par des tores { }. Pour l'image en couleurs, se sont des cylinder { }. On peut très bien imaginer aussi utiliser des sphere { }











Pour finir, deux liens qui ont largement inspiré ce travail :

• The Coding Train : #182 — Apollonian Gasket
• JSXGraph : Apollonian circle packing