Lubachevsky-Stillinger Circle Packing

// Lubachevsky-Stillinger 2D circle packing.
// N disks in a box. Radii inflate continuously while overlapping pairs
// are pushed apart along their line of centers. The growth rate slows
// as the packing approaches jamming, converging to a disordered jammed
// state with packing fraction phi ~ 0.83-0.85 (random close packing in 2D).

const N = 200;
const CELL_TARGET = 28; // approximate cell size for spatial hash
let circles;            // {x,y,vx,vy,r,hue}
let R;                  // common radius (all disks identical for clean RCP target)
let growRate;
let cols, rows, cell, grid;
let W, H;
let phiTarget = 0.86;   // hard cap; growth halts above this
let jammed;
let jamTimer;

function rebuildGrid() {
  cell = Math.max(CELL_TARGET, R * 2.4);
  cols = Math.max(1, Math.ceil(W / cell));
  rows = Math.max(1, Math.ceil(H / cell));
  const need = cols * rows;
  if (!grid || grid.length !== need) grid = new Array(need);
  for (let i = 0; i < need; i++) grid[i] = null;
  for (let i = 0; i < N; i++) {
    const c = circles[i];
    const gx = Math.min(cols - 1, Math.max(0, Math.floor(c.x / cell)));
    const gy = Math.min(rows - 1, Math.max(0, Math.floor(c.y / cell)));
    const k = gy * cols + gx;
    c.next = grid[k];
    grid[k] = c;
  }
}

function init({ width, height }) {
  W = width; H = height;
  // Start radii small enough that 200 placements at random rarely collide hard.
  const area = W * H;
  R = Math.sqrt((0.05 * area) / (Math.PI * N)); // start near phi ~ 0.05
  growRate = R * 0.35;                          // px / sec at start
  jammed = false;
  jamTimer = 0;
  circles = new Array(N);
  for (let i = 0; i < N; i++) {
    const a = Math.random() * Math.PI * 2;
    const s = 60 + Math.random() * 40;
    circles[i] = {
      id: i,
      x: R + Math.random() * (W - 2 * R),
      y: R + Math.random() * (H - 2 * R),
      vx: Math.cos(a) * s,
      vy: Math.sin(a) * s,
      hue: (i * 137.5) % 360,
      next: null,
    };
  }
}

function tick({ ctx, dt, width, height, time }) {
  if (width !== W || height !== H) { W = width; H = height; }
  if (dt > 0.05) dt = 0.05;

  // 1. Grow radii (Lubachevsky-Stillinger). Slow as we approach the cap.
  const phi = (N * Math.PI * R * R) / (W * H);
  let growThisFrame = 0;
  if (phi < phiTarget) {
    const slow = Math.max(0.02, (phiTarget - phi) / phiTarget);
    growThisFrame = growRate * slow * dt;
    R += growThisFrame;
  } else if (!jammed) {
    jamTimer += dt;
    if (jamTimer > 0.6) jammed = true;
  }

  rebuildGrid();

  // 2. Advect with small thermal velocity. Bounce off walls.
  for (let i = 0; i < N; i++) {
    const c = circles[i];
    c.x += c.vx * dt;
    c.y += c.vy * dt;
    if (c.x < R) { c.x = R; c.vx = Math.abs(c.vx); }
    else if (c.x > W - R) { c.x = W - R; c.vx = -Math.abs(c.vx); }
    if (c.y < R) { c.y = R; c.vy = Math.abs(c.vy); }
    else if (c.y > H - R) { c.y = H - R; c.vy = -Math.abs(c.vy); }
  }

  // 3. Resolve overlaps. Each contact: push apart along line of centers,
  // plus a velocity exchange (relative-normal swap) to dissipate energy.
  // The push compensates for the radius growth: we add 2*growThisFrame
  // to the equality constraint so contacts can form smoothly.
  const D = 2 * R;
  const D2 = D * D;
  let contacts = 0;
  for (let i = 0; i < N; i++) {
    const c = circles[i];
    const gx = Math.min(cols - 1, Math.max(0, Math.floor(c.x / cell)));
    const gy = Math.min(rows - 1, Math.max(0, Math.floor(c.y / cell)));
    for (let oy = -1; oy <= 1; oy++) {
      for (let ox = -1; ox <= 1; ox++) {
        const nx = gx + ox, ny = gy + oy;
        if (nx < 0 || ny < 0 || nx >= cols || ny >= rows) continue;
        let o = grid[ny * cols + nx];
        while (o) {
          if (o.id > c.id) { // each pair once
            const dx = o.x - c.x, dy = o.y - c.y;
            const d2 = dx * dx + dy * dy;
            if (d2 < D2 && d2 > 1e-9) {
              const d = Math.sqrt(d2);
              const overlap = D - d + 2 * growThisFrame;
              if (overlap > 0) {
                contacts++;
                const ux = dx / d, uy = dy / d;
                const push = overlap * 0.5;
                c.x -= ux * push; c.y -= uy * push;
                o.x += ux * push; o.y += uy * push;
                // Symmetric velocity reflection along normal (elastic, equal mass).
                const rvx = o.vx - c.vx, rvy = o.vy - c.vy;
                const vn = rvx * ux + rvy * uy;
                if (vn < 0) {
                  c.vx += vn * ux; c.vy += vn * uy;
                  o.vx -= vn * ux; o.vy -= vn * uy;
                }
              }
            }
          }
          o = o.next;
        }
      }
    }
  }

  // 4. Velocity bleed so the system settles instead of buzzing at jamming.
  const damp = jammed ? 0.85 : Math.max(0.985, 1 - 0.02 * phi);
  for (let i = 0; i < N; i++) {
    const c = circles[i];
    c.vx *= damp; c.vy *= damp;
    // Clamp tiny noise when jammed.
    if (jammed) {
      if (Math.abs(c.vx) < 0.5) c.vx = 0;
      if (Math.abs(c.vy) < 0.5) c.vy = 0;
    }
  }

  // 5. Draw.
  ctx.fillStyle = "#0a0c14";
  ctx.fillRect(0, 0, W, H);

  for (let i = 0; i < N; i++) {
    const c = circles[i];
    // Color by local contact count proxy = neighbors within 1.05 D.
    let near = 0;
    const gx = Math.min(cols - 1, Math.max(0, Math.floor(c.x / cell)));
    const gy = Math.min(rows - 1, Math.max(0, Math.floor(c.y / cell)));
    const R105 = D * 1.05, R1052 = R105 * R105;
    for (let oy = -1; oy <= 1; oy++) {
      for (let ox = -1; ox <= 1; ox++) {
        const nx = gx + ox, ny = gy + oy;
        if (nx < 0 || ny < 0 || nx >= cols || ny >= rows) continue;
        let o = grid[ny * cols + nx];
        while (o) {
          if (o !== c) {
            const dx = o.x - c.x, dy = o.y - c.y;
            if (dx * dx + dy * dy < R1052) near++;
          }
          o = o.next;
        }
      }
    }
    // Isostatic 2D jamming: z = 4 contacts for frictionless disks.
    const t = Math.min(1, near / 5);
    const hue = 210 - 180 * t; // blue (low) -> orange (high)
    const light = 38 + 22 * t;
    ctx.fillStyle = `hsl(${hue.toFixed(0)},75%,${light.toFixed(0)}%)`;
    ctx.beginPath();
    ctx.arc(c.x, c.y, R, 0, Math.PI * 2);
    ctx.fill();
    ctx.strokeStyle = "rgba(0,0,0,0.35)";
    ctx.lineWidth = 1;
    ctx.stroke();
  }

  // HUD.
  const phiNow = (N * Math.PI * R * R) / (W * H);
  ctx.fillStyle = "rgba(8,10,18,0.78)";
  ctx.fillRect(8, 8, 188, 64);
  ctx.strokeStyle = "rgba(255,255,255,0.18)";
  ctx.strokeRect(8.5, 8.5, 187, 63);
  ctx.fillStyle = "#e8ecf5";
  ctx.font = "13px monospace";
  ctx.textBaseline = "top";
  ctx.fillText(`phi = ${phiNow.toFixed(4)}`, 16, 14);
  ctx.fillText(`N   = ${N}`, 16, 31);
  ctx.fillText(jammed ? "state: jammed" : "state: inflating", 16, 48);

  // Progress bar toward phi ~ 0.84 reference.
  const ref = 0.84;
  ctx.fillStyle = "rgba(255,255,255,0.12)";
  ctx.fillRect(8, 78, 188, 6);
  ctx.fillStyle = phiNow >= ref ? "#f0a050" : "#5fb0ff";
  ctx.fillRect(8, 78, Math.min(188, 188 * (phiNow / ref)), 6);
}