FHP-I Lattice Gas

// FHP-I lattice gas (Frisch-Hasslacher-Pomeau, 1986). Each hex site
// holds up to 6 particles, one per 60-degree velocity. Step = collide
// (head-on pair rotates +/- 60; (0,2,4)<->(1,3,5)), then stream along
// each particle's direction. Mass + momentum conserved -> Navier-Stokes
// in the continuum limit. We render an EMA of the local momentum field
// as hue=direction, brightness=magnitude.

const GW = 192, GH = 112, N = GW * GH;
const STEPS_PER_FRAME = 2;
const DECAY = 0.97;
const TAU = 2 * Math.PI;

// Unit velocities for the six directions. y is screen-down so we flip
// the sin term: dir 1 = up-right, dir 5 = down-right, matching the NB
// table below.
const VX = new Float32Array(6), VY = new Float32Array(6);
for (let i = 0; i < 6; i++) { VX[i] = Math.cos(i*TAU/6); VY[i] = -Math.sin(i*TAU/6); }

// Offset-row hex neighbor table. NB[parity][dir] = [dx, dy].
const NB = [
  [[ 1, 0], [ 0,-1], [-1,-1], [-1, 0], [-1, 1], [ 0, 1]],
  [[ 1, 0], [ 1,-1], [ 0,-1], [-1, 0], [ 0, 1], [ 1, 1]],
];
const BIT = [1, 2, 4, 8, 16, 32];
const ALL = 63;

// Obstacle (circular cylinder) in lattice coords. OX/OY are mutable so
// the user can drag the cylinder around — see tick() for the input read.
let OX = GW * 0.32, OY = GH * 0.5;
const OR = Math.min(GW, GH) * 0.10, OR2 = OR * OR;
function isObs(x, y) { const dx=x-OX, dy=y-OY; return dx*dx+dy*dy <= OR2; }

let cur, nxt, avgVx, avgVy, off, octx, img, pix, stepCount;

const idx = (x, y) => y * GW + x;

// Two collision tables (random rotation choice for head-on pairs).
const CA = new Uint8Array(64), CB = new Uint8Array(64);
function buildTables() {
  for (let s = 0; s < 64; s++) { CA[s] = s; CB[s] = s; }
  for (const [a, b] of [[0,3],[1,4],[2,5]]) {
    const s = BIT[a] | BIT[b];
    CA[s] = BIT[(a+1)%6] | BIT[(b+1)%6];
    CB[s] = BIT[(a+5)%6] | BIT[(b+5)%6];
  }
  CA[BIT[0]|BIT[2]|BIT[4]] = BIT[1]|BIT[3]|BIT[5];
  CB[BIT[0]|BIT[2]|BIT[4]] = BIT[1]|BIT[3]|BIT[5];
  CA[BIT[1]|BIT[3]|BIT[5]] = BIT[0]|BIT[2]|BIT[4];
  CB[BIT[1]|BIT[3]|BIT[5]] = BIT[0]|BIT[2]|BIT[4];
}

function seed() {
  for (let i = 0; i < N; i++) {
    const y = (i / GW) | 0, x = i - y * GW;
    let s = 0;
    if (!isObs(x, y)) for (let d = 0; d < 6; d++) if (Math.random() < 0.30) s |= BIT[d];
    cur[i] = s;
  }
  avgVx.fill(0); avgVy.fill(0);
}

function collide() {
  for (let i = 0; i < N; i++) {
    cur[i] = ((i + stepCount) & 1) ? CA[cur[i]] : CB[cur[i]];
  }
}

function stream() {
  nxt.fill(0);
  for (let y = 0; y < GH; y++) {
    const par = y & 1, tab = NB[par];
    for (let x = 0; x < GW; x++) {
      const s = cur[idx(x, y)];
      if (s === 0) continue;
      if (isObs(x, y)) {
        let r = 0;
        for (let d = 0; d < 6; d++) if (s & BIT[d]) r |= BIT[(d+3)%6];
        nxt[idx(x, y)] |= r;
        continue;
      }
      for (let d = 0; d < 6; d++) {
        if (!(s & BIT[d])) continue;
        const nb = tab[d];
        let nx = x + nb[0], ny = y + nb[1];
        if (nx < 0) nx += GW; else if (nx >= GW) nx -= GW;
        if (ny < 0 || ny >= GH || isObs(nx, ny)) {
          nxt[idx(x, y)] |= BIT[(d+3)%6];
        } else {
          nxt[idx(nx, ny)] |= BIT[d];
        }
      }
    }
  }
  const t = cur; cur = nxt; nxt = t;
}

// Inject eastward flow at the inlet column.
function pump() {
  for (let y = 1; y < GH - 1; y++) {
    const i = idx(0, y);
    cur[i] |= BIT[0];
    if (Math.random() < 0.4) cur[i] |= BIT[1];
    if (Math.random() < 0.4) cur[i] |= BIT[5];
    cur[i] &= ~(BIT[2] | BIT[3] | BIT[4]);
  }
}

function accumulate() {
  for (let i = 0; i < N; i++) {
    const s = cur[i];
    let mx = 0, my = 0;
    if (s) for (let d = 0; d < 6; d++) if (s & BIT[d]) { mx += VX[d]; my += VY[d]; }
    avgVx[i] = DECAY * avgVx[i] + (1 - DECAY) * mx;
    avgVy[i] = DECAY * avgVy[i] + (1 - DECAY) * my;
  }
}

function render() {
  for (let y = 0; y < GH; y++) {
    for (let x = 0; x < GW; x++) {
      const i = idx(x, y), p = i * 4;
      if (isObs(x, y)) { pix[p]=30; pix[p+1]=30; pix[p+2]=38; pix[p+3]=255; continue; }
      const vx = avgVx[i], vy = avgVy[i];
      const m = Math.min(1, Math.sqrt(vx*vx + vy*vy) * 1.6);
      const hue = (Math.atan2(vy, vx) / TAU + 1) % 1;
      const v = 0.18 + 0.82 * m, s = 0.85;
      const h6 = hue * 6, c = v * s, xh = c * (1 - Math.abs((h6 % 2) - 1)), mm = v - c;
      let rr = 0, gg = 0, bb = 0;
      if (h6 < 1)      { rr = c; gg = xh; }
      else if (h6 < 2) { rr = xh; gg = c; }
      else if (h6 < 3) { gg = c; bb = xh; }
      else if (h6 < 4) { gg = xh; bb = c; }
      else if (h6 < 5) { rr = xh; bb = c; }
      else             { rr = c; bb = xh; }
      pix[p]   = ((rr + mm) * 255) | 0;
      pix[p+1] = ((gg + mm) * 255) | 0;
      pix[p+2] = ((bb + mm) * 255) | 0;
      pix[p+3] = 255;
    }
  }
  octx.putImageData(img, 0, 0);
}

function init({ width, height }) {
  cur = new Uint8Array(N);
  nxt = new Uint8Array(N);
  avgVx = new Float32Array(N);
  avgVy = new Float32Array(N);
  off = new OffscreenCanvas(GW, GH);
  octx = off.getContext("2d");
  img = octx.createImageData(GW, GH);
  pix = img.data;
  stepCount = 0;
  buildTables();
  seed();
  for (let k = 0; k < 30; k++) { pump(); collide(); stream(); accumulate(); stepCount++; }
  render();
}

function tick({ ctx, width, height, input }) {
  const scale = Math.min(width / GW, height / GH);
  const drawW = GW * scale, drawH = GH * scale;
  const ox = (width - drawW) * 0.5, oy = (height - drawH) * 0.5;

  // Drag the obstacle around: while mouse is held, snap OX/OY to the
  // mouse position mapped into lattice coords, clamped so the cylinder
  // stays fully inside the grid.
  if (input && input.mouseDown) {
    const lx = (input.mouseX - ox) / scale;
    const ly = (input.mouseY - oy) / scale;
    const margin = OR + 1;
    OX = Math.max(margin, Math.min(GW - margin, lx));
    OY = Math.max(margin, Math.min(GH - margin, ly));
  }

  for (let k = 0; k < STEPS_PER_FRAME; k++) {
    pump(); collide(); stream(); accumulate(); stepCount++;
  }
  render();

  ctx.fillStyle = "#0a0a10";
  ctx.fillRect(0, 0, width, height);
  ctx.imageSmoothingEnabled = true;
  ctx.drawImage(off, 0, 0, GW, GH, ox, oy, drawW, drawH);

  ctx.fillStyle = "rgba(0,0,0,0.55)";
  ctx.fillRect(8, 8, 200, 46);
  ctx.fillStyle = "#e8e8ee";
  ctx.font = "12px sans-serif";
  ctx.fillText(`FHP-I lattice gas  ${GW}x${GH} hex`, 14, 24);
  ctx.fillText(`step ${stepCount}   hue = mean flow dir`, 14, 40);
}