Galaxy Rotation: Visible Mass vs Dark Matter

// Galaxy rotation curves: Newtonian (visible-mass) prediction vs observed flat curve.
// Two face-on spiral galaxies side-by-side, each with ~125 stars in concentric rings.
// Left galaxy spins each star at v_newt(r); right galaxy at v_obs(r). Below them,
// the two v(r) curves are plotted continuously so the disagreement is visible.

const N_STARS = 250;          // total across both galaxies
const N_PER = N_STARS / 2;
const N_RINGS = 14;
const R_MAX = 1.0;            // sim units, scaled to pixel radius per panel
const R_BULGE = 0.18;         // bulge scale radius
const R_DISK = 0.32;          // exponential disk scale length

// Visible-mass profile M_vis(r): bulge (Hernquist-ish) + exponential disk.
// Normalized so M_vis(R_MAX) ~ 1.
function mVis(r) {
  const mBulge = (r * r) / Math.pow(r + R_BULGE, 2);            // 0..1, saturates fast
  // Exponential disk enclosed mass: 1 - (1 + r/Rd) * exp(-r/Rd).
  const u = r / R_DISK;
  const mDisk = 1 - (1 + u) * Math.exp(-u);
  return 0.35 * mBulge + 0.65 * mDisk;
}

// Newtonian orbital speed from visible mass: v_n(r) = sqrt(M_vis(r) / r).
// Add a small softening to avoid divide-by-zero at r=0.
function vNewt(r) {
  const rs = Math.max(r, 0.02);
  return Math.sqrt(mVis(rs) / rs);
}

// Observed flat rotation: rises through the bulge, then flat plateau.
// v_obs(r) = V_FLAT * (r / sqrt(r^2 + R_CORE^2)).
const V_FLAT = 1.05;
const R_CORE = 0.12;
function vObs(r) {
  return V_FLAT * r / Math.sqrt(r * r + R_CORE * R_CORE);
}

let W, H;
let stars;            // Float32Array: [r, theta, ring, brightness, size] per star
let vNewtMax, vObsMax;
let vCurveCache;      // precomputed v(r) samples for the plot
const CURVE_SAMPLES = 96;

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

  stars = new Float32Array(N_STARS * 5);

  // Distribute stars: ~30% in the bulge (small r, high density), ~70% in the disk
  // with exponential falloff. Same distribution mirrored across both panels.
  for (let g = 0; g < 2; g++) {
    for (let i = 0; i < N_PER; i++) {
      const idx = g * N_PER + i;
      let r;
      if (i < N_PER * 0.30) {
        // Bulge: r ~ uniform(0, R_BULGE * 1.5) but pulled inward
        r = Math.pow(Math.random(), 1.4) * R_BULGE * 1.8 + 0.02;
      } else {
        // Disk: exponential. Sample u = -ln(1-rand) * Rd, clipped.
        r = -Math.log(Math.max(1e-3, Math.random())) * R_DISK * 0.55 + R_BULGE * 0.6;
        if (r > R_MAX) r = R_MAX * (0.9 + Math.random() * 0.1);
      }
      const theta = Math.random() * Math.PI * 2;
      const ring = Math.floor(r * N_RINGS / R_MAX);
      const bright = 0.55 + Math.random() * 0.45;
      const size = r < R_BULGE ? 1.4 + Math.random() * 0.8 : 0.9 + Math.random() * 0.7;
      stars[idx * 5 + 0] = r;
      stars[idx * 5 + 1] = theta;
      stars[idx * 5 + 2] = ring;
      stars[idx * 5 + 3] = bright;
      stars[idx * 5 + 4] = size;
    }
  }

  // Precompute velocity curves and their maxima for plot normalization.
  vCurveCache = new Float32Array(CURVE_SAMPLES * 2);
  vNewtMax = 0;
  vObsMax = 0;
  for (let i = 0; i < CURVE_SAMPLES; i++) {
    const r = (i + 1) / CURVE_SAMPLES * R_MAX;
    const vn = vNewt(r);
    const vo = vObs(r);
    vCurveCache[i * 2 + 0] = vn;
    vCurveCache[i * 2 + 1] = vo;
    if (vn > vNewtMax) vNewtMax = vn;
    if (vo > vObsMax) vObsMax = vo;
  }

  ctx.fillStyle = '#02030a';
  ctx.fillRect(0, 0, W, H);
}

function drawGalaxy(ctx, cx, cy, panelR, panelW, panelH, glowHue) {
  // Soft disk glow.
  const grad = ctx.createRadialGradient(cx, cy, 0, cx, cy, panelR * 1.15);
  grad.addColorStop(0, `hsla(${glowHue}, 80%, 78%, 0.30)`);
  grad.addColorStop(0.18, `hsla(${glowHue}, 70%, 60%, 0.18)`);
  grad.addColorStop(0.55, `hsla(${glowHue + 20}, 60%, 40%, 0.07)`);
  grad.addColorStop(1, 'hsla(0, 0%, 0%, 0)');
  ctx.fillStyle = grad;
  ctx.beginPath();
  ctx.arc(cx, cy, panelR * 1.15, 0, Math.PI * 2);
  ctx.fill();

  // Central bulge core.
  const core = ctx.createRadialGradient(cx, cy, 0, cx, cy, panelR * 0.22);
  core.addColorStop(0, `hsla(${glowHue + 5}, 100%, 92%, 0.95)`);
  core.addColorStop(0.4, `hsla(${glowHue + 10}, 90%, 75%, 0.55)`);
  core.addColorStop(1, 'hsla(0, 0%, 0%, 0)');
  ctx.fillStyle = core;
  ctx.beginPath();
  ctx.arc(cx, cy, panelR * 0.22, 0, Math.PI * 2);
  ctx.fill();
}

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

  // Slight trail / fade for motion blur on stars.
  ctx.globalCompositeOperation = 'source-over';
  ctx.fillStyle = 'rgba(2, 3, 10, 0.32)';
  ctx.fillRect(0, 0, W, H);

  // Layout: two galaxy panels in the top 60%, v(r) plot in the bottom 40%.
  const topH = H * 0.62;
  const botH = H - topH;
  const panelW = W / 2;
  const cyTop = topH * 0.50;
  const panelR = Math.min(panelW, topH) * 0.40;
  const cxL = panelW * 0.5;
  const cxR = panelW * 1.5;

  // Background star field (a few sparse stars per panel for context).
  ctx.fillStyle = 'rgba(180, 200, 255, 0.45)';
  // Don't redraw bg stars every frame — they'd flicker against the trail fade.
  // Instead skip; the panel glow gives enough atmosphere.

  drawGalaxy(ctx, cxL, cyTop, panelR, panelW, topH, 45);   // warm amber
  drawGalaxy(ctx, cxR, cyTop, panelR, panelW, topH, 200);  // cool blue

  // Advance star angles and draw.
  // We pick an angular-speed scale so that a star at r = 0.5 in the OBSERVED
  // panel completes a full revolution in ~12 seconds. omega = v/r.
  const TIME_SCALE = (2 * Math.PI / 12) * 0.5 / (vObs(0.5) / 0.5);

  ctx.globalCompositeOperation = 'lighter';

  for (let g = 0; g < 2; g++) {
    const cx = g === 0 ? cxL : cxR;
    const useObs = g === 1;
    for (let i = 0; i < N_PER; i++) {
      const idx = (g * N_PER + i) * 5;
      const r = stars[idx + 0];
      let theta = stars[idx + 1];
      const bright = stars[idx + 3];
      const size = stars[idx + 4];

      const v = useObs ? vObs(r) : vNewt(r);
      const omega = (v / Math.max(r, 0.02)) * TIME_SCALE;
      theta += omega * dt;
      if (theta > Math.PI * 2) theta -= Math.PI * 2;
      stars[idx + 1] = theta;

      const px = cx + Math.cos(theta) * r * panelR;
      const py = cyTop + Math.sin(theta) * r * panelR;

      // Bulge stars warmer-white, disk stars bluish-white.
      const hue = r < R_BULGE ? 45 + Math.random() * 10 : 215 + Math.random() * 20;
      const lum = r < R_BULGE ? 85 : 80;
      ctx.fillStyle = `hsla(${hue}, 60%, ${lum}%, ${(bright * 0.85).toFixed(3)})`;
      ctx.beginPath();
      ctx.arc(px, py, size, 0, Math.PI * 2);
      ctx.fill();
    }
  }

  ctx.globalCompositeOperation = 'source-over';

  // --- Labels above the panels ---
  ctx.fillStyle = 'rgba(255, 220, 160, 0.92)';
  ctx.font = 'bold 13px sans-serif';
  ctx.textAlign = 'center';
  ctx.fillText('Newtonian prediction (visible mass only)', cxL, 18);
  ctx.fillStyle = 'rgba(160, 210, 255, 0.92)';
  ctx.fillText('Observed (real galaxies)', cxR, 18);

  ctx.fillStyle = 'rgba(255, 220, 160, 0.65)';
  ctx.font = '11px sans-serif';
  ctx.fillText('outer stars lag — v ∝ 1/√r', cxL, topH - 8);
  ctx.fillStyle = 'rgba(160, 210, 255, 0.65)';
  ctx.fillText('outer stars keep pace — v ≈ const', cxR, topH - 8);

  // --- Bottom panel: v(r) plot ---
  // Clear the plot area cleanly so the curves don't ghost.
  ctx.fillStyle = '#02030a';
  ctx.fillRect(0, topH, W, botH);

  const plotPadL = 56;
  const plotPadR = 24;
  const plotPadT = 16;
  const plotPadB = 28;
  const plotX0 = plotPadL;
  const plotY0 = topH + plotPadT;
  const plotW = W - plotPadL - plotPadR;
  const plotH = botH - plotPadT - plotPadB;

  // Plot frame.
  ctx.strokeStyle = 'rgba(120, 130, 160, 0.55)';
  ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.moveTo(plotX0, plotY0);
  ctx.lineTo(plotX0, plotY0 + plotH);
  ctx.lineTo(plotX0 + plotW, plotY0 + plotH);
  ctx.stroke();

  // Axis labels.
  ctx.fillStyle = 'rgba(190, 200, 220, 0.85)';
  ctx.font = '11px sans-serif';
  ctx.textAlign = 'center';
  ctx.fillText('radius r  (galactocentric distance)', plotX0 + plotW / 2, plotY0 + plotH + 18);
  ctx.save();
  ctx.translate(14, plotY0 + plotH / 2);
  ctx.rotate(-Math.PI / 2);
  ctx.textAlign = 'center';
  ctx.fillText('orbital speed v(r)', 0, 0);
  ctx.restore();

  // Normalize plot to the larger of the two maxima.
  const vMax = Math.max(vNewtMax, vObsMax) * 1.08;

  const xOf = (i) => plotX0 + (i / (CURVE_SAMPLES - 1)) * plotW;
  const yOf = (v) => plotY0 + plotH - (v / vMax) * plotH;

  // Newtonian curve.
  ctx.strokeStyle = 'rgba(255, 200, 110, 0.95)';
  ctx.lineWidth = 2;
  ctx.beginPath();
  for (let i = 0; i < CURVE_SAMPLES; i++) {
    const vn = vCurveCache[i * 2 + 0];
    const x = xOf(i);
    const y = yOf(vn);
    if (i === 0) ctx.moveTo(x, y);
    else ctx.lineTo(x, y);
  }
  ctx.stroke();

  // Observed curve.
  ctx.strokeStyle = 'rgba(140, 200, 255, 0.95)';
  ctx.lineWidth = 2;
  ctx.beginPath();
  for (let i = 0; i < CURVE_SAMPLES; i++) {
    const vo = vCurveCache[i * 2 + 1];
    const x = xOf(i);
    const y = yOf(vo);
    if (i === 0) ctx.moveTo(x, y);
    else ctx.lineTo(x, y);
  }
  ctx.stroke();

  // Shaded gap = inferred dark-matter contribution.
  ctx.fillStyle = 'rgba(180, 140, 255, 0.18)';
  ctx.beginPath();
  // Top edge = observed curve
  for (let i = 0; i < CURVE_SAMPLES; i++) {
    const vo = vCurveCache[i * 2 + 1];
    const x = xOf(i);
    const y = yOf(vo);
    if (i === 0) ctx.moveTo(x, y);
    else ctx.lineTo(x, y);
  }
  // Back along Newtonian curve.
  for (let i = CURVE_SAMPLES - 1; i >= 0; i--) {
    const vn = vCurveCache[i * 2 + 0];
    const x = xOf(i);
    const y = yOf(vn);
    ctx.lineTo(x, y);
  }
  ctx.closePath();
  ctx.fill();

  // Legend.
  const legX = plotX0 + plotW - 200;
  const legY = plotY0 + 6;
  ctx.font = '11px sans-serif';
  ctx.textAlign = 'left';

  ctx.strokeStyle = 'rgba(255, 200, 110, 0.95)';
  ctx.lineWidth = 2;
  ctx.beginPath();
  ctx.moveTo(legX, legY + 6);
  ctx.lineTo(legX + 22, legY + 6);
  ctx.stroke();
  ctx.fillStyle = 'rgba(255, 220, 160, 0.95)';
  ctx.fillText('Newtonian (visible mass)', legX + 28, legY + 9);

  ctx.strokeStyle = 'rgba(140, 200, 255, 0.95)';
  ctx.beginPath();
  ctx.moveTo(legX, legY + 22);
  ctx.lineTo(legX + 22, legY + 22);
  ctx.stroke();
  ctx.fillStyle = 'rgba(160, 210, 255, 0.95)';
  ctx.fillText('Observed (flat rotation)', legX + 28, legY + 25);

  ctx.fillStyle = 'rgba(200, 170, 255, 0.95)';
  ctx.fillRect(legX, legY + 34, 22, 8);
  ctx.fillStyle = 'rgba(210, 180, 255, 0.95)';
  ctx.fillText('inferred dark matter', legX + 28, legY + 41);
}