Maxwell-Boltzmann Thermalization

// Maxwell-Boltzmann thermalization in 2D.
// 200 equal-mass hard disks start at the SAME speed in random directions.
// Pairwise elastic collisions redistribute KE until the speed distribution
// matches the analytic 2D Maxwell-Boltzmann P(v) = (v/kT) exp(-v^2/2kT)
// with kT read off from <KE>/N (m=1, k_B=1).

const N = 200, R = 4, V0 = 140, SUBSTEPS = 3, BINS = 32, HIST_VMAX = 360;
const HIST_DECAY = 0.985;

let W = 0, H = 0, simW = 0, simH = 0, plotX = 0, plotW = 0;
let px, py, vx, vy;
let head, nxt, cellSize, cols, rows;
let hist, frameHist;
let collisions = 0;

function layout() {
  const split = W >= 720 ? 0.58 : 0.5;
  simW = Math.max(120, Math.floor(W * split) - 8);
  simH = H;
  plotX = simW + 8;
  plotW = Math.max(80, W - plotX - 8);
}

function init({ width, height }) {
  W = width; H = height; layout();
  px = new Float32Array(N); py = new Float32Array(N);
  vx = new Float32Array(N); vy = new Float32Array(N);
  const c0 = Math.ceil(Math.sqrt(N * simW / Math.max(simH, 1)));
  const r0 = Math.ceil(N / c0);
  const sx = simW / (c0 + 1), sy = simH / (r0 + 1);
  let k = 0;
  for (let j = 0; j < r0 && k < N; j++) for (let i = 0; i < c0 && k < N; i++) {
    px[k] = sx * (i + 1) + (Math.random() - 0.5) * Math.min(sx, sy) * 0.3;
    py[k] = sy * (j + 1) + (Math.random() - 0.5) * Math.min(sx, sy) * 0.3;
    const a = Math.random() * Math.PI * 2;
    vx[k] = Math.cos(a) * V0; vy[k] = Math.sin(a) * V0;
    k++;
  }
  cellSize = R * 2.4;
  cols = Math.max(1, Math.ceil(simW / cellSize));
  rows = Math.max(1, Math.ceil(simH / cellSize));
  head = new Int32Array(cols * rows);
  nxt = new Int32Array(N);
  hist = new Float32Array(BINS);
  frameHist = new Float32Array(BINS);
  collisions = 0;
}

function collide(i, j) {
  const dx = px[j] - px[i], dy = py[j] - py[i];
  const d2 = dx * dx + dy * dy;
  const Rs = 2 * R, R2 = Rs * Rs;
  if (d2 >= R2 || d2 < 1e-8) return;
  const d = Math.sqrt(d2), nx = dx / d, ny = dy / d;
  const half = (Rs - d) * 0.5;
  px[i] -= nx * half; py[i] -= ny * half;
  px[j] += nx * half; py[j] += ny * half;
  const vn = (vx[j] - vx[i]) * nx + (vy[j] - vy[i]) * ny;
  if (vn >= 0) return;
  vx[i] += vn * nx; vy[i] += vn * ny;
  vx[j] -= vn * nx; vy[j] -= vn * ny;
  collisions++;
}

const NB_DX = [1, -1, 0, 1], NB_DY = [0, 1, 1, 1];

function step(h) {
  for (let i = 0; i < N; i++) {
    let x = px[i] + vx[i] * h, y = py[i] + vy[i] * h;
    if (x < R) { x = R; vx[i] = Math.abs(vx[i]); }
    else if (x > simW - R) { x = simW - R; vx[i] = -Math.abs(vx[i]); }
    if (y < R) { y = R; vy[i] = Math.abs(vy[i]); }
    else if (y > simH - R) { y = simH - R; vy[i] = -Math.abs(vy[i]); }
    px[i] = x; py[i] = y;
  }
  for (let c = 0; c < head.length; c++) head[c] = -1;
  for (let i = 0; i < N; i++) {
    const cx = Math.min(cols - 1, Math.max(0, Math.floor(px[i] / cellSize)));
    const cy = Math.min(rows - 1, Math.max(0, Math.floor(py[i] / cellSize)));
    const c = cy * cols + cx;
    nxt[i] = head[c]; head[c] = i;
  }
  for (let cy = 0; cy < rows; cy++) for (let cx = 0; cx < cols; cx++) {
    const c = cy * cols + cx;
    for (let i = head[c]; i !== -1; i = nxt[i]) {
      for (let j = nxt[i]; j !== -1; j = nxt[j]) collide(i, j);
      for (let k = 0; k < 4; k++) {
        const ax = cx + NB_DX[k], ay = cy + NB_DY[k];
        if (ax < 0 || ax >= cols || ay < 0 || ay >= rows) continue;
        for (let j = head[ay * cols + ax]; j !== -1; j = nxt[j]) collide(i, j);
      }
    }
  }
}

function meanKE() {
  let s = 0;
  for (let i = 0; i < N; i++) s += vx[i] * vx[i] + vy[i] * vy[i];
  return 0.5 * s / N;
}

function buildHist() {
  for (let b = 0; b < BINS; b++) frameHist[b] = 0;
  const dv = HIST_VMAX / BINS;
  for (let i = 0; i < N; i++) {
    let b = Math.floor(Math.hypot(vx[i], vy[i]) / dv);
    if (b >= BINS) b = BINS - 1;
    frameHist[b]++;
  }
  for (let b = 0; b < BINS; b++) {
    const pdf = frameHist[b] / (N * dv);
    hist[b] = hist[b] * HIST_DECAY + pdf * (1 - HIST_DECAY);
  }
}

function speedColor(spd) {
  const t = Math.min(1, spd / (V0 * 2.5));
  return `hsl(${(220 * (1 - t)).toFixed(0)},85%,60%)`;
}

function drawHist(ctx) {
  const x = plotX, y = 8, w = plotW, h = simH - 16;
  if (w < 60) return;
  ctx.fillStyle = 'rgba(0,0,0,0.55)'; ctx.fillRect(x, y, w, h);
  ctx.strokeStyle = 'rgba(255,255,255,0.25)';
  ctx.strokeRect(x + 0.5, y + 0.5, w - 1, h - 1);

  const kT = Math.max(1e-3, meanKE());
  const samples = 96, dv = HIST_VMAX / samples;
  const curve = new Float32Array(samples + 1);
  let pdfMax = 0;
  for (let i = 0; i <= samples; i++) {
    const v = i * dv;
    const p = (v / kT) * Math.exp(-(v * v) / (2 * kT));
    curve[i] = p; if (p > pdfMax) pdfMax = p;
  }
  let yMax = pdfMax;
  for (let b = 0; b < BINS; b++) if (hist[b] > yMax) yMax = hist[b];
  yMax *= 1.15;

  const padL = 28, padR = 8, padT = 22, padB = 22;
  const ax = x + padL, ay = y + padT, aw = w - padL - padR, ah = h - padT - padB;
  ctx.strokeStyle = 'rgba(255,255,255,0.3)';
  ctx.beginPath();
  ctx.moveTo(ax, ay); ctx.lineTo(ax, ay + ah); ctx.lineTo(ax + aw, ay + ah);
  ctx.stroke();

  const binDv = HIST_VMAX / BINS;
  ctx.fillStyle = 'rgba(255,180,90,0.55)';
  for (let b = 0; b < BINS; b++) {
    const x0 = ax + (b * binDv / HIST_VMAX) * aw;
    const x1 = ax + ((b + 1) * binDv / HIST_VMAX) * aw;
    const bh = Math.min(ah, (hist[b] / yMax) * ah);
    ctx.fillRect(x0 + 1, ay + ah - bh, Math.max(1, x1 - x0 - 2), bh);
  }

  ctx.strokeStyle = '#7cc4ff'; ctx.lineWidth = 1.6;
  ctx.beginPath();
  for (let i = 0; i <= samples; i++) {
    const cx = ax + (i * dv / HIST_VMAX) * aw;
    const cy = ay + ah - (curve[i] / yMax) * ah;
    if (i === 0) ctx.moveTo(cx, cy); else ctx.lineTo(cx, cy);
  }
  ctx.stroke();

  const v0x = ax + (V0 / HIST_VMAX) * aw;
  ctx.strokeStyle = 'rgba(255,255,255,0.35)';
  ctx.setLineDash([3, 3]);
  ctx.beginPath(); ctx.moveTo(v0x, ay); ctx.lineTo(v0x, ay + ah); ctx.stroke();
  ctx.setLineDash([]);

  ctx.fillStyle = 'rgba(255,255,255,0.85)';
  ctx.font = '11px system-ui, sans-serif'; ctx.textAlign = 'left';
  ctx.fillText('P(v)  vs  v', x + 8, y + 14);
  ctx.textAlign = 'center';
  ctx.fillText('v (px/s)', x + w / 2, y + h - 4);
  ctx.textAlign = 'left';
  ctx.fillStyle = '#7cc4ff'; ctx.fillText('— analytic', x + w - 92, y + 14);
  ctx.fillStyle = 'rgba(255,180,90,0.9)'; ctx.fillText('▮ measured', x + w - 92, y + 28);
  ctx.fillStyle = 'rgba(255,255,255,0.75)';
  ctx.font = '10px system-ui, sans-serif'; ctx.textAlign = 'center';
  ctx.fillText('v₀', v0x, ay - 4);
}

function tick({ ctx, dt, width, height }) {
  if (width !== W || height !== H) { W = width; H = height; layout(); }
  const h = Math.min(dt, 1 / 30) / SUBSTEPS;
  for (let s = 0; s < SUBSTEPS; s++) step(h);
  buildHist();

  ctx.fillStyle = '#070912'; ctx.fillRect(0, 0, W, H);
  ctx.strokeStyle = 'rgba(255,255,255,0.18)'; ctx.lineWidth = 1;
  ctx.strokeRect(0.5, 0.5, simW - 1, simH - 1);
  for (let i = 0; i < N; i++) {
    ctx.fillStyle = speedColor(Math.hypot(vx[i], vy[i]));
    ctx.beginPath(); ctx.arc(px[i], py[i], R, 0, Math.PI * 2); ctx.fill();
  }
  drawHist(ctx);

  const kT = meanKE();
  ctx.fillStyle = 'rgba(0,0,0,0.55)'; ctx.fillRect(8, 8, 188, 50);
  ctx.fillStyle = '#fff';
  ctx.font = '12px monospace'; ctx.textAlign = 'left'; ctx.textBaseline = 'alphabetic';
  ctx.fillText(`N=${N}  collisions=${collisions}`, 14, 24);
  ctx.fillText(`kT = <KE>/N = ${kT.toFixed(0)}`, 14, 42);
}