Brownian Motion from Molecular Kicks

// Brownian motion: heavy disk pummeled by 500 light molecules.
// Elastic collisions transfer momentum; heavy particle's path is a random walk.
// Inset plots mean-squared displacement vs time — should grow ~linearly.

const N = 500;             // number of light molecules
const R_HEAVY = 20;        // heavy particle radius
const M_HEAVY = 4000;      // heavy mass
const R_LIGHT = 1.8;       // light radius
const M_LIGHT = 1;         // light mass
const V_LIGHT = 220;       // light typical speed (px/s)
const SUBSTEPS = 2;

let W = 0, H = 0;
let lx, ly, lvx, lvy;      // light particle SoA buffers
let hx, hy, hvx, hvy;      // heavy state
let x0, y0;                // heavy start point
let msdT;                  // ring buffer of times (s)
let msdD;                  // ring buffer of squared displacement
let msdHead = 0, msdCount = 0;
const MSD_CAP = 512;
let tNow = 0;

function init({ width, height }) {
  W = width; H = height;
  lx = new Float32Array(N);
  ly = new Float32Array(N);
  lvx = new Float32Array(N);
  lvy = new Float32Array(N);
  hx = W * 0.5; hy = H * 0.5;
  hvx = 0; hvy = 0;
  x0 = hx; y0 = hy;
  for (let i = 0; i < N; i++) {
    let x, y, tries = 0;
    do {
      x = R_LIGHT + Math.random() * (W - 2 * R_LIGHT);
      y = R_LIGHT + Math.random() * (H - 2 * R_LIGHT);
      tries++;
    } while (tries < 10 && (x - hx) * (x - hx) + (y - hy) * (y - hy) < (R_HEAVY + R_LIGHT + 2) * (R_HEAVY + R_LIGHT + 2));
    lx[i] = x; ly[i] = y;
    const ang = Math.random() * Math.PI * 2;
    const s = V_LIGHT * (0.7 + Math.random() * 0.6);
    lvx[i] = Math.cos(ang) * s;
    lvy[i] = Math.sin(ang) * s;
  }
  msdT = new Float32Array(MSD_CAP);
  msdD = new Float32Array(MSD_CAP);
  msdHead = 0; msdCount = 0;
  tNow = 0;
}

function step(h) {
  // light molecules: ballistic + walls
  for (let i = 0; i < N; i++) {
    let x = lx[i] + lvx[i] * h;
    let y = ly[i] + lvy[i] * h;
    if (x < R_LIGHT) { x = R_LIGHT; lvx[i] = Math.abs(lvx[i]); }
    else if (x > W - R_LIGHT) { x = W - R_LIGHT; lvx[i] = -Math.abs(lvx[i]); }
    if (y < R_LIGHT) { y = R_LIGHT; lvy[i] = Math.abs(lvy[i]); }
    else if (y > H - R_LIGHT) { y = H - R_LIGHT; lvy[i] = -Math.abs(lvy[i]); }
    lx[i] = x; ly[i] = y;
  }

  // heavy particle: ballistic + walls
  hx += hvx * h; hy += hvy * h;
  if (hx < R_HEAVY) { hx = R_HEAVY; hvx = Math.abs(hvx); }
  else if (hx > W - R_HEAVY) { hx = W - R_HEAVY; hvx = -Math.abs(hvx); }
  if (hy < R_HEAVY) { hy = R_HEAVY; hvy = Math.abs(hvy); }
  else if (hy > H - R_HEAVY) { hy = H - R_HEAVY; hvy = -Math.abs(hvy); }

  // light <-> heavy: elastic disk-disk collision
  const Rsum = R_HEAVY + R_LIGHT;
  const R2 = Rsum * Rsum;
  const mSum = M_HEAVY + M_LIGHT;
  for (let i = 0; i < N; i++) {
    const dx = lx[i] - hx;
    const dy = ly[i] - hy;
    const d2 = dx * dx + dy * dy;
    if (d2 < R2 && d2 > 1e-6) {
      const d = Math.sqrt(d2);
      const nx = dx / d, ny = dy / d;
      // de-overlap (split by inverse mass so heavy barely moves)
      const overlap = Rsum - d;
      const wL = M_HEAVY / mSum;
      const wH = M_LIGHT / mSum;
      lx[i] += nx * overlap * wL;
      ly[i] += ny * overlap * wL;
      hx   -= nx * overlap * wH;
      hy   -= ny * overlap * wH;
      // relative velocity along normal
      const rvx = lvx[i] - hvx;
      const rvy = lvy[i] - hvy;
      const vn = rvx * nx + rvy * ny;
      if (vn < 0) {
        // elastic impulse: J = -2 vn / (1/m1 + 1/m2)
        const J = -2 * vn / (1 / M_LIGHT + 1 / M_HEAVY);
        lvx[i] += J * nx / M_LIGHT;
        lvy[i] += J * ny / M_LIGHT;
        hvx   -= J * nx / M_HEAVY;
        hvy   -= J * ny / M_HEAVY;
      }
    }
  }
}

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

  // handle clicks: teleport the heavy tracer and restart MSD tracking
  const clicks = input.consumeClicks();
  if (clicks && clicks.length) {
    const c = clicks[clicks.length - 1];
    hx = Math.min(W - R_HEAVY, Math.max(R_HEAVY, c.x));
    hy = Math.min(H - R_HEAVY, Math.max(R_HEAVY, c.y));
    hvx = 0; hvy = 0;
    x0 = hx; y0 = hy;
    msdHead = 0; msdCount = 0;
    tNow = 0;
  }

  const h = Math.min(dt, 1 / 30) / SUBSTEPS;
  for (let s = 0; s < SUBSTEPS; s++) step(h);
  tNow += Math.min(dt, 1 / 30);

  const ddx = hx - x0, ddy = hy - y0;
  const d2 = ddx * ddx + ddy * ddy;
  msdT[msdHead] = tNow; msdD[msdHead] = d2;
  msdHead = (msdHead + 1) % MSD_CAP;
  if (msdCount < MSD_CAP) msdCount++;

  ctx.fillStyle = '#070912';
  ctx.fillRect(0, 0, W, H);

  // faint molecules
  ctx.fillStyle = 'rgba(140,170,220,0.55)';
  const s = R_LIGHT * 2;
  for (let i = 0; i < N; i++) ctx.fillRect(lx[i] - R_LIGHT, ly[i] - R_LIGHT, s, s);

  // heavy particle glow + body
  const grad = ctx.createRadialGradient(hx, hy, 0, hx, hy, R_HEAVY * 2.4);
  grad.addColorStop(0, 'rgba(255,210,120,0.9)');
  grad.addColorStop(1, 'rgba(255,210,120,0)');
  ctx.fillStyle = grad;
  ctx.beginPath(); ctx.arc(hx, hy, R_HEAVY * 2.4, 0, Math.PI * 2); ctx.fill();
  ctx.fillStyle = '#ffd66b';
  ctx.beginPath(); ctx.arc(hx, hy, R_HEAVY, 0, Math.PI * 2); ctx.fill();
  ctx.strokeStyle = '#fff7d6'; ctx.lineWidth = 1.5; ctx.stroke();

  // origin cross-hair
  ctx.strokeStyle = 'rgba(255,255,255,0.35)'; ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.moveTo(x0 - 6, y0); ctx.lineTo(x0 + 6, y0);
  ctx.moveTo(x0, y0 - 6); ctx.lineTo(x0, y0 + 6);
  ctx.stroke();

  // MSD inset (bottom-right)
  const iw = Math.min(180, Math.max(120, W * 0.32));
  const ih = Math.min(110, Math.max(70, H * 0.28));
  const ix = W - iw - 8, iy = H - ih - 8;
  ctx.fillStyle = 'rgba(0,0,0,0.7)'; ctx.fillRect(ix, iy, iw, ih);
  ctx.strokeStyle = 'rgba(255,255,255,0.35)'; ctx.strokeRect(ix + 0.5, iy + 0.5, iw - 1, ih - 1);

  let dMax = 1, tMax = Math.max(1, tNow);
  for (let k = 0; k < msdCount; k++) if (msdD[k] > dMax) dMax = msdD[k];

  ctx.strokeStyle = '#7cc4ff'; ctx.lineWidth = 1.4;
  ctx.beginPath();
  for (let k = 0; k < msdCount; k++) {
    const idx = (msdHead - msdCount + k + MSD_CAP) % MSD_CAP;
    const px = ix + 4 + (iw - 8) * (msdT[idx] / tMax);
    const py = iy + ih - 4 - (ih - 12) * (msdD[idx] / dMax);
    if (k === 0) ctx.moveTo(px, py); else ctx.lineTo(px, py);
  }
  ctx.stroke();

  // linear reference line through origin & current point
  if (msdCount > 8 && tNow > 0) {
    const slope = d2 / tNow;
    ctx.strokeStyle = 'rgba(255,210,120,0.55)';
    ctx.setLineDash([3, 3]); ctx.beginPath();
    ctx.moveTo(ix + 4, iy + ih - 4);
    ctx.lineTo(ix + iw - 4, iy + ih - 4 - (ih - 12) * Math.min(1, slope * tMax / dMax));
    ctx.stroke(); ctx.setLineDash([]);
  }

  ctx.fillStyle = 'rgba(255,255,255,0.85)';
  ctx.font = '10px system-ui, sans-serif';
  ctx.fillText('MSD vs t', ix + 6, iy + 12);
  ctx.fillText(`t=${tNow.toFixed(1)}s`, ix + 6, iy + ih - 6);

  ctx.font = '12px system-ui, sans-serif';
  ctx.fillText(`N=${N} molecules  M/m=${M_HEAVY}/${M_LIGHT}`, 10, 18);
  ctx.fillText(`<r^2>/t = ${(tNow > 0 ? d2 / tNow : 0).toFixed(1)} px^2/s`, 10, 34);
}