Maxwell–Boltzmann Speed Distribution

// 300 gas particles in a box, elastic wall bounces (and lightweight
// particle-particle elastic collisions on a coarse spatial hash).
// Build a live histogram of speeds; overlay the analytic Maxwell-Boltzmann
// distribution f(v) ∝ v^2 exp(-m v^2 / (2 k_B T)) at the temperature
// inferred from the current mean kinetic energy. Press H to heat
// (multiply speeds by 1.15), C to cool (divide), R to reset.

const N = 300;
const R = 3;
const NBINS = 40;
const VMAX_DISPLAY = 600;   // px/s — top of histogram axis

let W = 0, H = 0;
let px, py, vx, vy;
let bins;
let cellSize, cellsX, cellsY, cellHead, cellNext;
const INITIAL_SIGMA = 140; // px/s — sets initial temperature

function init({ width, height }) {
  W = width; H = height;
  px = new Float32Array(N);
  py = new Float32Array(N);
  vx = new Float32Array(N);
  vy = new Float32Array(N);
  // Maxwell-Boltzmann initial: 2D Gaussian velocity components
  // (so |v| follows Rayleigh — already MB in 2D).
  const sigma = INITIAL_SIGMA;
  for (let i = 0; i < N; i++) {
    px[i] = R + Math.random() * (W - 2 * R);
    py[i] = R + Math.random() * (H - 2 * R);
    // Box-Muller
    const u1 = Math.max(1e-9, Math.random());
    const u2 = Math.random();
    const u3 = Math.max(1e-9, Math.random());
    const u4 = Math.random();
    vx[i] = sigma * Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2);
    vy[i] = sigma * Math.sqrt(-2 * Math.log(u3)) * Math.cos(2 * Math.PI * u4);
  }
  bins = new Float32Array(NBINS);
  cellSize = R * 4;
}

function rescale(factor) {
  for (let i = 0; i < N; i++) { vx[i] *= factor; vy[i] *= factor; }
}

function step(h) {
  // ballistic + walls
  for (let i = 0; i < N; i++) {
    let x = px[i] + vx[i] * h;
    let y = py[i] + vy[i] * h;
    if (x < R) { x = R; vx[i] = Math.abs(vx[i]); }
    else if (x > W - R) { x = W - R; vx[i] = -Math.abs(vx[i]); }
    if (y < R) { y = R; vy[i] = Math.abs(vy[i]); }
    else if (y > H - R) { y = H - R; vy[i] = -Math.abs(vy[i]); }
    px[i] = x; py[i] = y;
  }

  // pair collisions via uniform grid
  cellsX = Math.max(1, Math.ceil(W / cellSize));
  cellsY = Math.max(1, Math.ceil(H / cellSize));
  const nc = cellsX * cellsY;
  if (!cellHead || cellHead.length !== nc) cellHead = new Int32Array(nc);
  if (!cellNext || cellNext.length !== N) cellNext = new Int32Array(N);
  cellHead.fill(-1);
  for (let i = 0; i < N; i++) {
    const cx = Math.min(cellsX - 1, Math.max(0, (px[i] / cellSize) | 0));
    const cy = Math.min(cellsY - 1, Math.max(0, (py[i] / cellSize) | 0));
    const c = cy * cellsX + cx;
    cellNext[i] = cellHead[c];
    cellHead[c] = i;
  }
  const R2 = (2 * R) * (2 * R);
  for (let cy = 0; cy < cellsY; cy++) {
    for (let cx = 0; cx < cellsX; cx++) {
      const c = cy * cellsX + cx;
      for (let dy = 0; dy <= 1; dy++) {
        for (let dx = -1; dx <= 1; dx++) {
          if (dy === 0 && dx < 0) continue;
          const nx2 = cx + dx, ny2 = cy + dy;
          if (nx2 < 0 || nx2 >= cellsX || ny2 < 0 || ny2 >= cellsY) continue;
          const c2 = ny2 * cellsX + nx2;
          for (let i = cellHead[c]; i !== -1; i = cellNext[i]) {
            for (let j = (c === c2) ? cellNext[i] : cellHead[c2]; j !== -1; j = cellNext[j]) {
              const ddx = px[i] - px[j], ddy = py[i] - py[j];
              const d2 = ddx * ddx + ddy * ddy;
              if (d2 < R2 && d2 > 1e-6) {
                const d = Math.sqrt(d2);
                const nxn = ddx / d, nyn = ddy / d;
                // de-overlap equally
                const overlap = (2 * R - d) * 0.5;
                px[i] += nxn * overlap; py[i] += nyn * overlap;
                px[j] -= nxn * overlap; py[j] -= nyn * overlap;
                const rvx = vx[i] - vx[j], rvy = vy[i] - vy[j];
                const vn = rvx * nxn + rvy * nyn;
                if (vn < 0) {
                  vx[i] -= vn * nxn; vy[i] -= vn * nyn;
                  vx[j] += vn * nxn; vy[j] += vn * nyn;
                }
              }
            }
          }
        }
      }
    }
  }
}

function tick({ ctx, dt, width, height, input }) {
  if (width !== W || height !== H) { W = width; H = height; }
  const h = Math.min(dt, 1 / 30);
  step(h);

  if (input.justPressed('h') || input.justPressed('H') ||
      input.justPressed('ArrowUp')) rescale(1.15);
  if (input.justPressed('c') || input.justPressed('C') ||
      input.justPressed('ArrowDown')) rescale(1 / 1.15);
  if (input.justPressed('r') || input.justPressed('R')) init({ width: W, height: H });

  // touch/click controls: top band heats, bottom band cools, middle resets
  const clicks = input.consumeClicks ? input.consumeClicks() : [];
  for (let k = 0; k < clicks.length; k++) {
    const cy = clicks[k].y;
    if (cy < H * 0.4) rescale(1.06);
    else if (cy > H * 0.6) rescale(0.95);
    else init({ width: W, height: H });
  }

  // measure mean kinetic energy (m=1) => <v^2> = 2 k_B T / m for 2D
  // so k_B T = <v^2>/2. In our units we just use a "T-proxy" = <v^2>/2.
  let v2sum = 0;
  let vmax = 0;
  bins.fill(0);
  for (let i = 0; i < N; i++) {
    const sp = Math.hypot(vx[i], vy[i]);
    if (sp > vmax) vmax = sp;
    v2sum += sp * sp;
    const b = Math.min(NBINS - 1, (sp / VMAX_DISPLAY * NBINS) | 0);
    bins[b] += 1;
  }
  const meanV2 = v2sum / N;
  const kT = 0.5 * meanV2;          // k_B*T in our units (m=1, 2D)
  const vRMS = Math.sqrt(meanV2);
  const vMostProb = Math.sqrt(kT);  // 2D: v_p = sqrt(k_B T / m)
  const vMean = Math.sqrt(Math.PI * kT / 2);  // 2D mean speed

  // background
  ctx.fillStyle = '#0c0f17';
  ctx.fillRect(0, 0, W, H);

  // particle area = top half
  const topH = Math.floor(H * 0.58);
  ctx.fillStyle = '#070912';
  ctx.fillRect(0, 0, W, topH);

  // colour by speed
  for (let i = 0; i < N; i++) {
    const sp = Math.hypot(vx[i], vy[i]);
    const t = Math.min(1, sp / VMAX_DISPLAY);
    const r8 = (60 + 195 * t) | 0;
    const g8 = (160 - 80 * t) | 0;
    const b8 = (255 - 200 * t) | 0;
    ctx.fillStyle = `rgb(${r8},${g8},${b8})`;
    const yy = px[i]; // unused — clarity only
    if (py[i] < topH) {
      ctx.beginPath();
      ctx.arc(px[i], py[i], R, 0, Math.PI * 2);
      ctx.fill();
    } else {
      // wrap into the box (we constrain motion to top portion)
      py[i] = topH - R - 1; vy[i] = -Math.abs(vy[i]);
    }
  }

  // confine motion: lower wall at topH
  for (let i = 0; i < N; i++) {
    if (py[i] > topH - R) { py[i] = topH - R; vy[i] = -Math.abs(vy[i]); }
  }

  // histogram region
  const hx0 = 50, hy0 = topH + 30;
  const hw  = W - 80;
  const hh  = H - hy0 - 40;
  ctx.fillStyle = 'rgba(0,0,0,0.55)';
  ctx.fillRect(hx0 - 6, hy0 - 6, hw + 12, hh + 18);

  // bar bins (normalised so the histogram has unit area)
  let binMax = 0;
  for (let b = 0; b < NBINS; b++) if (bins[b] > binMax) binMax = bins[b];
  // analytic MB max we'll overlay; scale both to same vertical max
  // 2D MB pdf: f(v) = (m / kT) v exp(-m v^2 / (2 kT))   for m=1
  const dv = VMAX_DISPLAY / NBINS;
  // peak of analytic pdf is at v = sqrt(kT), value = sqrt(1/(e*kT))
  const pdfPeak = 1 / Math.sqrt(Math.E * Math.max(1e-6, kT));
  // we want histogram bar to match: counts_b ~ N * pdf(v_b) * dv
  // so scale-pdf-to-bin-units by N*dv
  const yScale = (hh - 4) / Math.max(binMax, N * pdfPeak * dv * 1.05);

  for (let b = 0; b < NBINS; b++) {
    const v = (b + 0.5) * dv;
    const bw = hw / NBINS;
    const barH = bins[b] * yScale;
    const t = Math.min(1, v / VMAX_DISPLAY);
    const r8 = (60 + 195 * t) | 0;
    const g8 = (160 - 80 * t) | 0;
    const b8 = (255 - 200 * t) | 0;
    ctx.fillStyle = `rgba(${r8},${g8},${b8},0.7)`;
    ctx.fillRect(hx0 + b * bw, hy0 + hh - barH, bw - 1, barH);
  }

  // analytic MB curve overlay (2D)
  ctx.strokeStyle = '#ffffff';
  ctx.lineWidth = 1.8;
  ctx.beginPath();
  for (let i = 0; i <= 100; i++) {
    const v = (i / 100) * VMAX_DISPLAY;
    const pdf = (1 / Math.max(1e-6, kT)) * v * Math.exp(-v * v / (2 * Math.max(1e-6, kT)));
    const yv = N * pdf * dv * yScale;
    const xp = hx0 + (v / VMAX_DISPLAY) * hw;
    const yp = hy0 + hh - yv;
    if (i === 0) ctx.moveTo(xp, yp); else ctx.lineTo(xp, yp);
  }
  ctx.stroke();

  // axes / ticks
  ctx.strokeStyle = 'rgba(255,255,255,0.4)';
  ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.moveTo(hx0, hy0); ctx.lineTo(hx0, hy0 + hh); ctx.lineTo(hx0 + hw, hy0 + hh);
  ctx.stroke();
  ctx.fillStyle = 'rgba(255,255,255,0.7)';
  ctx.font = '10px system-ui, sans-serif';
  for (let k = 0; k <= 5; k++) {
    const v = (k / 5) * VMAX_DISPLAY;
    const xp = hx0 + (k / 5) * hw;
    ctx.fillRect(xp, hy0 + hh, 1, 4);
    ctx.fillText(v.toFixed(0), xp - 8, hy0 + hh + 14);
  }
  ctx.fillText('speed |v| (px/s)', hx0 + hw / 2 - 50, hy0 + hh + 28);

  // markers
  function marker(v, color, label) {
    const xp = hx0 + (v / VMAX_DISPLAY) * hw;
    ctx.strokeStyle = color;
    ctx.setLineDash([3, 3]);
    ctx.beginPath();
    ctx.moveTo(xp, hy0); ctx.lineTo(xp, hy0 + hh);
    ctx.stroke();
    ctx.setLineDash([]);
    ctx.fillStyle = color;
    ctx.fillText(label, xp + 3, hy0 + 11);
  }
  marker(vMostProb, '#7cffb2', `v_p=${vMostProb.toFixed(0)}`);
  marker(vMean,     '#ffd66b', `<v>=${vMean.toFixed(0)}`);
  marker(vRMS,      '#ff8e8e', `v_rms=${vRMS.toFixed(0)}`);

  // header
  ctx.fillStyle = 'rgba(0,0,0,0.6)';
  ctx.fillRect(0, 0, W, 26);
  ctx.fillStyle = '#ffffff';
  ctx.font = 'bold 13px system-ui, sans-serif';
  ctx.fillText(`N=${N} particles   k_B T = <v²>/2 = ${kT.toFixed(0)}   (m=1)`,
    10, 18);
  ctx.font = '11px system-ui, sans-serif';
  ctx.fillStyle = 'rgba(255,255,255,0.75)';
  ctx.fillText('H = heat  ·  C = cool  ·  R = reset', W - 200, 18);

  // touch hints — small labels along the left edge for tap zones
  ctx.font = '10px system-ui, sans-serif';
  ctx.fillStyle = 'rgba(255,180,140,0.85)';
  ctx.fillText('tap top to heat', 8, H * 0.4 - 6);
  ctx.fillStyle = 'rgba(255,255,255,0.7)';
  ctx.fillText('tap middle to reset', 8, H * 0.5 + 3);
  ctx.fillStyle = 'rgba(140,200,255,0.85)';
  ctx.fillText('tap bottom to cool', 8, H * 0.6 + 12);
}