/************************
    apollonian.asy
    May 18, 2009
    see http://www.ams.org/featurecolumn/archive/kissing.html
 ***********************/

size(2.7 inch,0);

real Infinity=100000.;

// *****************   objects
struct kcircle {
  pair O; real r; real k; int serial;

  void operator init(explicit pair O, explicit real r, explicit real k) {
    this.O=O;
    this.r=r;
    this.k=k;
  }

  path path() {return circle(O,r);}
  void draw(pen p = currentpen) {
      draw(circle(O,r),p);
  }
  void label(pen p = currentpen) {
    if (k>0)
      label(scale(25*r)*string(k),O);
  }
  void count(int i) {
    serial = i;
  }
}


struct kcirclelist {
  kcircle circlist[]; int count;

  void operator init(explicit kcircle C) {
    this.circlist[0]=C;
    this.count=1;
  }

  void append(kcircle C) {
    circlist[count] = C;
    C.count(count);
    count = count + 1;
  }
  
  void draw(bool showlabel=true) {
    for (int i=0; i< count; i+=1) {
      circlist[i].draw();
      if (showlabel)
	circlist[i].label();
    }
  }  


}

// *****************   working functions

real kfun(real k1,real k2,real k3,real k4) {
  //    The Descartes formula for getting the curvature of
  //    the partner C4 tangent circle.
  return(2*k1+2*k2+2*k3-k4);
}

kcircle circfun(kcircle c1,kcircle c2,kcircle c3,kcircle c4) {
  //    returns the partner kcircle of c4 with respect to the
  //    the first three circles.
  real newk = kfun(c1.k,c2.k,c3.k,c4.k);
  real newx = (2*c1.O.x*c1.k+2*c2.O.x*c2.k+2*c3.O.x*c3.k-c4.O.x*c4.k)/newk;
  real newy = (2*c1.O.y*c1.k+2*c2.O.y*c2.k+2*c3.O.y*c3.k-c4.O.y*c4.k)/newk;
  real newr;
  if (newk > 0) {
    newr = 1/newk;
  }
  else {
    if (newk < 0) {
      newr = -1/newk;
    }
    else newr = Infinity;
  }
  return(kcircle((newx, newy), newr, newk));
}


void add_circs(kcircle c1,kcircle c2,kcircle c3, kcircle c4,
	       kcirclelist circlist, int cutoff = 300) {
  /*
    Find the partner of c4 through c1,c2,c3.  Add the result to circlist,
    then (if the circle is big enough) keep the recursion going.
  */
  kcircle newcirc = circfun(c1, c2, c3, c4);
  if (newcirc.k < cutoff) {
    circlist.append(newcirc);
    add_circs(newcirc, c1, c2, c3, circlist, cutoff);
    add_circs(newcirc, c2, c3, c1, circlist, cutoff);
    add_circs(newcirc, c3, c1, c2, circlist, cutoff);
  }
}



/************************************
Make usage changes down here.
 ************************************/

// *****************   opening Descartes configuration

kcircle zst1=kcircle((0.,0.),1/2,-2.);
kcircle zst2=kcircle((1/6,0.),1/3,3.);
kcircle zst3=kcircle((-1/3,0.),1/6,6.);
kcircle zst4=kcircle((-3/14,2/7),1/7,7.);



void main(int cutoff = 300) {
  // main routine
  kcirclelist circlist=kcirclelist(zst1); 
  circlist.append(zst2);
  circlist.append(zst3);
  circlist.append(zst4);

  add_circs(zst1,zst2,zst3,zst4,circlist,cutoff);
  add_circs(zst2,zst3,zst4,zst1,circlist,cutoff);
  add_circs(zst3,zst4,zst1,zst2,circlist,cutoff);
  add_circs(zst4,zst1,zst2,zst3,circlist,cutoff);

  circlist.draw();              

  shipout(bbox(0.25cm, p=nullpen));  // put margins around final circle
}

main();