On a world Turing

I was recently one of the luck 10,000 when I stumbled across these implementations of a Reaction-Diffusion system in APL & BQN. Investigating more, it turns out these ideas date back at least as far as Alan Turing. Karl Sims’s particularly useful tutorial describes the Gray-Scott model. In this, the system describes the interactions of two “chemicals” \(A\) and \(B\) according to the following pair of differential equations:

\[\frac{\partial A}{\partial t} = D_A\nabla^2A - AB^2 + f(1 - A)\] \[\frac{\partial B}{\partial t} = D_B\nabla^2B + AB^2 - (k + f)B \]

There are a number of great resources all over the web covering different ways of modeling this in code.

Demo

Borrowing some ideas from the aforementioned sources, here’s a little demo implemented in Javascript. Adjust the parameters to your heart’s content, then hit play/stop & reset as desired!

Code plz

My Javascript code is short enough to include here; even unminified, it’s only 75 lines after formatting with deno.¹

"use strict";
const N = 256;
const A = new Float64Array(N * N), B = new Float64Array(N * N);
const canvas = document.getElementById("canvas");
canvas.width = canvas.height = N;
const ctx = canvas.getContext("2d");
const pixels = new ImageData(N, N);
const buttons = {
  reset: document.getElementById("reset"),
  render: document.getElementById("render"),
};
let running = false;
const splitmix = (seed) => {
  let state = BigInt(seed);
  const mask = (1n << 64n) - 1n;
  return () => {
    state = (state + 0x9e3779b97f4a7c15n) & mask;
    let z = state;
    z = (z ^ (z >> 30n)) * 0xbf58476d1ce4e5b9n & mask;
    z = (z ^ (z >> 27n)) * 0x94d049bb133111ebn & mask;
    return Number(z ^ (z >> 31n) & mask) / 2**64;
  };
};
const init = (seed) => {
  const N2 = Math.floor(N / 2), R = Math.floor(N / 10);
  const rng = splitmix(seed);
  for (let x = 0; x < N; x++) {
    for (let y = 0; y < N; y++) {
      const i = x + y * N;
      if (x >= N2 - R && x < N2 + R && y >= N2 - R && y < N2 + R) {
        A[i] = 0.5;
        B[i] = 0.25;
      } else {
        A[i] = 1 - (B[i] = 0.25 * rng());
      }
    }
  }
  draw();
};
const laplacian = (m, i) => {
  return (m[i - N] ?? 0) + (m[i + 1] ?? 0) + (m[i + N] ?? 0) + (m[i - 1] ?? 0) -
    4 * m[i];
};
const step = (da, db, f, k) => {
  for (let x = 0; x < N; x++) {
    for (let y = 0; y < N; y++) {
      const i = x + y * N;
      const reaction = A[i] * B[i] * B[i];
      A[i] += da * laplacian(A, i) - reaction + f * (1 - A[i]);
      B[i] += db * laplacian(B, i) + reaction - (k + f) * B[i];
    }
  }
};
const draw = () => {
  for (let i = 0; i < N * N; i++) {
    const v = 16 * Math.pow(B[i], 2);
    pixels.data[4 * i] = Math.floor(2 * v);
    pixels.data[4 * i + 1] = Math.floor(127 * v);
    pixels.data[4 * i + 2] = Math.floor(141 * v);
    pixels.data[4 * i + 3] = 255;
  }
  ctx.putImageData(pixels, 0, 0);
};
const loop = () => {
  step(Number(da.value), Number(db.value), Number(feed.value), Number(kill.value));
  draw();
  if (running) requestAnimationFrame(loop);
};
init(Number(seed.value));
buttons.reset.addEventListener("click", () => { init(Number(seed.value)); });
buttons.render.addEventListener("click", () => {
  running = !running;
  buttons.render.innerHTML = running ? "Stop" : "Play";
  if (running) loop();
});