settings.render=2;
settings.outformat="pdf";

// -----------------------------------------------------------------------------
import three;
import solids;
import animate;

// -----------------------------------------------------------------------------
real[] circumsphere(triple p1, triple p2, triple p3, triple p4){
  real a = orient(p1,p2,p3,p4);
  real q1 = p1.x*p1.x + p1.y*p1.y + p1.z*p1.z;
  real q2 = p2.x*p2.x + p2.y*p2.y + p2.z*p2.z;
  real q3 = p3.x*p3.x + p3.y*p3.y + p3.z*p3.z;
  real q4 = p4.x*p4.x + p4.y*p4.y + p4.z*p4.z;
  real Dx = orient((q1,p1.y,p1.z),(q2,p2.y,p2.z),(q3,p3.y,p3.z),(q4,p4.y,p4.z));
  real Dy = orient((q1,p1.x,p1.z),(q2,p2.x,p2.z),(q3,p3.x,p3.z),(q4,p4.x,p4.z));
  real Dz = orient((q1,p1.x,p1.y),(q2,p2.x,p2.y),(q3,p3.x,p3.y),(q4,p4.x,p4.y));
  triple center = 0.5/a * (Dx,-Dy,Dz);
  real r = length(p1-center);
  real[] out = {center.x,center.y,center.z,r};
  return out;
}

// -----------------------------------------------------------------------------
triple inversion(triple pnt, real r, triple center){
  triple omega = (2*r, 0, 0);
  real k = 3*r*r;
  triple v = pnt - omega - center;
  return omega + center - k/(v.x*v.x+v.y*v.y+v.z*v.z)*v;
}

// -----------------------------------------------------------------------------
real[] oneSphere(triple center, real r, real beta){
  triple pnt = center + (r*cos(beta)*2/3, r*sin(beta)*2/3, 0);
  real thRadius = r/3;
  triple p1 = pnt + (thRadius,0,0);
  triple p2 = pnt + (0,thRadius,0);
  triple p3 = pnt + (-thRadius,0,0);
  triple p4 = pnt + (0,0,thRadius);
  triple q1 = inversion(p1, r, center);
  triple q2 = inversion(p2, r, center);
  triple q3 = inversion(p3, r, center);
  triple q4 = inversion(p4, r, center);
  real[] cs = circumsphere(q1,q2,q3,q4);
  cs[0] = cs[0] - 4*r;
  return cs;
}

// -----------------------------------------------------------------------------
real[][] hexlet(triple center, real r, int k, bool clockwise){
  real xi = k*pi/90;
  if(clockwise){
    xi = -xi;
  }
  real[][] out = new real[6][4];
  for(int i=0; i<6; ++i){
    real[] s = oneSphere(center, r, i*pi/3+xi);
    for(int j=0; j<4; ++j){
      out[i][j] = s[j];
    }
  }
  return out;
}

real[][] hexlet2(int n, int k, bool clockwise){
  if(n==1){
    return hexlet((0,0,0), 1, k, clockwise);
  }
  real[][] hx = hexlet2(n-1, k, !clockwise);
  real[][] out = new real [0][4];
  for(int i=0; i<6^(n-1); ++i){
    real s[] = hx[i];
    real[][] hxnew = hexlet((s[0],s[1],s[2]), s[3], k, !clockwise);
    out.append(hxnew);
  }
  return out;
}

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

animation anim;

int n = 2;
for(int k=0; k<30; ++k){
  erase();
  real[][] hx = hexlet2(n, k, true);
  for(int i=0; i<6^n; ++i){
    real[] s = hx[i];
    draw(surface(sphere((s[0],s[1],s[2]), s[3])), red);
  }
  anim.add();
}

anim.movie();