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])
  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
}

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

# main function (recursive) ####
Apollony <- function(C1, C2, C3, n, colors=FALSE, ...){
  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, ...)
  #print("x")
  if(n==1) return(invisible())
  Apollony(Cnew,C2,C3,n-1, colors=colors)
  Apollony(C1,Cnew,C3,n-1, colors=colors)
  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(0,sqrt(3)/2), asp=1)
drawCircle(C1); drawCircle(C2); drawCircle(C3)
n <- 7
Colors <- rainbow(n)
Apollony(C1, C2, C3, n, colors=TRUE)