settings.render = 4;
settings.outformat = "eps";

// import modules --------------------------------------------------------------
import graph;
import three;
import solids;

// overall settings ------------------------------------------------------------
currentprojection = orthographic(0,0,6);
viewportmargin = (10,10);
size(10cm);
currentlight = ((3,3,0));

// -----------------------------------------------------------------------------
struct circle {
  pair center;
  real radius;
}

// radius of inner Soddy circle ------------------------------------------------
real soddyRadius(real r1, real r2, real r3){
  return 1/(1/r1+1/r2+1/r3+2*sqrt(1/r1/r2+1/r2/r3+1/r3/r1));
}

// center of inner Soddy circle ------------------------------------------------
pair soddyCenter(pair A, pair B, pair C){
  pair AB = B-A;
  pair AC = C-A;
  real a = length(C-B); real b = length(AC); real c = length(AB);
  // area of triangle ABC:
  real Delta = 0.5*abs(AB.x*AC.y-AB.y*AC.x);
  // triangular coordinates:
  real tc1 = 1 + 2*Delta/(a*(b+c-a));
  real tc2 = 1 + 2*Delta/(b*(c+a-b));
  real tc3 = 1 + 2*Delta/(c*(a+b-c));
  //
  real den = a*tc1 + b*tc2 + c*tc3;
  real k1 = a*tc1/den; real k2 = b*tc2/den; real k3 = c*tc3/den;
  return k1*A + k2*B + k3*C;
}

// inner Soddy circle ----------------------------------------------------------
circle soddyCircle1(circle C1, circle C2, circle C3){
  circle out;
  out.radius = soddyRadius(C1.radius, C2.radius, C3.radius);
  out.center = soddyCenter(C1.center, C2.center, C3.center);
  return out;
}

// case of an exterior center C3 -----------------------------------------------
circle soddyCircle2(circle C1, circle C2, circle C3){
  real r1 = C1.radius; real r2 = C2.radius; real r3 = -C3.radius;
  real r = soddyRadius(r1,r2,r3);
  pair z1 = C1.center; pair z2 = C2.center; pair z3 = C3.center;
  pair term1 = 1/r1*z1 + 1/r2*z2 + 1/r3*z3;
  pair term2 = 2*sqrt(1/r1/r2*z1*z2 + 1/r2/r3*z2*z3 + 1/r1/r3*z1*z3);
  pair center1 = r*(term1-term2);
  pair center2 = r*(term1+term2);
  real d1 = (center1.x-z3.x)^2 + (center1.y-z3.y)^2;
  real d2 = (center2.x-z3.x)^2 + (center2.y-z3.y)^2;
  circle out;
  out.radius = r;
  if(d1>d2){
    out.center = center1;
  }else{
    out.center = center2;
  }
  return out;
}

// starting circles ------------------------------------------------------------
circle C1; C1.center = (1,-1/sqrt(3)); C1.radius = 1;
circle C2; C2.center = (-1,-1/sqrt(3)); C2.radius = 1;
circle C3; C3.center = (0,sqrt(3)-1/sqrt(3)); C3.radius = 1;
circle C4; C4.center = (0,0); C4.radius = 2+soddyRadius(1,1,1);

// main function ---------------------------------------------------------------
void Apollony(circle C1, circle C2, circle C3, int n, bool exterior){
  circle Cnew;
  if(exterior){
    Cnew = soddyCircle2(C1,C2,C3);
  }else{
    Cnew = soddyCircle1(C1,C2,C3);
  }
  draw(surface(sphere((Cnew.center.x,Cnew.center.y,0), Cnew.radius)), red);
  if(n>1){
    Apollony(Cnew, C2, C3, n-1, exterior);
    Apollony(C1, Cnew, C3, n-1, exterior);
    Apollony(C1, C2, Cnew, n-1, false);
  }
}

// plot ------------------------------------------------------------------------
int n = 6;
Apollony(C1, C2, C3, n, false);
Apollony(C1, C2, C4, n, true);
Apollony(C1, C3, C4, n, true);
Apollony(C2, C3, C4, n, true);