Three-Species Trophic Cascade

// Three-species generalized Lotka-Volterra food chain: grass -> rabbits -> foxes.
// Left panel: scrolling time series of G, R, F.
// Right panel: oblique-projected 3D phase trajectory in (G,R,F) space.

const R_G = 1.0;     // grass intrinsic growth
const K   = 1.0;     // grass carrying capacity
const A_GR = 1.2;    // grass -> rabbit grazing rate
const E_GR = 0.7;    // grass-to-rabbit conversion efficiency
const A_RF = 1.0;    // rabbit -> fox predation rate
const E_RF = 0.5;    // rabbit-to-fox conversion efficiency
const D_R  = 0.4;    // rabbit death rate
const D_F_MIN = 0.18;
const D_F_MAX = 0.55;
let   D_F  = 0.3;    // fox death rate (mouseY-controlled)

const DT_SIM = 0.01;
const SUBSTEPS = 4;

const TS_MAX = 720;   // time-series buffer (frames worth of samples)
const PHASE_MAX = 3000;

let W, H;
let state;            // [G, R, F]
let ts;               // Float32Array, 3*TS_MAX
let tsHead, tsCount;
let phase;            // Float32Array, 3*PHASE_MAX
let phHead, phCount;
let simTime;
let lastPerturbT;
let perturbFlash;

// Scratch slots for rk4 — inlined to avoid per-step allocations.
let _k1G = 0, _k1R = 0, _k1F = 0;
let _k2G = 0, _k2R = 0, _k2F = 0;
let _k3G = 0, _k3R = 0, _k3F = 0;
let _k4G = 0, _k4R = 0, _k4F = 0;

function _setDerivs(G, R, F, target) {
  const dG = R_G * G * (1 - G / K) - A_GR * G * R;
  const dR = E_GR * A_GR * G * R - A_RF * R * F - D_R * R;
  const dF = E_RF * A_RF * R * F - D_F * F;
  if (target === 1) { _k1G = dG; _k1R = dR; _k1F = dF; }
  else if (target === 2) { _k2G = dG; _k2R = dR; _k2F = dF; }
  else if (target === 3) { _k3G = dG; _k3R = dR; _k3F = dF; }
  else { _k4G = dG; _k4R = dR; _k4F = dF; }
}

function rk4Step(s, h) {
  const G = s[0], R = s[1], F = s[2];
  _setDerivs(G, R, F, 1);
  _setDerivs(G + 0.5 * h * _k1G, R + 0.5 * h * _k1R, F + 0.5 * h * _k1F, 2);
  _setDerivs(G + 0.5 * h * _k2G, R + 0.5 * h * _k2R, F + 0.5 * h * _k2F, 3);
  _setDerivs(G + h * _k3G,       R + h * _k3R,       F + h * _k3F,       4);
  let nG = G + (h / 6) * (_k1G + 2 * _k2G + 2 * _k3G + _k4G);
  let nR = R + (h / 6) * (_k1R + 2 * _k2R + 2 * _k3R + _k4R);
  let nF = F + (h / 6) * (_k1F + 2 * _k2F + 2 * _k3F + _k4F);
  if (nG < 0) nG = 0;
  if (nR < 0) nR = 0;
  if (nF < 0) nF = 0;
  s[0] = nG; s[1] = nR; s[2] = nF;
  return s;
}

function pushTS(G, R, F) {
  ts[tsHead * 3]     = G;
  ts[tsHead * 3 + 1] = R;
  ts[tsHead * 3 + 2] = F;
  tsHead = (tsHead + 1) % TS_MAX;
  if (tsCount < TS_MAX) tsCount++;
}

function pushPhase(G, R, F) {
  phase[phHead * 3]     = G;
  phase[phHead * 3 + 1] = R;
  phase[phHead * 3 + 2] = F;
  phHead = (phHead + 1) % PHASE_MAX;
  if (phCount < PHASE_MAX) phCount++;
}

function init({ width, height }) {
  W = width;
  H = height;

  state = [0.6, 0.3, 0.2];
  ts = new Float32Array(TS_MAX * 3);
  tsHead = 0; tsCount = 0;
  phase = new Float32Array(PHASE_MAX * 3);
  phHead = 0; phCount = 0;
  simTime = 0;
  lastPerturbT = -1e9;
  perturbFlash = 0;

  // burn in so we land on/near the attractor
  for (let i = 0; i < 800; i++) {
    rk4Step(state, DT_SIM);
    if (i % 4 === 0) pushPhase(state[0], state[1], state[2]);
  }
}

// Project (G, R, F) into 2D using a simple oblique projection.
// X = G + 0.45*F,  Y = R + 0.30*F  (then mapped into panel pixel space).
// Writes result into module-scoped _projX/_projY to avoid per-call allocations.
let _projX = 0, _projY = 0;
function project(G, R, F, panelX, panelY, panelW, panelH) {
  const sX = G + 0.45 * F;
  const sY = R + 0.30 * F;
  const SX_MAX = 1.7;
  const SY_MAX = 1.7;
  _projX = panelX + 30 + (sX / SX_MAX) * (panelW - 50);
  _projY = panelY + panelH - 30 - (sY / SY_MAX) * (panelH - 50);
}

function drawAxes3D(ctx, panelX, panelY, panelW, panelH) {
  ctx.strokeStyle = '#2a3344';
  ctx.lineWidth = 1;
  // origin
  project(0, 0, 0, panelX, panelY, panelW, panelH); const ox = _projX, oy = _projY;
  // G axis -> (1,0,0)
  project(1, 0, 0, panelX, panelY, panelW, panelH); const gx = _projX, gy = _projY;
  // R axis -> (0,1,0)
  project(0, 1, 0, panelX, panelY, panelW, panelH); const rx = _projX, ry = _projY;
  // F axis -> (0,0,1)
  project(0, 0, 1, panelX, panelY, panelW, panelH); const fx = _projX, fy = _projY;

  ctx.beginPath();
  ctx.moveTo(ox, oy); ctx.lineTo(gx, gy);
  ctx.moveTo(ox, oy); ctx.lineTo(rx, ry);
  ctx.moveTo(ox, oy); ctx.lineTo(fx, fy);
  ctx.stroke();

  ctx.font = '11px monospace';
  ctx.fillStyle = '#7fdc7f'; ctx.fillText('G', gx + 4, gy + 4);
  ctx.fillStyle = '#f0a040'; ctx.fillText('R', rx - 12, ry - 2);
  ctx.fillStyle = '#d04848'; ctx.fillText('F', fx + 4, fy - 4);
}

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

  // mouseY scrubs fox death rate
  if (input && typeof input.mouseY === 'number' && input.mouseY >= 0) {
    const t = Math.max(0, Math.min(1, input.mouseY / H));
    // top of screen = low death rate (predators thrive), bottom = high death rate
    D_F = D_F_MIN + (1 - t) * (D_F_MAX - D_F_MIN);
  }

  // clicks cull 50% of foxes
  if (input && input.consumeClicks) {
    const clicks = input.consumeClicks();
    if (clicks && clicks.length > 0) {
      state[2] *= 0.5;
      lastPerturbT = simTime;
      perturbFlash = 1.0;
    }
  }

  // integrate
  const stepsPerFrame = Math.max(1, Math.min(8, Math.round(SUBSTEPS * (dt / (1 / 60)))));
  for (let i = 0; i < stepsPerFrame; i++) {
    rk4Step(state, DT_SIM);
    simTime += DT_SIM;
  }
  pushTS(state[0], state[1], state[2]);
  pushPhase(state[0], state[1], state[2]);

  // bg
  ctx.fillStyle = '#0a0e14';
  ctx.fillRect(0, 0, W, H);

  // layout
  const margin = 14;
  const leftW = Math.floor((W - margin * 3) * 0.55);
  const rightW = W - margin * 3 - leftW;
  const PX_L = margin, PY_L = 36, PW_L = leftW, PH_L = H - 56;
  const PX_R = margin * 2 + leftW, PY_R = 36, PW_R = rightW, PH_R = H - 56;

  // headline title
  ctx.fillStyle = '#cfd6e0';
  ctx.font = '13px monospace';
  ctx.fillText('Three-Species Trophic Cascade', margin, 20);

  // === LEFT PANEL: scrolling time series ===
  ctx.strokeStyle = '#1a2030';
  ctx.strokeRect(PX_L, PY_L, PW_L, PH_L);

  // gridlines
  ctx.strokeStyle = '#141a26';
  ctx.lineWidth = 1;
  for (let i = 1; i < 4; i++) {
    const yy = PY_L + (i / 4) * PH_L;
    ctx.beginPath();
    ctx.moveTo(PX_L, yy);
    ctx.lineTo(PX_L + PW_L, yy);
    ctx.stroke();
  }

  // y-axis label
  ctx.fillStyle = '#6a7388';
  ctx.font = '10px monospace';
  ctx.fillText('1.5', PX_L + 4, PY_L + 11);
  ctx.fillText('0',   PX_L + 4, PY_L + PH_L - 4);
  ctx.fillText('time →', PX_L + PW_L - 56, PY_L + PH_L + 14);

  const Y_MAX = 1.5;
  const sx = (i, n) => PX_L + 4 + ((i) / (TS_MAX - 1)) * (PW_L - 8);
  const sy = (v) => PY_L + PH_L - 4 - Math.min(1, v / Y_MAX) * (PH_L - 14);

  if (tsCount > 1) {
    const startIdx = (tsHead - tsCount + TS_MAX) % TS_MAX;
    const drawSeries = (offset, color) => {
      ctx.strokeStyle = color;
      ctx.lineWidth = 1.6;
      ctx.beginPath();
      for (let i = 0; i < tsCount; i++) {
        const idx = (startIdx + i) % TS_MAX;
        const v = ts[idx * 3 + offset];
        const px = sx(i + (TS_MAX - tsCount), TS_MAX);
        const py = sy(v);
        if (i === 0) ctx.moveTo(px, py); else ctx.lineTo(px, py);
      }
      ctx.stroke();
    };
    drawSeries(0, '#7fdc7f');  // grass
    drawSeries(1, '#f0a040');  // rabbits
    drawSeries(2, '#d04848');  // foxes
  }

  // legend (top-right of left panel)
  const legendX = PX_L + PW_L - 110;
  const legendY = PY_L + 8;
  ctx.font = '11px monospace';
  ctx.fillStyle = '#7fdc7f'; ctx.fillRect(legendX, legendY, 10, 2); ctx.fillText('grass G', legendX + 14, legendY + 4);
  ctx.fillStyle = '#f0a040'; ctx.fillRect(legendX, legendY + 12, 10, 2); ctx.fillText('rabbits R', legendX + 14, legendY + 16);
  ctx.fillStyle = '#d04848'; ctx.fillRect(legendX, legendY + 24, 10, 2); ctx.fillText('foxes F',  legendX + 14, legendY + 28);

  // current values + dF readout at bottom of left panel
  ctx.font = '11px monospace';
  ctx.fillStyle = '#7fdc7f'; ctx.fillText('G=' + state[0].toFixed(3), PX_L + 8, PY_L + PH_L - 28);
  ctx.fillStyle = '#f0a040'; ctx.fillText('R=' + state[1].toFixed(3), PX_L + 100, PY_L + PH_L - 28);
  ctx.fillStyle = '#d04848'; ctx.fillText('F=' + state[2].toFixed(3), PX_L + 192, PY_L + PH_L - 28);
  ctx.fillStyle = '#cfd6e0'; ctx.fillText('d_F = ' + D_F.toFixed(3), PX_L + 8, PY_L + PH_L - 12);

  // === RIGHT PANEL: 3D phase trajectory ===
  ctx.strokeStyle = '#1a2030';
  ctx.strokeRect(PX_R, PY_R, PW_R, PH_R);
  ctx.fillStyle = '#6a7388';
  ctx.font = '10px monospace';
  ctx.fillText('phase space  (G, R, F)', PX_R + 8, PY_R + 14);

  drawAxes3D(ctx, PX_R, PY_R, PW_R, PH_R);

  if (phCount > 1) {
    const startIdx = (phHead - phCount + PHASE_MAX) % PHASE_MAX;
    ctx.globalCompositeOperation = 'lighter';
    ctx.lineWidth = 1.2;
    ctx.lineCap = 'round';

    let prevPX = -1, prevPY = -1;
    for (let i = 0; i < phCount; i++) {
      const idx = (startIdx + i) % PHASE_MAX;
      const G = phase[idx * 3];
      const R = phase[idx * 3 + 1];
      const F = phase[idx * 3 + 2];
      project(G, R, F, PX_R, PY_R, PW_R, PH_R);
      const px = _projX, py = _projY;
      if (i > 0) {
        const t = i / phCount;
        const alpha = 0.04 + t * 0.55;
        // colour by which species dominates locally
        const denom = G + R + F + 1e-6;
        const cr = Math.round(80 + 175 * (F / denom));
        const cg = Math.round(80 + 175 * (G / denom));
        const cb = Math.round(80 + 130 * (R / denom));
        ctx.strokeStyle = `rgba(${cr},${cg},${cb},${alpha.toFixed(3)})`;
        ctx.beginPath();
        ctx.moveTo(prevPX, prevPY);
        ctx.lineTo(px, py);
        ctx.stroke();
      }
      prevPX = px; prevPY = py;
    }

    // current point
    ctx.globalCompositeOperation = 'source-over';
    project(state[0], state[1], state[2], PX_R, PY_R, PW_R, PH_R);
    const cpx = _projX, cpy = _projY;
    ctx.fillStyle = 'rgba(255,230,140,0.25)';
    ctx.beginPath(); ctx.arc(cpx, cpy, 10, 0, Math.PI * 2); ctx.fill();
    ctx.fillStyle = '#ffe680';
    ctx.beginPath(); ctx.arc(cpx, cpy, 3.4, 0, Math.PI * 2); ctx.fill();
  }
  ctx.globalCompositeOperation = 'source-over';

  // perturbation flash overlay on right panel
  if (perturbFlash > 0) {
    ctx.fillStyle = `rgba(208, 72, 72, ${(perturbFlash * 0.18).toFixed(3)})`;
    ctx.fillRect(PX_R, PY_R, PW_R, PH_R);
    perturbFlash = Math.max(0, perturbFlash - dt * 1.5);
  }

  // footer hint
  ctx.fillStyle = '#6a7388';
  ctx.font = '11px monospace';
  ctx.fillText('click: cull 50% of foxes  ·  mouse Y: tune fox death rate', margin, H - 8);
}