Cercle tangent à trois cercles 2 à 2 tangents

library(plotrix) # for the draw.circle function

# radius of the inner Soddy circle ####
soddyRadius <- function(r1,r2,r3){
  1/(1/r1+1/r2+1/r3+2*sqrt(1/r1/r2+1/r2/r3+1/r3/r1))
}

# center of the inner Soddy circle (equal detour point) ####
soddyCenter <- function(A,B,C){
  measure <- function(M,N){
    sqrt(c(crossprod(N-M)))
  }
  a <- measure(B,C)
  b <- measure(A,C)
  c <- measure(A,B)
  AB <- B-A
  AC <- C-A
  Delta <- 0.5*abs(AB[1]*AC[2]-AB[2]*AC[1]) # area of triangle ABC
  tc <- 1 + 2*Delta/c(a*(b+c-a), b*(c+a-b), c*(a+b-c)) # triangular coordinates
  den <- c(tcrossprod(t(tc), c(a,b,c)))
  k1 <- a*tc[1]/den
  k2 <- b*tc[2]/den
  k3 <- c*tc[3]/den
  k1*A + k2*B + k3*C
}

# case of an exterior circle C3 (using Descartes' theorems) ####
soddyCircles <- function(C1, C2, C3){
  r1 <- C1$radius; r2 <- C2$radius; r3 <- -C3$radius
  soddyR <- soddyRadius(r1,r2,r3)
  z1 <- C1$center[1] + 1i*C1$center[2]
  z2 <- C2$center[1] + 1i*C2$center[2]
  z3 <- C3$center[1] + 1i*C3$center[2]
  term1 <- 1/r1*z1+1/r2*z2+1/r3*z3
  term2 <- 2*sqrt(1/r1/r2*z1*z2+1/r2/r3*z2*z3+1/r1/r3*z1*z3)
  soddyCenter1 <- (term1 - term2)*soddyR
  soddyCenter2 <- (term1 + term2)*soddyR
  soddyCenter1 <- c(Re(soddyCenter1), Im(soddyCenter1))
  soddyCenter2 <- c(Re(soddyCenter2), Im(soddyCenter2))
  sols <- list(
    list(center=soddyCenter1, radius=soddyR),
    list(center=soddyCenter2, radius=soddyR)
  )
  d1 <- c(crossprod(soddyCenter1-C3$center))
  d2 <- c(crossprod(soddyCenter2-C3$center))
  if(d1>d2) sols[[1]] else sols[[2]]
}

# draw a circle ####
drawCircle <- function(C, ...){
  draw.circle(C$center[1], C$center[2], C$radius, ...)
}

# starting circles ####
C1 <- list(center=c(1,0),radius=1)
C2 <- list(center=c(-1,0),radius=1)
C3 <- list(center=c(0,sqrt(3)),radius=1)
C4 <- list(center = soddyCenter(C1$center,C2$center,C3$center), 
           radius = 2 + soddyRadius(C1$radius,C2$radius,C3$radius))
  
# main function (recursive) ####
Apollony <- function(C1, C2, C3, n, colors=FALSE, exterior=FALSE, ...){
  if(exterior){
    Cnew <- soddyCircles(C1,C2,C3)
  }else{
    soddyCenter <- soddyCenter(C1$center,C2$center,C3$center)
    radius <- soddyRadius(C1$radius,C2$radius,C3$radius)
    Cnew <- list(center=soddyCenter, radius=radius)
  }
  drawCircle(Cnew, col=if(colors) Colors[n] else NA, ...)
  if(n==1) return(invisible())
  Apollony(Cnew,C2,C3,n-1, colors=colors, exterior = exterior)
  Apollony(C1,Cnew,C3,n-1, colors=colors, exterior = exterior)
  Apollony(C1,C2,Cnew,n-1, colors=colors)
}

# Apollony fractal ####
par(mar=c(0,0,0,0))
plot(0, 0, type="n", axes=FALSE, xlab=NA, ylab=NA, ylim=c(-1.6,2.74), asp=1)
drawCircle(C1); drawCircle(C2); drawCircle(C3); drawCircle(C4)
n <- 5
Apollony(C1, C2, C3, n)
Apollony(C1, C2, C4, n, exterior=TRUE)
Apollony(C1, C3, C4, n, exterior=TRUE)
Apollony(C2, C3, C4, n, exterior=TRUE)