Apollonian Gasket

// Apollonian Gasket — Descartes Circle Theorem, depth-by-depth growth.
// Curvature k = 1/r (negative for the enclosing circle).
// k4 = k1+k2+k3 +/- 2*sqrt(k1*k2 + k2*k3 + k1*k3)
// In complex centers, with b_i = k_i*z_i:
//   b4 = b1+b2+b3 +/- 2*sqrt(b1*b2 + b2*b3 + b1*b3)

let W, H, circles, frontier, depth, accum, lastSpawn;
const MAX_DEPTH = 9;

function cAdd(a, b) { return { x: a.x + b.x, y: a.y + b.y }; }
function cSub(a, b) { return { x: a.x - b.x, y: a.y - b.y }; }
function cMul(a, b) { return { x: a.x * b.x - a.y * b.y, y: a.x * b.y + a.y * b.x }; }
function cScale(a, s) { return { x: a.x * s, y: a.y * s }; }
function cSqrt(a) {
  const r = Math.hypot(a.x, a.y);
  const sgn = a.y < 0 ? -1 : 1;
  return { x: Math.sqrt((r + a.x) * 0.5), y: sgn * Math.sqrt((r - a.x) * 0.5) };
}

function mk(k, z, d) { return { k, z, r: Math.abs(1 / k), d }; }

function tangent(a, b) {
  const d = Math.hypot(a.z.x - b.z.x, a.z.y - b.z.y);
  if (a.k < 0 || b.k < 0) return Math.abs(d - Math.abs(a.r - b.r)) < 0.6;
  return Math.abs(d - (a.r + b.r)) < 0.6;
}

function startGasket() {
  const cx = W * 0.5, cy = H * 0.5;
  const R = Math.min(W, H) * 0.46;
  const C0 = mk(-1 / R, { x: cx, y: cy }, 0);

  // Three equal mutually tangent inner circles inscribed in C0:
  // r_inner = R * sqrt(3) / (sqrt(3) + 2)
  const ri = R * Math.sqrt(3) / (Math.sqrt(3) + 2);
  const dC = R - ri;
  const inners = [];
  for (let i = 0; i < 3; i++) {
    const a = -Math.PI / 2 + (i * 2 * Math.PI) / 3;
    inners.push(mk(1 / ri, { x: cx + Math.cos(a) * dC, y: cy + Math.sin(a) * dC }, 1));
  }
  circles = [C0, ...inners];

  // Initial tangent triples to recurse on.
  frontier = [
    [1, 2, 3],       // three inners
    [0, 1, 2], [0, 2, 3], [0, 1, 3], // outer + pair of inners
  ];
  depth = 1;
}

function descartesNew(c1, c2, c3) {
  const k1 = c1.k, k2 = c2.k, k3 = c3.k;
  const sum = k1 + k2 + k3;
  const rad = k1 * k2 + k2 * k3 + k1 * k3;
  if (rad < 0) return null;
  const sq = Math.sqrt(rad);
  const b1 = cScale(c1.z, k1), b2 = cScale(c2.z, k2), b3 = cScale(c3.z, k3);
  const bsum = cAdd(cAdd(b1, b2), b3);
  const inside = cAdd(cAdd(cMul(b1, b2), cMul(b2, b3)), cMul(b1, b3));
  const root = cSqrt(inside);
  const cands = [
    { k: sum + 2 * sq, b: cAdd(bsum, cScale(root, 2)) },
    { k: sum + 2 * sq, b: cSub(bsum, cScale(root, 2)) },
    { k: sum - 2 * sq, b: cAdd(bsum, cScale(root, 2)) },
    { k: sum - 2 * sq, b: cSub(bsum, cScale(root, 2)) },
  ];
  let best = null;
  for (const cd of cands) {
    if (!isFinite(cd.k) || cd.k <= 0) continue;
    const r = 1 / cd.k;
    if (r < 0.6) continue;
    const z = cScale(cd.b, 1 / cd.k);
    const nc = { k: cd.k, z, r, d: 0 };
    if (!tangent(nc, c1) || !tangent(nc, c2) || !tangent(nc, c3)) continue;
    if (!best || r < best.r) best = nc;
  }
  return best;
}

function spawnDepth() {
  const next = [];
  const added = [];
  for (const [a, b, c] of frontier) {
    const nc = descartesNew(circles[a], circles[b], circles[c]);
    if (!nc) continue;
    let dup = false;
    for (let i = 0; i < circles.length; i++) {
      const ex = circles[i];
      if (Math.abs(ex.r - nc.r) < 0.6 && Math.hypot(ex.z.x - nc.z.x, ex.z.y - nc.z.y) < 0.6) { dup = true; break; }
    }
    if (dup) continue;
    nc.d = depth + 1;
    const idx = circles.length + added.length;
    added.push(nc);
    next.push([a, b, idx], [a, c, idx], [b, c, idx]);
  }
  for (const nc of added) circles.push(nc);
  frontier = next.filter(t => Math.min(circles[t[0]].r, circles[t[1]].r, circles[t[2]].r) > 1.2);
  depth++;
}

function depthColor(d) {
  const h = (d * 53 + 200) % 360;
  return `hsl(${h} 80% 60%)`;
}

function drawScene(ctx) {
  ctx.fillStyle = "#070914";
  ctx.fillRect(0, 0, W, H);
  const outer = circles[0];
  ctx.beginPath();
  ctx.arc(outer.z.x, outer.z.y, outer.r, 0, Math.PI * 2);
  ctx.fillStyle = "#141a34";
  ctx.fill();

  const order = circles.map((_, i) => i).sort((i, j) => {
    const a = circles[i], b = circles[j];
    return a.d !== b.d ? a.d - b.d : b.r - a.r;
  });
  for (const i of order) {
    const c = circles[i];
    if (c.r < 0.8) continue;
    ctx.beginPath();
    ctx.arc(c.z.x, c.z.y, c.r, 0, Math.PI * 2);
    if (i === 0) {
      ctx.strokeStyle = "rgba(220,230,255,0.9)";
      ctx.lineWidth = 2;
      ctx.stroke();
    } else {
      ctx.fillStyle = depthColor(c.d);
      ctx.globalAlpha = Math.max(0.2, 0.9 - c.d * 0.07);
      ctx.fill();
      ctx.globalAlpha = 1;
      ctx.strokeStyle = "rgba(0,0,0,0.55)";
      ctx.lineWidth = c.r > 30 ? 1.5 : 1;
      ctx.stroke();
    }
  }
}

function drawHUD(ctx) {
  ctx.fillStyle = "rgba(0,0,0,0.55)";
  ctx.fillRect(10, 10, 230, 56);
  ctx.fillStyle = "#fff";
  ctx.font = "13px ui-monospace, monospace";
  ctx.fillText("depth: " + Math.min(depth, MAX_DEPTH) + " / " + MAX_DEPTH, 18, 30);
  ctx.fillText("circles: " + circles.length, 18, 50);
}

function init({ width, height }) {
  W = width; H = height;
  accum = 0;
  lastSpawn = 0;
  startGasket();
  spawnDepth();
}

function tick({ ctx, dt, width, height }) {
  if (width !== W || height !== H) {
    W = width; H = height;
    startGasket();
    spawnDepth();
    accum = 0;
    lastSpawn = 0;
  }
  accum += dt || 0.016;
  if (depth < MAX_DEPTH && accum - lastSpawn > 0.55) {
    spawnDepth();
    lastSpawn = accum;
  }
  if (depth >= MAX_DEPTH && accum - lastSpawn > 3.5) {
    startGasket();
    spawnDepth();
    accum = 0;
    lastSpawn = 0;
  }
  drawScene(ctx);
  drawHUD(ctx);
}