Consideriamo 3 cerchi tangenti a due a due e di centri rispettivamente A(x₁,y₁), B(x₂, y₂), C(x₃,y₃) e di raggi r₁, r₂, r₃.
Vogliamo calcolare centro D(x₄,y₄) e raggio r₄ del cerchio interno tangente a tutti e tre i cerchi.
Calcoliamo r₄ mediante la seguente formula, dovuta a Frederick Soddy (biografia su Wikipedia):

Nella formula scegliamo il segno + davanti alla radice, infatti il segno - si riferisce al cerchio esterno tangente ai tre cerchi.

Sappiamo inoltre che, nel caso di due cerchi tangenti esternamente vale la proprietą che la distanza fra i centri č uguale alla somma dei raggi.

Possiamo dunque scrivere (abbiamo omesso le graffe):

1. (x₁ - x₄)² + (y₁ - y₄)² = (r₁ + r₄)²
2. (x₂ - x₄)² + (y₂ - y₄)² = (r₂ + r₄)²
3. (x₃ - x₄)² + (y₃ - y₄)² = (r₃ + r₄)²

Sottraendo membro a membro la B. dalla A. e  la C. dalla B. e svolgendo i calcoli otteniamo (abbiamo omesso le graffe):

-2x₄(x₁-x₂)-2y₄(y₁-y₂) = r₁²-r₂²+2r₄ (r₁-r₂)-x₁²+x₂²-y₁²+y₂²

-2x₄(x₂-x₃)-2y₄(y₂-y₃) = r₂²-r₃²+2r₄ (r₂-r₃)-x₂²+x₃²-y₂²+y₃²

Si tratta di risolvere un sistema lineare di 2 equazioni in 2 incognite (x₄ e y₄), compito che lasciamo volentieri al computer.

CODICE PER IL FRATTALE

Le tre circonferenze, C_A, C_B, C_C, mutuamente tangenti, tutte di raggio √3/2, hanno centro rispettivamente in
A(1,0)
B(-1/2, +√3/2)
C(-1/2, -√3/2).
Tali punti sono i vertici di un triangolo equilatero inscritto in un cerchio di raggio 1.
Disegniamo le tre circonferenze.
Prepariamo una procedura ricorsiva, di nome apollo, alla quale dovremo passare ordinatamente, come parametri:

1. il numero di passi n
2. le coordinate dei centri ed i raggi delle tre circonferenze xa, ya, ra, xb, yb, rb, xc, yc, rc

Se n = 0 allora
  disegna il cerchio di raggio ra e centro xa e ya
   altrimenti
     raggio = (ra * rb * rc) / (ra * rb + rb * rc + ra * rc + 2 * Sqr(ra * rb * rc * (ra + rb + rc)))
     a1S = -2 * (xb - xc)
     b1S = -2 * (yb - yc)
     c1S = rb ^ 2 - rc ^ 2 + 2 * raggio * (rb - rc) + xc ^ 2 - xb ^ 2 + yc ^ 2 - yb ^ 2
     a2S = -2 * (xa - xb)
     b2S = -2 * (ya - yb)
     c2S = ra ^ 2 - rb ^ 2 + 2 * raggio * (ra - rb) + xb ^ 2 - xa ^ 2 + yb ^ 2 - ya ^ 2
    Se (a1S * b2S - a2S * b1S <> 0) allora
       x0 = (b2S * c1S - b1S * c2S) / (a1S * b2S - a2S * b1S)
       y0 = (a1S * c2S - a2S * c1S) / (a1S * b2S - a2S * b1S)
   Fine se
   xaT = x0
   yaT = y0
   raT = raggio
   xbT = x0
   ybT = y0
   rbT = raggio
   xcT = x0
   ycT = y0
   rcT = raggio
   Chiama apollo(n - 1, xa, ya, ra, xb, yb, rb, xcT, ycT, rcT)
   Chiama apollo(n - 1, xa, ya, ra, xbT, ybT, rbT, xc, yc, rc)
   Chiama apollo(n - 1, xaT, yaT, raT, xb, yb, rb, xc, yc, rc)
 Fine se

Torna al frattale di Apollonio