100 Pendulums Diverging

// 100 double pendulums released from nearly-identical initial angles
// (delta-theta = 1e-6 between adjacent ones). For ~8 seconds they move as one
// ribbon, then chaos shreds them apart.
//
// Performance notes:
// - State lives in flat Float64Array (4*N), not per-pendulum objects, so the
//   RK4 kernel is a tight loop with predictable memory access.
// - Scratch buffers (k1..k4, tmp) are preallocated; no allocations per frame.
// - Same dynamics across all pendulums — the per-pendulum work is just
//   evaluating sin/cos and a few multiplies.

const N = 100;
const G = 9.81;
const L1 = 1.0, L2 = 1.0;
const M1 = 1.0, M2 = 1.0;
const DT = 0.005;       // physics step (s)
const SUBSTEPS = 6;     // → 30 ms / frame max physics work; targets 60 fps
const DELTA = 1e-6;     // angular separation between adjacent pendulums

let state, k1, k2, k3, k4, tmp;          // Float64Array(4*N)
let pivotX, pivotY, scale;
let elapsed = 0;
let energyPreset = 0;                    // 0 = high (2.0,2.0), 1 = lower

// Reusable per-frame scratch for divergence calc
let bobX, bobY;                          // Float64Array(N)

function presetAngles(preset) {
  if (preset === 1) return { th1: 1.2, th2: 1.6 }; // lower energy — slower onset
  return { th1: 2.0, th2: 2.0 };                    // default: chaos within seconds
}

function reset() {
  const { th1, th2 } = presetAngles(energyPreset);
  for (let i = 0; i < N; i++) {
    const o = i * 4;
    state[o]     = th1 + i * DELTA;   // θ1
    state[o + 1] = th2;               // θ2
    state[o + 2] = 0.0;               // ω1
    state[o + 3] = 0.0;               // ω2
  }
  elapsed = 0;
}

// Standard double-pendulum Lagrangian (same form as the existing
// double-pendulum.js sim, kept inline so the RK4 inner loop is monomorphic).
// Writes derivatives into `out` at offset `o`.
function derivsInto(src, o, out) {
  const a1 = src[o];
  const a2 = src[o + 1];
  const w1 = src[o + 2];
  const w2 = src[o + 3];
  const d = a1 - a2;
  const sd = Math.sin(d);
  const cd = Math.cos(d);
  const denomBase = (2 * M1 + M2 - M2 * Math.cos(2 * d));
  const den1 = L1 * denomBase;
  const den2 = L2 * denomBase;
  const dw1 = (-G * (2 * M1 + M2) * Math.sin(a1)
    - M2 * G * Math.sin(a1 - 2 * a2)
    - 2 * sd * M2 * (w2 * w2 * L2 + w1 * w1 * L1 * cd)) / den1;
  const dw2 = (2 * sd * (w1 * w1 * L1 * (M1 + M2)
    + G * (M1 + M2) * Math.cos(a1)
    + w2 * w2 * L2 * M2 * cd)) / den2;
  out[o]     = w1;
  out[o + 1] = w2;
  out[o + 2] = dw1;
  out[o + 3] = dw2;
}

function rk4Step(dt) {
  const len = 4 * N;
  // k1 = f(s)
  for (let i = 0; i < N; i++) derivsInto(state, i * 4, k1);
  // tmp = s + 0.5*dt*k1
  for (let i = 0; i < len; i++) tmp[i] = state[i] + 0.5 * dt * k1[i];
  for (let i = 0; i < N; i++) derivsInto(tmp, i * 4, k2);
  // tmp = s + 0.5*dt*k2
  for (let i = 0; i < len; i++) tmp[i] = state[i] + 0.5 * dt * k2[i];
  for (let i = 0; i < N; i++) derivsInto(tmp, i * 4, k3);
  // tmp = s + dt*k3
  for (let i = 0; i < len; i++) tmp[i] = state[i] + dt * k3[i];
  for (let i = 0; i < N; i++) derivsInto(tmp, i * 4, k4);
  // s += dt/6 * (k1 + 2k2 + 2k3 + k4)
  const c = dt / 6;
  for (let i = 0; i < len; i++) {
    state[i] += c * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]);
  }
}

function init({ canvas, ctx, width, height }) {
  state = new Float64Array(4 * N);
  k1    = new Float64Array(4 * N);
  k2    = new Float64Array(4 * N);
  k3    = new Float64Array(4 * N);
  k4    = new Float64Array(4 * N);
  tmp   = new Float64Array(4 * N);
  bobX  = new Float64Array(N);
  bobY  = new Float64Array(N);

  pivotX = width / 2;
  pivotY = height * 0.28;
  // Both rods extended downward span ~2*scale; want it to fit under the pivot.
  scale = Math.min(width * 0.32, height * 0.34);
  reset();
}

function tick({ ctx, dt, width, height, input }) {
  // Adapt to canvas size changes between frames.
  if (Math.abs(pivotX - width / 2) > 0.5) {
    pivotX = width / 2;
    pivotY = height * 0.28;
    scale = Math.min(width * 0.32, height * 0.34);
  }

  const clicks = input.consumeClicks();
  if (clicks.length > 0) reset();

  if (input.justPressed && input.justPressed("1")) { energyPreset = 1; reset(); }
  if (input.justPressed && input.justPressed("2")) { energyPreset = 0; reset(); }

  // Fixed substep count — physics dt is constant regardless of frame rate so
  // all 100 pendulums see the same trajectory each run. (The ribbon "reveal"
  // moment must be reproducible.)
  for (let s = 0; s < SUBSTEPS; s++) rk4Step(DT);
  elapsed += SUBSTEPS * DT;

  // ---- compute bob-2 positions ----
  let sumX = 0, sumY = 0;
  for (let i = 0; i < N; i++) {
    const o = i * 4;
    const a1 = state[o], a2 = state[o + 1];
    const s1 = Math.sin(a1), c1 = Math.cos(a1);
    const s2 = Math.sin(a2), c2 = Math.cos(a2);
    const x1 = s1 * L1 * scale;
    const y1 = c1 * L1 * scale;
    const x2 = x1 + s2 * L2 * scale;
    const y2 = y1 + c2 * L2 * scale;
    bobX[i] = x2;
    bobY[i] = y2;
    sumX += x2;
    sumY += y2;
  }
  const meanX = sumX / N, meanY = sumY / N;
  let varAcc = 0;
  for (let i = 0; i < N; i++) {
    const dx = bobX[i] - meanX;
    const dy = bobY[i] - meanY;
    varAcc += dx * dx + dy * dy;
  }
  const stddev = Math.sqrt(varAcc / N) / scale; // normalize to rod units

  // ---- render ----
  ctx.fillStyle = "rgba(7, 8, 16, 1)";
  ctx.fillRect(0, 0, width, height);

  // soft radial vignette behind the pivot for depth
  const grd = ctx.createRadialGradient(pivotX, pivotY, 0, pivotX, pivotY, Math.max(width, height) * 0.7);
  grd.addColorStop(0, "rgba(28, 22, 58, 0.55)");
  grd.addColorStop(1, "rgba(5, 5, 12, 0)");
  ctx.fillStyle = grd;
  ctx.fillRect(0, 0, width, height);

  // Alpha balances readability of the ribbon (early) vs. visibility of
  // individual strands (late). Slightly lower alpha when they're coherent.
  const ribbonAlpha = 0.55;
  ctx.lineCap = "round";
  ctx.lineJoin = "round";

  for (let i = 0; i < N; i++) {
    const o = i * 4;
    const a1 = state[o], a2 = state[o + 1];
    const x1 = pivotX + Math.sin(a1) * L1 * scale;
    const y1 = pivotY + Math.cos(a1) * L1 * scale;
    const x2 = x1 + Math.sin(a2) * L2 * scale;
    const y2 = y1 + Math.cos(a2) * L2 * scale;

    const hue = (i / N) * 360;
    ctx.strokeStyle = `hsla(${hue}, 90%, 60%, ${ribbonAlpha})`;
    ctx.lineWidth = 1.1;
    ctx.beginPath();
    ctx.moveTo(pivotX, pivotY);
    ctx.lineTo(x1, y1);
    ctx.lineTo(x2, y2);
    ctx.stroke();

    // tiny bob dot — saturated, low alpha
    ctx.fillStyle = `hsla(${hue}, 95%, 65%, 0.85)`;
    ctx.beginPath();
    ctx.arc(x2, y2, 2.2, 0, Math.PI * 2);
    ctx.fill();
  }

  // shared pivot
  ctx.fillStyle = "rgba(220, 225, 255, 1)";
  ctx.beginPath();
  ctx.arc(pivotX, pivotY, 4, 0, Math.PI * 2);
  ctx.fill();

  // ---- HUD ----
  ctx.font = "12px ui-monospace, monospace";
  ctx.fillStyle = "rgba(220, 230, 255, 0.85)";
  ctx.fillText(`t = ${elapsed.toFixed(2)} s`, 14, 22);
  ctx.fillText(`divergence σ = ${stddev.toExponential(2)} L`, 14, 40);
  ctx.fillText(`N = ${N}   Δθ₀ = 1e-6 rad`, 14, 58);

  // Phase callout: pre-divergence vs chaotic
  const phase = stddev < 0.05 ? "coherent ribbon"
              : stddev < 0.5  ? "separating"
              : "fully chaotic";
  ctx.fillStyle = stddev < 0.05 ? "hsla(200, 80%, 70%, 0.9)"
                  : stddev < 0.5 ? "hsla(40, 90%, 65%, 0.95)"
                  : "hsla(0, 85%, 65%, 0.95)";
  ctx.fillText(phase, 14, 76);

  ctx.fillStyle = "rgba(200, 210, 240, 0.5)";
  ctx.fillText("click to reset · 1/2 = lower/higher energy", 14, height - 14);
}