Kuramoto Oscillators

const N = 100;
const K_MAX = 6.0;
const DT_SIM = 0.02;
const SUBSTEPS = 2;
const TWO_PI = Math.PI * 2;

let theta;
let omega;
let sinBuf;
let cosBuf;
let K;
let kSmooth;
let rHist;
let rHistIdx;
let W;
let H;

function gaussian() {
  // Box-Muller — natural frequencies ω_i from N(0, 1).
  let u = 0;
  let v = 0;
  while (u === 0) u = Math.random();
  while (v === 0) v = Math.random();
  return Math.sqrt(-2 * Math.log(u)) * Math.cos(TWO_PI * v);
}

function init({ ctx, width, height }) {
  W = width;
  H = height;

  theta = new Float32Array(N);
  omega = new Float32Array(N);
  sinBuf = new Float32Array(N);
  cosBuf = new Float32Array(N);
  for (let i = 0; i < N; i++) {
    theta[i] = Math.random() * TWO_PI;
    omega[i] = gaussian();
  }

  K = 1.5;
  kSmooth = 1.5;
  rHist = new Float32Array(240);
  rHistIdx = 0;

  ctx.fillStyle = '#06080d';
  ctx.fillRect(0, 0, W, H);
}

function step(dtSim, k) {
  // Mean-field form: dθ_i/dt = ω_i + K · ( sin(ψ - θ_i) · r )
  // where r e^{iψ} = (1/N) Σ e^{iθ_j}. Avoids the O(N^2) double-sum.
  let sx = 0;
  let sy = 0;
  for (let i = 0; i < N; i++) {
    sx += Math.cos(theta[i]);
    sy += Math.sin(theta[i]);
  }
  const cx = sx / N;
  const cy = sy / N;
  const r = Math.sqrt(cx * cx + cy * cy);
  const psi = Math.atan2(cy, cx);

  for (let i = 0; i < N; i++) {
    const dth = omega[i] + k * r * Math.sin(psi - theta[i]);
    theta[i] += dth * dtSim;
    if (theta[i] > TWO_PI) theta[i] -= TWO_PI;
    else if (theta[i] < 0) theta[i] += TWO_PI;
  }
  return r;
}

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

  // Mouse X scrubs coupling K across [0, K_MAX]. Default to 1.5 when offscreen.
  let targetK = kSmooth;
  if (input && input.mouseX >= 0 && input.mouseX <= W) {
    targetK = (input.mouseX / W) * K_MAX;
  }
  kSmooth += (targetK - kSmooth) * Math.min(1, dt * 6);
  K = kSmooth;

  let r = 0;
  const steps = Math.max(1, Math.min(8, Math.round(SUBSTEPS * (dt / (1 / 60)))));
  for (let s = 0; s < steps; s++) {
    r = step(DT_SIM, K);
  }

  rHist[rHistIdx] = r;
  rHistIdx = (rHistIdx + 1) % rHist.length;

  // Fade prior frame for a subtle motion trail.
  ctx.globalCompositeOperation = 'source-over';
  ctx.fillStyle = 'rgba(6, 8, 13, 0.55)';
  ctx.fillRect(0, 0, W, H);

  const cx = W * 0.5;
  const cy = H * 0.5;
  const radius = Math.min(W, H) * 0.36;

  // Reference circle.
  ctx.strokeStyle = 'rgba(120, 140, 180, 0.25)';
  ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.arc(cx, cy, radius, 0, TWO_PI);
  ctx.stroke();

  // Oscillators — hue from natural frequency ω, position on unit circle by θ.
  ctx.globalCompositeOperation = 'lighter';
  for (let i = 0; i < N; i++) {
    const x = cx + Math.cos(theta[i]) * radius;
    const y = cy + Math.sin(theta[i]) * radius;
    const hue = 200 + omega[i] * 50; // fast → red, slow → blue
    ctx.fillStyle = `hsla(${hue.toFixed(0)}, 90%, 65%, 0.9)`;
    ctx.beginPath();
    ctx.arc(x, y, 3.5, 0, TWO_PI);
    ctx.fill();
  }
  ctx.globalCompositeOperation = 'source-over';

  // Order parameter arrow: r e^{iψ}.
  let sx = 0;
  let sy = 0;
  for (let i = 0; i < N; i++) {
    sx += Math.cos(theta[i]);
    sy += Math.sin(theta[i]);
  }
  const mx = sx / N;
  const my = sy / N;
  const rNow = Math.sqrt(mx * mx + my * my);
  const psi = Math.atan2(my, mx);
  const ax = cx + Math.cos(psi) * radius * rNow;
  const ay = cy + Math.sin(psi) * radius * rNow;

  ctx.strokeStyle = 'rgba(255, 220, 100, 0.85)';
  ctx.lineWidth = 2.2;
  ctx.beginPath();
  ctx.moveTo(cx, cy);
  ctx.lineTo(ax, ay);
  ctx.stroke();
  ctx.fillStyle = 'rgba(255, 230, 140, 0.95)';
  ctx.beginPath();
  ctx.arc(ax, ay, 4, 0, TWO_PI);
  ctx.fill();

  // r-history strip along the bottom.
  const stripH = 40;
  const stripY = H - stripH - 6;
  const stripX = 12;
  const stripW = W - 24;
  ctx.fillStyle = 'rgba(20, 26, 38, 0.7)';
  ctx.fillRect(stripX, stripY, stripW, stripH);
  ctx.strokeStyle = 'rgba(120, 140, 180, 0.35)';
  ctx.strokeRect(stripX + 0.5, stripY + 0.5, stripW - 1, stripH - 1);
  ctx.strokeStyle = 'rgba(255, 220, 100, 0.85)';
  ctx.lineWidth = 1.5;
  ctx.beginPath();
  const L = rHist.length;
  for (let i = 0; i < L; i++) {
    const idx = (rHistIdx + i) % L;
    const v = rHist[idx];
    const px = stripX + (i / (L - 1)) * stripW;
    const py = stripY + stripH - v * stripH;
    if (i === 0) ctx.moveTo(px, py);
    else ctx.lineTo(px, py);
  }
  ctx.stroke();

  // HUD.
  ctx.fillStyle = 'rgba(220, 230, 245, 0.95)';
  ctx.font = '13px ui-monospace, monospace';
  ctx.textBaseline = 'top';
  ctx.fillText(`K = ${K.toFixed(2)}`, 12, 10);
  ctx.fillText(`r = ${rNow.toFixed(3)}`, 12, 28);
  ctx.fillText(`N = ${N}`, 12, 46);

  ctx.fillStyle = 'rgba(160, 175, 200, 0.8)';
  ctx.font = '11px ui-monospace, monospace';
  ctx.textAlign = 'right';
  ctx.fillText('move cursor → scrub K', W - 12, 10);
  ctx.fillText(rNow > 0.7 ? 'synchronized' : rNow > 0.3 ? 'partial lock' : 'incoherent', W - 12, 26);
  ctx.textAlign = 'left';
}