Belousov-Zhabotinsky Reaction-Diffusion

// Belousov-Zhabotinsky reaction-diffusion CA on a 200x150 grid.
// Three coupled concentrations A, B, C with cyclic activator-inhibitor
// kinetics: A activates B, B activates C, C suppresses A. Each species
// diffuses through a 9-point stencil (self + 8 neighbors averaged) so
// local imbalances propagate outward as travelling fronts. The cyclic
// coupling forces those fronts to chase one another, breaking into
// spirals and target patterns characteristic of BZ chemistry.
// Render: additive RGB — R from A, G from B, B from C.

const GW = 200, GH = 150;
const N = GW * GH;

let A, B, C, A2, B2, C2;
let off, octx, img, pix;
let acc = 0;
const stepMs = 1000 / 45;

function seedSpiral(cx, cy, r) {
  // Stamp a disc of high-A, low-B/C cells to nucleate a fresh travelling wave.
  for (let y = -r; y <= r; y++) {
    for (let x = -r; x <= r; x++) {
      if (x * x + y * y > r * r) continue;
      const xi = ((cx + x) % GW + GW) % GW;
      const yi = ((cy + y) % GH + GH) % GH;
      const i = yi * GW + xi;
      A[i] = 1.0;
      B[i] = 0.0;
      C[i] = 0.0;
    }
  }
}

function seed() {
  for (let i = 0; i < N; i++) {
    A[i] = Math.random();
    B[i] = Math.random();
    C[i] = Math.random();
  }
  // a few hot spots to bias spiral nucleation
  for (let k = 0; k < 6; k++) {
    const cx = (Math.random() * GW) | 0;
    const cy = (Math.random() * GH) | 0;
    const r = 4 + ((Math.random() * 6) | 0);
    seedSpiral(cx, cy, r);
  }
}

function init({ ctx }) {
  A = new Float32Array(N);
  B = new Float32Array(N);
  C = new Float32Array(N);
  A2 = new Float32Array(N);
  B2 = new Float32Array(N);
  C2 = new Float32Array(N);
  seed();

  off = new OffscreenCanvas(GW, GH);
  octx = off.getContext("2d");
  img = octx.createImageData(GW, GH);
  pix = img.data;
  for (let i = 3; i < pix.length; i += 4) pix[i] = 255;

  // pre-steps so first frame has visible structure (kinetics are gentle now)
  for (let s = 0; s < 40; s++) step();
}

function step() {
  // Cyclic BZ-like kinetics. A feeds B, B feeds C, C inhibits A.
  // Local diffusion via 9-point average (Moore neighborhood).
  // Tuned for stable rotating spirals (vs chaotic flicker): gentle reaction
  // rates relative to diffusion, small dt so the explicit Euler step never
  // overshoots the [0,1] clamp.
  const kAB = 0.45;  // A -> B conversion strength when B present
  const kBC = 0.45;  // B -> C conversion strength when C present
  const kCA = 0.45;  // C decay catalysed by A absence
  const fA  = 0.03;  // A feed
  const dA  = 0.02;  // A natural decay
  const dB  = 0.035; // B natural decay
  const dC  = 0.035; // C natural decay
  const diff = 0.45; // fraction of value replaced by neighborhood mean
  const dt  = 0.3;

  for (let y = 0; y < GH; y++) {
    const ym = (y - 1 + GH) % GH;
    const yp = (y + 1) % GH;
    const row  = y  * GW;
    const rowm = ym * GW;
    const rowp = yp * GW;
    for (let x = 0; x < GW; x++) {
      const xm = (x - 1 + GW) % GW;
      const xp = (x + 1) % GW;
      const i = row + x;

      // Moore-neighborhood means (8 neighbors averaged).
      const aN = (A[rowm + xm] + A[rowm + x] + A[rowm + xp] +
                  A[row  + xm]               + A[row  + xp] +
                  A[rowp + xm] + A[rowp + x] + A[rowp + xp]) * 0.125;
      const bN = (B[rowm + xm] + B[rowm + x] + B[rowm + xp] +
                  B[row  + xm]               + B[row  + xp] +
                  B[rowp + xm] + B[rowp + x] + B[rowp + xp]) * 0.125;
      const cN = (C[rowm + xm] + C[rowm + x] + C[rowm + xp] +
                  C[row  + xm]               + C[row  + xp] +
                  C[rowp + xm] + C[rowp + x] + C[rowp + xp]) * 0.125;

      const a = A[i] * (1 - diff) + aN * diff;
      const b = B[i] * (1 - diff) + bN * diff;
      const c = C[i] * (1 - diff) + cN * diff;

      // Reactions: cyclic activator → next species, last suppresses first.
      const rAB = kAB * a * b;     // A consumed, B produced
      const rBC = kBC * b * c;     // B consumed, C produced
      const rCA = kCA * c * (1 - a); // C consumed when A is low

      let na = a + dt * (fA * (1 - a) - rAB - dA * a + 0.3 * rCA);
      let nb = b + dt * (rAB - rBC - dB * b);
      let nc = c + dt * (rBC - rCA - dC * c);

      if (na < 0) na = 0; else if (na > 1) na = 1;
      if (nb < 0) nb = 0; else if (nb > 1) nb = 1;
      if (nc < 0) nc = 0; else if (nc > 1) nc = 1;

      A2[i] = na;
      B2[i] = nb;
      C2[i] = nc;
    }
  }
  const ta = A; A = A2; A2 = ta;
  const tb = B; B = B2; B2 = tb;
  const tc = C; C = C2; C2 = tc;
}

function tick({ ctx, dt, width, height, input }) {
  // Click-to-seed: each click stamps a small high-A disc at the click
  // location, nucleating a new travelling wave / spiral.
  const clicks = input.consumeClicks();
  if (clicks && clicks.length) {
    for (let k = 0; k < clicks.length; k++) {
      const cl = clicks[k];
      const cx = Math.round((cl.x / width) * GW);
      const cy = Math.round((cl.y / height) * GH);
      seedSpiral(cx, cy, 6);
    }
  }

  acc += Math.min(0.05, dt) * 1000;
  let n = 0;
  while (acc >= stepMs && n < 4) { step(); acc -= stepMs; n++; }

  // Additive RGB from concentrations A, B, C.
  for (let i = 0; i < N; i++) {
    const j = i << 2;
    let r = A[i] * 320;
    let g = B[i] * 320;
    let b = C[i] * 320;
    if (r > 255) r = 255;
    if (g > 255) g = 255;
    if (b > 255) b = 255;
    pix[j]     = r;
    pix[j + 1] = g;
    pix[j + 2] = b;
  }
  octx.putImageData(img, 0, 0);
  ctx.imageSmoothingEnabled = false;
  ctx.drawImage(off, 0, 0, width, height);
}