guts/web/js/world/spline.js

// world/spline.js (Lane A) — the math core of the canal. Deliberately imports NOTHING but
// core/rng.js: no three.js, so `node web/js/world/spline.js --selfcheck` runs it headless in
// qa.sh. index.js wraps these plain {x,y,z} numbers into THREE objects at the contract
// boundary (docs/TECH.md §THE WORLD CONTRACT).
//
// Model
// -----
// The centreline is a graph over the -Z axis: P(u) = (X(u), Y(u), -u). Two consequences we
// rely on everywhere: it can never self-intersect, and `u ≈ -pos.z` is a free, close seed for
// project(). X/Y are sums of three seeded low-frequency sines, smoothstep-blended across
// segment joins so the tangent never kinks.
//
// Amplitudes are NOT authored in units — they're solved from a curvature budget. A tube bends
// through itself when the centreline's turn radius drops below the tube radius, so `curviness`
// means "fraction of the tightest bend this pipe's radius can survive":
//     kappa_budget(s) = curviness(s) / (PINCH_SAFETY * baseRadius(s))
// split across axes and octaves by fixed weights, then inverted per octave (A = w*kappa/f^2).
// A wide pipe therefore straightens itself automatically and Lane C never has to think about
// curvature. selfcheck measures the result (`pinchRatio`) rather than trusting this argument —
// it caught the first version of this file, where amplitude was authored in units and a
// curviness-0.95 join into a radius-16 segment folded the tube shut.
//
// u is NOT arclength (|dP/du| >= 1). Everything public is in arclength `s`, mapped through
// LUTs built by one forward march.

import { createRng } from '../core/rng.js';

// --- tiny vec3 (plain objects; index.js converts) ---------------------------------------
const V = (x = 0, y = 0, z = 0) => ({ x, y, z });
const sub = (a, b) => V(a.x - b.x, a.y - b.y, a.z - b.z);
const dot = (a, b) => a.x * b.x + a.y * b.y + a.z * b.z;
const len = (a) => Math.sqrt(dot(a, a));
const norm = (a) => { const l = len(a) || 1; return V(a.x / l, a.y / l, a.z / l); };
const cross = (a, b) => V(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
const lerpV = (a, b, t) => V(a.x + (b.x - a.x) * t, a.y + (b.y - a.y) * t, a.z + (b.z - a.z) * t);
const clamp = (x, lo, hi) => (x < lo ? lo : x > hi ? hi : x);
const lerp = (a, b, t) => a + (b - a) * t;
const smoothstep = (e0, e1, x) => { const t = clamp((x - e0) / (e1 - e0), 0, 1); return t * t * (3 - 2 * t); };

// --- tuning -----------------------------------------------------------------------------
// OMEGA is the peristaltic contraction *rate*, global to the whole canal and constant in
// time — physiologically it's how often the muscle fires, which doesn't change because you
// crossed into a wider pipe. The wavenumber is what varies: k(s) = OMEGA / flow(s), so a
// wave crest travels at exactly the local flow speed and Lane B can surf one. Phase is the
// integral K(s) = ∫k ds (phaseAt), which stays continuous across flow changes — plain
// `k*s - w*t` does not, and tears the wave at every segment join.
// 3.08 = 0.22 rad/unit x 14 units/s, i.e. the stub's esophagus wave, preserved exactly.
export const OMEGA = 3.08;

const DU = 0.5;        // march step in curve parameter
const DS = 0.5;        // frame/phase LUT spacing in arclength
const NOISE_TAB = 1024;
const RADIUS_WAVELENGTH = 55;   // units per radius-fbm octave-0 cycle
const PINCH_SAFETY = 2.2;       // min turn radius = this x base radius (>2 covers wobble peaks)

// Blend widths (units) for smoothstep-crossfading segment params across a join. Radius and
// flow are local scalars — a tight 12u transition reads as a sphincter. Curviness feeds
// centreline *amplitude*, which has a long lever arm: ramping it over a short span is itself a
// hard lateral swerve, i.e. curvature. It gets a wide, gentle ramp.
// `curvBase` is baseRadius again, but read through the wide ramp: it divides the curvature
// budget, so a 12u step in it would swerve the centreline just as hard as curviness would.
const BLEND = { curviness: 60, curvBase: 60, base: 12, wobble: 12, flow: 12 };

// [wavelength, share of the axis curvature budget]. Y is tamer than X on purpose: a canal that
// writhes vertically as hard as it does laterally would swing the parallel-transport frame
// around and make Lane B's chase cam seasick. Axis shares are set so that even with every
// octave peaking together, sqrt(X^2 + Y^2) lands exactly on the budget.
const OCT_X = [[240, 0.5], [110, 0.3], [60, 0.2]];
const OCT_Y = [[190, 0.5], [85, 0.3], [45, 0.2]];
const AXIS_X = 0.8, AXIS_Y = 0.6;

/** Solve each octave's amplitude from its share of a unit curvature budget: kappa ~ A*f^2. */
function octaveCoefs(octaves, axisShare) {
  return octaves.map(([wavelength, weight]) => {
    const f = (2 * Math.PI) / wavelength;
    return { f, coef: (weight * axisShare) / (f * f) };
  });
}

const GET = {
  curviness: (g) => (typeof g.curviness === 'number' ? g.curviness : 0.5),
  base: (g) => (g.radius && typeof g.radius.base === 'number' ? g.radius.base : 10),
  wobble: (g) => (g.radius && typeof g.radius.wobble === 'number' ? g.radius.wobble : 0.2),
  flow: (g) => (typeof g.flow === 'number' ? g.flow : 10),
};
GET.curvBase = GET.base;

/**
 * Build the canal math for a level. Deterministic: same `level.seed` (or injected rng) =>
 * byte-identical output, forever. See world.hash().
 * @param {object} level  level JSON (TECH.md §Level data schema v0)
 * @param {function} [rng]  core/rng.js createRng() instance; defaults to level.seed
 */
export function buildSpline(level, rng) {
  if (!level || !Array.isArray(level.segments) || level.segments.length === 0)
    throw new Error('[world] level.segments is required and must be non-empty');
  const R = rng || createRng((level.seed ?? 0) >>> 0);

  // --- segment table (declared arclength boundaries; C's events index this same s) -------
  const spans = [];
  let acc = 0;
  for (const seg of level.segments) {
    const length = typeof seg.length === 'number' && seg.length > 0 ? seg.length : 100;
    spans.push({ s0: acc, s1: acc + length, seg });
    acc += length;
  }
  const L = acc;

  function idxAt(s) {
    // linear scan: levels have a handful of segments, and a branchy binary search here would
    // cost more than it saves.
    for (let i = 0; i < spans.length; i++) if (s < spans[i].s1) return i;
    return spans.length - 1;
  }

  /** Segment param at s, smoothstep-blended across joins so nothing steps discontinuously. */
  function paramAt(s, key) {
    const get = GET[key], w = BLEND[key];
    const i = idxAt(s);
    const sp = spans[i];
    const v = get(sp.seg);
    if (i > 0 && s < sp.s0 + w)
      return lerp(get(spans[i - 1].seg), v, smoothstep(sp.s0 - w, sp.s0 + w, s));
    if (i < spans.length - 1 && s > sp.s1 - w)
      return lerp(v, get(spans[i + 1].seg), smoothstep(sp.s1 - w, sp.s1 + w, s));
    return v;
  }

  // --- seeded noise ---------------------------------------------------------------------
  const tab = new Float32Array(NOISE_TAB);
  { const r = R('world.radius'); for (let i = 0; i < NOISE_TAB; i++) tab[i] = r() * 2 - 1; }
  function vnoise(x) {
    const i = Math.floor(x), f = x - i;
    const a = tab[i & (NOISE_TAB - 1)], b = tab[(i + 1) & (NOISE_TAB - 1)];
    return lerp(a, b, f * f * (3 - 2 * f));
  }
  function fbm1(x) {
    let sum = 0, amp = 1, freq = 1, normSum = 0;
    for (let o = 0; o < 3; o++) { sum += amp * vnoise(x * freq + o * 17.3); normSum += amp; amp *= 0.5; freq *= 2.03; }
    return sum / normSum;              // [-1, 1]
  }
  const radiusAt = (s) => paramAt(s, 'base') * (1 + paramAt(s, 'wobble') * fbm1(s / RADIUS_WAVELENGTH));

  // --- centreline -----------------------------------------------------------------------
  const cx = octaveCoefs(OCT_X, AXIS_X), cy = octaveCoefs(OCT_Y, AXIS_Y);
  const phX = cx.map(() => R('world.spline')() * Math.PI * 2);
  const phY = cy.map(() => R('world.spline')() * Math.PI * 2);

  /** How hard the canal is allowed to bend at s: curviness as a fraction of the pinch limit. */
  const curvBudgetAt = (s) => paramAt(s, 'curviness') / (PINCH_SAFETY * Math.max(1, paramAt(s, 'curvBase')));

  /** `budget` is kappa (1/units), not a length — octave coefs turn it into lateral offsets. */
  function centreWith(u, budget) {
    let x = 0, y = 0;
    for (let i = 0; i < cx.length; i++) x += budget * cx[i].coef * Math.sin(u * cx[i].f + phX[i]);
    for (let i = 0; i < cy.length; i++) y += budget * cy[i].coef * Math.sin(u * cy[i].f + phY[i]);
    return V(x, y, -u);
  }

  // --- the march: u -> s, and the curvature budget sampled along it ----------------------
  // The budget is a function of arclength (that's where segments are declared), but the curve
  // is evaluated in u — so we resolve it during the march, where s is already known, and
  // cache the schedule per u-step. One pass, no circularity.
  const uToS = [0];
  const ampArr = [curvBudgetAt(0)];
  {
    let s = 0, prev = centreWith(0, ampArr[0]), n = 0;
    const MAXN = Math.ceil((L / DU) * 1.5) + 16;
    while (s < L && n < MAXN) {
      n++;
      const amp = curvBudgetAt(s);             // s lags one step: deterministic, sub-step error
      const p = centreWith(n * DU, amp);
      s += len(sub(p, prev));
      uToS.push(s); ampArr.push(amp);
      prev = p;
    }
  }
  const uMax = (uToS.length - 1) * DU;
  const ampAtU = (u) => {
    const t = clamp(u / DU, 0, ampArr.length - 1), i = Math.min(Math.floor(t), ampArr.length - 2);
    return lerp(ampArr[i], ampArr[i + 1], t - i);
  };
  const centre = (u) => centreWith(u, ampAtU(u));
  const sOfU = (u) => {
    const t = clamp(u / DU, 0, uToS.length - 1), i = Math.min(Math.floor(t), uToS.length - 2);
    return lerp(uToS[i], uToS[i + 1], t - i);
  };

  // uniform s -> u table, inverted by a monotone merge-walk (uToS is strictly increasing)
  const M = Math.ceil(L / DS) + 1;
  const sToU = new Float64Array(M);
  {
    let j = 0;
    for (let i = 0; i < M; i++) {
      const target = Math.min(i * DS, L);
      while (j < uToS.length - 2 && uToS[j + 1] < target) j++;
      const s0 = uToS[j], s1 = uToS[j + 1];
      const f = s1 > s0 ? (target - s0) / (s1 - s0) : 0;
      sToU[i] = (j + f) * DU;
    }
  }
  const uOfS = (s) => {
    const t = clamp(s, 0, L) / DS, i = Math.min(Math.floor(t), M - 2);
    return lerp(sToU[i], sToU[i + 1], t - i);
  };

  // --- frames: parallel transport, NOT Frenet -------------------------------------------
  // Frenet's normal flips through inflection points; on a canal that writhes that shows up as
  // the camera snapping upside-down mid-corridor (LANE_A_WORLD.md calls this out by name).
  // Parallel transport carries the previous normal forward by the minimal rotation, so it's
  // history-dependent => must be integrated once here and looked up, never recomputed ad hoc.
  const posArr = new Float64Array(M * 3);
  const tanArr = new Float64Array(M * 3);
  const norArr = new Float64Array(M * 3);
  {
    const H = 0.25;
    const tangentAtU = (u) => {
      const a = centre(Math.max(0, u - H)), b = centre(Math.min(uMax, u + H));
      return norm(sub(b, a));
    };
    const put = (arr, i, v) => { arr[i * 3] = v.x; arr[i * 3 + 1] = v.y; arr[i * 3 + 2] = v.z; };

    let prevT = tangentAtU(uOfS(0));
    // seed the frame from world up; the canal is a graph over -Z so it can never be vertical
    let n = norm(sub(V(0, 1, 0), V(prevT.x * prevT.y, prevT.y * prevT.y, prevT.z * prevT.y)));
    put(posArr, 0, centre(uOfS(0))); put(tanArr, 0, prevT); put(norArr, 0, n);

    for (let i = 1; i < M; i++) {
      const u = uOfS(Math.min(i * DS, L));
      const p = centre(u), t = tangentAtU(u);
      const axis = cross(prevT, t);
      const sinA = len(axis), cosA = dot(prevT, t);
      if (sinA > 1e-9) {                       // Rodrigues: rotate n by the frame's own turn
        const k = norm(axis), angle = Math.atan2(sinA, cosA);
        const c = Math.cos(angle), sN = Math.sin(angle), kv = cross(k, n), kd = dot(k, n);
        n = V(n.x * c + kv.x * sN + k.x * kd * (1 - c),
              n.y * c + kv.y * sN + k.y * kd * (1 - c),
              n.z * c + kv.z * sN + k.z * kd * (1 - c));
      }
      n = norm(V(n.x - t.x * dot(n, t), n.y - t.y * dot(n, t), n.z - t.z * dot(n, t)));  // re-orthogonalize
      put(posArr, i, p); put(tanArr, i, t); put(norArr, i, n);
      prevT = t;
    }
  }
  const at = (arr, i) => V(arr[i * 3], arr[i * 3 + 1], arr[i * 3 + 2]);

  // --- peristalsis phase: K(s) = ∫ k dx, trapezoid on the same grid ----------------------
  const kAt = (s) => OMEGA / Math.max(1, paramAt(s, 'flow'));
  const phaseArr = new Float64Array(M);
  for (let i = 1; i < M; i++)
    phaseArr[i] = phaseArr[i - 1] + 0.5 * (kAt((i - 1) * DS) + kAt(i * DS)) * DS;
  const phaseAt = (s) => {
    const t = clamp(s, 0, L) / DS, i = Math.min(Math.floor(t), M - 2);
    return lerp(phaseArr[i], phaseArr[i + 1], t - i);
  };

  /** Frame at arclength s. pos/tan/nor/bin + geometric radius. */
  function frameAt(s) {
    const t = clamp(s, 0, L) / DS, i = Math.min(Math.floor(t), M - 2), f = t - i;
    const pos = lerpV(at(posArr, i), at(posArr, i + 1), f);
    const tan = norm(lerpV(at(tanArr, i), at(tanArr, i + 1), f));
    let nor = lerpV(at(norArr, i), at(norArr, i + 1), f);
    nor = norm(V(nor.x - tan.x * dot(nor, tan), nor.y - tan.y * dot(nor, tan), nor.z - tan.z * dot(nor, tan)));
    return { pos, tan, nor, bin: cross(tan, nor), radius: radiusAt(s) };
  }

  /**
   * World -> spline space. Seeded from u = -pos.z (exact for the centreline's own axis, so
   * the seed error is only the lateral tilt), then fixed-point: s += (p - C(s))·T(s).
   * Converges at rate ~ curvature*rho (<0.5 in tube mode). Step-clamped so a point far off
   * the centreline (arena mode) walks instead of exploding.
   */
  function project(p, hint) {
    let s = typeof hint === 'number' ? clamp(hint, 0, L) : sOfU(clamp(-p.z, 0, uMax));
    for (let i = 0; i < 6; i++) {
      const f = frameAt(s);
      const d = dot(sub(p, f.pos), f.tan);
      const next = clamp(s + clamp(d, -8, 8), 0, L);
      if (Math.abs(next - s) < 1e-4) { s = next; break; }
      s = next;
    }
    const f = frameAt(s);
    const rel = sub(p, f.pos);
    const along = dot(rel, f.tan);
    return {
      s,
      theta: Math.atan2(dot(rel, f.bin), dot(rel, f.nor)),
      rho: Math.sqrt(Math.max(0, dot(rel, rel) - along * along)),   // true cross-section radius
    };
  }

  /** Golden determinism hash over sampled frames + radius + wave phase. qa.sh compares. */
  function hash() {
    let h = 2166136261 >>> 0;
    const push = (x) => { h ^= Math.round(x * 1000) & 0xffff; h = Math.imul(h, 16777619); };
    for (let s = 0; s <= L; s += 10) {
      const f = frameAt(s);
      push(f.pos.x); push(f.pos.y); push(f.pos.z); push(f.radius);
      push(f.nor.x); push(f.nor.y); push(f.nor.z); push(phaseAt(s));
    }
    return (h >>> 0).toString(16);
  }

  /** Measured geometry health — the pinch guard is the one that bites. */
  function stats() {
    let maxCurvature = 0, sAtMax = 0, maxRadius = 0, minRadius = Infinity;
    for (let i = 1; i < M - 1; i++) {
      const k = len(sub(at(tanArr, i + 1), at(tanArr, i - 1))) / (2 * DS);
      if (k > maxCurvature) { maxCurvature = k; sAtMax = i * DS; }
    }
    for (let s = 0; s <= L; s += 1) {
      const r = radiusAt(s);
      if (r > maxRadius) maxRadius = r;
      if (r < minRadius) minRadius = r;
    }
    const minTurnRadius = maxCurvature > 0 ? 1 / maxCurvature : Infinity;
    return {
      length: L, gridPoints: M, uMax,
      maxCurvature, sAtMaxCurvature: sAtMax, minTurnRadius,
      minRadius, maxRadius,
      pinchRatio: minTurnRadius / maxRadius,   // <1 => tube folds through itself on a bend
    };
  }

  return {
    length: L, spans, uMax,
    centre, sOfU, uOfS, ampAtU,
    radiusAt, frameAt, project, phaseAt, kAt, paramAt,
    omega: OMEGA,
    hash, stats,
    segmentAt: (s) => spans[idxAt(clamp(s, 0, L))].seg,
  };
}

// ---------------------------------------------------------------------------------------
// headless selfcheck: `node web/js/world/spline.js --selfcheck`   (qa.sh gate)
// Browsers never reach this (no `process`).
// ---------------------------------------------------------------------------------------
async function selfcheck() {
  const { FIXTURE } = await import('./_fixture.js');
  let fail = 0;
  const ok = (name, cond, detail = '') => {
    if (cond) console.log(`  \x1b[32m✓\x1b[0m ${name}${detail ? ' — ' + detail : ''}`);
    else { fail++; console.log(`  \x1b[31m✗\x1b[0m ${name}${detail ? ' — ' + detail : ''}`); }
  };
  console.log('world/spline selfcheck (fixture: %s)', FIXTURE.id);

  const a = buildSpline(FIXTURE);
  const b = buildSpline(FIXTURE);
  const c = buildSpline({ ...FIXTURE, seed: FIXTURE.seed + 1 });
  const st = a.stats();

  ok('determinism: same seed => same hash', a.hash() === b.hash(), a.hash());
  ok('determinism: different seed => different hash', a.hash() !== c.hash(), `${a.hash()} vs ${c.hash()}`);

  const declared = FIXTURE.segments.reduce((n, s) => n + s.length, 0);
  ok('arclength: length == declared total', Math.abs(a.length - declared) < 1e-6, `${a.length}`);

  let mono = true, prev = -1;
  for (let s = 0; s <= a.length; s += 1) { const u = a.uOfS(s); if (u < prev) mono = false; prev = u; }
  ok('arclength: s -> u monotone', mono);

  let maxRt = 0;
  for (let s = 0; s <= a.length; s += 3) maxRt = Math.max(maxRt, Math.abs(a.sOfU(a.uOfS(s)) - s));
  ok('arclength: s -> u -> s round-trip < 0.05', maxRt < 0.05, `max err ${maxRt.toExponential(2)}`);

  let orth = 0, unit = 0;
  for (let s = 0; s <= a.length; s += 2) {
    const f = a.frameAt(s);
    orth = Math.max(orth, Math.abs(f.tan.x * f.nor.x + f.tan.y * f.nor.y + f.tan.z * f.nor.z));
    unit = Math.max(unit, Math.abs(Math.hypot(f.nor.x, f.nor.y, f.nor.z) - 1), Math.abs(Math.hypot(f.bin.x, f.bin.y, f.bin.z) - 1));
  }
  ok('frames: orthonormal (tan·nor ~ 0, |nor| = |bin| = 1)', orth < 1e-9 && unit < 1e-9,
     `tan·nor ${orth.toExponential(1)}, unit err ${unit.toExponential(1)}`);

  // parallel transport must not spin the frame up: total twist stays bounded & smooth
  let maxTwist = 0;
  for (let s = 1; s <= a.length; s += 1) {
    const p = a.frameAt(s - 1).nor, q = a.frameAt(s).nor;
    maxTwist = Math.max(maxTwist, Math.acos(Math.min(1, Math.abs(p.x * q.x + p.y * q.y + p.z * q.z))));
  }
  ok('frames: no roll snap (max per-unit normal turn < 0.2 rad)', maxTwist < 0.2, `${maxTwist.toFixed(4)} rad`);

  ok('geometry: tube does not pinch (minTurnRadius / maxRadius > 1.5)', st.pinchRatio > 1.5,
     `turn R ${st.minTurnRadius.toFixed(1)} vs tube R ${st.maxRadius.toFixed(1)} => ${st.pinchRatio.toFixed(2)}x`);

  // project() round-trip: place a point at a known (s, theta, rho), recover it
  let maxSerr = 0, maxTerr = 0, maxRerr = 0;
  for (let s = 5; s < a.length - 5; s += 7) {
    for (const theta of [0, 1.3, -2.4, 3.0]) {
      const rho = a.radiusAt(s) * 0.6;
      const f = a.frameAt(s);
      const p = {
        x: f.pos.x + (f.nor.x * Math.cos(theta) + f.bin.x * Math.sin(theta)) * rho,
        y: f.pos.y + (f.nor.y * Math.cos(theta) + f.bin.y * Math.sin(theta)) * rho,
        z: f.pos.z + (f.nor.z * Math.cos(theta) + f.bin.z * Math.sin(theta)) * rho,
      };
      const q = a.project(p);
      maxSerr = Math.max(maxSerr, Math.abs(q.s - s));
      let dt = Math.abs(q.theta - theta) % (Math.PI * 2);
      if (dt > Math.PI) dt = Math.PI * 2 - dt;
      maxTerr = Math.max(maxTerr, dt);
      maxRerr = Math.max(maxRerr, Math.abs(q.rho - rho));
    }
  }
  ok('project: recovers s within 0.25', maxSerr < 0.25, `max ${maxSerr.toFixed(4)}`);
  ok('project: recovers theta within 0.02 rad', maxTerr < 0.02, `max ${maxTerr.toFixed(5)}`);
  ok('project: recovers rho within 0.05', maxRerr < 0.05, `max ${maxRerr.toFixed(4)}`);

  // wave: phase monotone, and crest speed == local flow (the whole point of K(s))
  let phaseMono = true;
  for (let s = 1; s <= a.length; s += 1) if (a.phaseAt(s) <= a.phaseAt(s - 1)) phaseMono = false;
  ok('wave: phase K(s) strictly increasing', phaseMono);
  const kEso = a.kAt(50), kSto = a.kAt(a.length - 20);
  ok('wave: k == OMEGA/flow (esophagus flow 14 => k 0.22, matches stub)', Math.abs(kEso - 0.22) < 0.001, `k=${kEso.toFixed(4)}`);
  ok('wave: crest speed == local flow in each biome',
     Math.abs(OMEGA / kEso - 14) < 0.01 && Math.abs(OMEGA / kSto - 4) < 0.01,
     `${(OMEGA / kEso).toFixed(2)} u/s eso, ${(OMEGA / kSto).toFixed(2)} u/s stomach`);

  console.log(`  stats: L=${st.length} grid=${st.gridPoints} maxK=${st.maxCurvature.toFixed(4)} (s=${st.sAtMaxCurvature}) radius=${st.minRadius.toFixed(1)}..${st.maxRadius.toFixed(1)} hash=${a.hash()}`);
  console.log(fail === 0 ? '\x1b[32mworld/spline: OK\x1b[0m' : `\x1b[31mworld/spline: ${fail} FAILED\x1b[0m`);
  return fail;
}

if (typeof process !== 'undefined' && process.argv && process.argv.includes('--selfcheck')) {
  selfcheck().then((f) => process.exit(f === 0 ? 0 : 1));
}