HardYards/web/world/js/sail.js

/**
 * sail.js — shade sail cloth simulation, corner loads, hardware failure.  [Lane B]
 *
 * A 3D verlet cloth on a bilinear patch between 4 anchors. Wind pressure is
 * applied per FACE, not per node, which is the whole point: a twisted (hypar)
 * sail turns most of its faces edge-on to the wind and sheds load, while a flat
 * one presents every face square-on and catches everything. That difference is
 * the game's thesis and it is asserted in sail.selftest.js.
 *
 * Units are SI throughout: metres, kilograms, seconds, newtons. Corner loads
 * come out in real newtons and hardware ratings are real working load limits,
 * so a 5x5 m sail in a 30 m/s gust genuinely puts ~5 kN on a corner — which is
 * genuinely why real shade sails use 3 kN+ shackles.
 *
 * This module deliberately does NOT import three.js. The sim core is plain
 * typed arrays so it runs headless under node (see sail.selftest.js) before any
 * renderer exists, stays allocation-free in the hot loop, and can be replayed
 * bit-for-bit. three.js is pulled in lazily by createSailView() only.
 */

// ---------- sim tunables ----------
const SIM_DT = 1 / 60;        // sim always steps at a fixed rate; step() accumulates
const MAX_SUBSTEPS = 5;       // spiral-of-death guard when the frame hitches
const RELAX_ITERS = 5;        // FABRIC_K is calibrated against this; changing it rescales loads
const GRAVITY = -9.81;

// ---------- aerodynamics ----------
// 0.5 * air density (1.225) * flat-plate drag coefficient (~1.4).
// Newtons per m^2 of face area per (m/s)^2 of normal-on airflow.
const PRESSURE_COEFF = 0.86;
const TANGENT_COEFF = 0.02;   // skin friction dragging along the face
const MAX_NORMAL_SPEED = 45;  // clamp on the normal-on component, m/s — stability in extreme gusts

// ---------- fabric ----------
const FABRIC_DENSITY = 0.32;  // kg/m^2, typical knitted shade cloth
// Axial stiffness of one grid spring, N/m — roughly E*t*width/length for
// knitted HDPE mesh. Fed to the solver as a compliance (1/k), not used to
// convert stretch into force: see _measureLoads for why that distinction is
// the whole ballgame.
const FABRIC_K = 100000;
const K_COMPRESS = 0.08;      // cloth resists stretch hard, compression barely (from prototype)
const K_BEND = 0.04;
const COMP_STRETCH = 1 / FABRIC_K;
const COMP_COMPRESS = 1 / (FABRIC_K * K_COMPRESS);
const COMP_BEND = 1 / (FABRIC_K * K_BEND);
const VEL_DAMP = 0.995;       // light; relative-wind drag supplies the real damping

// ---------- failure ----------
const OVERLOAD_SECS = 0.4;    // prototype: 0.4 s sustained overload before it lets go
const OVERLOAD_RECOVER = 2.0; // prototype: overload timer bleeds off at 2x
const LOAD_TAU = 0.11;        // load meter smoothing time constant, s

/**
 * Hardware tiers. Costs are the prototype's economy verbatim; ratings are
 * retuned from the prototype's arbitrary 9/19/40 into real newtons, preserving
 * the same relative spread. `rating` is a working load limit in N.
 */
export const HARDWARE = [
  { name: 'carabiner', cost: 5, rating: 1200, color: '#e2b04a' },
  { name: 'shackle', cost: 15, rating: 3200, color: '#c8d2d8' },
  { name: 'rated shackle', cost: 30, rating: 6500, color: '#7ee0ff' },
];

export const TENSION_MIN = 0.6;
export const TENSION_MAX = 1.4;
const TRIM_MIN = 0.85;
const TRIM_MAX = 1.15;

const clamp = (v, lo, hi) => (v < lo ? lo : v > hi ? hi : v);

/**
 * Order 4 anchors into a non-self-intersecting ring by angle around their
 * centroid, projected onto the ground plane. Ported from the prototype's
 * orderRing; without it, picking corners in a silly order knots the sail.
 */
export function orderRing(anchors) {
  const n = anchors.length;
  let cx = 0, cz = 0;
  for (const a of anchors) { cx += a.pos.x; cz += a.pos.z; }
  cx /= n; cz /= n;
  return [...anchors].sort(
    (a, b) => Math.atan2(a.pos.z - cz, a.pos.x - cx) - Math.atan2(b.pos.z - cz, b.pos.x - cx)
  );
}

export class SailRig {
  /**
   * @param {object}   opts
   * @param {Array}    opts.anchors  world.anchors — [{id, pos:{x,y,z}, type, sway(t)->{x,y,z}}]
   *                                 sway(t) returns an OFFSET to add to pos, not an absolute position.
   * @param {number}   opts.gridN    nodes per side (default 10)
   * @param {number}   opts.porosity 0 = solid membrane, ~0.3 = knitted shade cloth (blows through, less load)
   */
  constructor({ anchors = [], gridN = 10, porosity = 0 } = {}) {
    this.anchors = anchors;
    this.N = gridN;
    this.porosity = porosity;
    this.corners = [];
    this.events = [];
    this.tension = 1.0;
    this.t = 0;
    this.rigged = false;
    this._acc = 0;
    // scratch, reused every face to keep the hot loop allocation-free
    this._probe = { x: 0, y: 0, z: 0 };
  }

  /**
   * Rig the sail across 4 anchors.
   * @param {string[]} anchorIds  4 anchor ids; reordered into a ring internally
   * @param {object[]} hwChoices  hardware per anchor id, same order as anchorIds
   * @param {number}   tension    0.6 (loose, flogs) .. 1.4 (drum tight, shock-loads)
   */
  attach(anchorIds, hwChoices, tension = 1.0) {
    if (anchorIds.length !== 4) throw new Error(`sail needs exactly 4 corners, got ${anchorIds.length}`);

    const picked = anchorIds.map((id) => {
      const a = this.anchors.find((x) => x.id === id);
      if (!a) throw new Error(`unknown anchor "${id}"`);
      return a;
    });
    const hwById = new Map(anchorIds.map((id, i) => [id, hwChoices[i] || HARDWARE[0]]));

    const ring = orderRing(picked);
    this.tension = clamp(tension, TENSION_MIN, TENSION_MAX);
    this.corners = ring.map((a) => ({
      anchorId: a.id,
      anchor: a,
      hw: hwById.get(a.id),
      load: 0,
      peakLoad: 0,
      overload: 0,
      broken: false,
      trim: 1.0,
      loadVec: { x: 0, y: 0, z: 0 }, // reaction direction, not just magnitude — see _measureLoads
    }));

    this._build(ring);
    this.rigged = true;
    return this;
  }

  _build(ring) {
    const N = this.N;
    const nodeCount = N * N;
    this.pos = new Float64Array(nodeCount * 3);
    this.prev = new Float64Array(nodeCount * 3);
    this.force = new Float64Array(nodeCount * 3);
    this.invMass = new Float64Array(nodeCount);

    // Bilinear patch across the 4 corners. Because the anchors sit at different
    // heights, this initial surface is already a hypar — the sim just relaxes it.
    const [c0, c1, c2, c3] = ring.map((a) => a.pos);
    for (let v = 0; v < N; v++) {
      for (let u = 0; u < N; u++) {
        const fu = u / (N - 1), fv = v / (N - 1);
        const i = (v * N + u) * 3;
        for (let k = 0; k < 3; k++) {
          const ax = ['x', 'y', 'z'][k];
          const top = (1 - fu) * c0[ax] + fu * c1[ax];
          const bot = (1 - fu) * c3[ax] + fu * c2[ax];
          this.pos[i + k] = this.prev[i + k] = (1 - fv) * top + fv * bot;
        }
      }
    }

    const idx = (u, v) => v * N + u;
    this.cornerIdx = [idx(0, 0), idx(N - 1, 0), idx(N - 1, N - 1), idx(0, N - 1)];

    // springs: structural + shear carry load; bend only resists folding
    this.springs = [];
    const link = (a, b, kind) => {
      const ax = a * 3, bx = b * 3;
      const dx = this.pos[bx] - this.pos[ax];
      const dy = this.pos[bx + 1] - this.pos[ax + 1];
      const dz = this.pos[bx + 2] - this.pos[ax + 2];
      this.springs.push({ a, b, restBase: Math.hypot(dx, dy, dz), rest: 0, kind });
    };
    for (let v = 0; v < N; v++) {
      for (let u = 0; u < N; u++) {
        if (u < N - 1) link(idx(u, v), idx(u + 1, v), 'struct');
        if (v < N - 1) link(idx(u, v), idx(u, v + 1), 'struct');
        if (u < N - 1 && v < N - 1) {
          link(idx(u, v), idx(u + 1, v + 1), 'shear');
          link(idx(u + 1, v), idx(u, v + 1), 'shear');
        }
        if (u < N - 2) link(idx(u, v), idx(u + 2, v), 'bend');
        if (v < N - 2) link(idx(u, v), idx(u, v + 2), 'bend');
      }
    }

    // XPBD Lagrange multipliers, one per spring, reset every substep
    this.lambda = new Float64Array(this.springs.length);

    // Springs meeting each corner, kept as {spring, index} so the load meter can
    // look up each one's multiplier. Bend springs are included: the hardware
    // physically carries every element that touches it, and leaving them out
    // under-reports the reaction and breaks the statics balance.
    this.cornerSprings = this.cornerIdx.map((ci) =>
      this.springs
        .map((s, si) => ({ s, si }))
        .filter(({ s }) => s.a === ci || s.b === ci)
    );

    // triangles: wind acts per face, and coverage raycasts against these
    this.tris = new Uint16Array((N - 1) * (N - 1) * 6);
    let ti = 0;
    for (let v = 0; v < N - 1; v++) {
      for (let u = 0; u < N - 1; u++) {
        const a = idx(u, v), b = idx(u + 1, v), c = idx(u + 1, v + 1), d = idx(u, v + 1);
        this.tris[ti++] = a; this.tris[ti++] = b; this.tris[ti++] = c;
        this.tris[ti++] = a; this.tris[ti++] = c; this.tris[ti++] = d;
      }
    }

    // grid-space proximity of every node to each corner, for per-corner trim
    this._cornerWeight = [];
    for (let k = 0; k < 4; k++) {
      const cu = [0, N - 1, N - 1, 0][k], cv = [0, 0, N - 1, N - 1][k];
      const w = new Float64Array(nodeCount);
      for (let v = 0; v < N; v++) {
        for (let u = 0; u < N; u++) {
          const dist = Math.hypot(u - cu, v - cv) / (N - 1);
          w[idx(u, v)] = Math.max(0, 1 - dist);
        }
      }
      this._cornerWeight.push(w);
    }

    this.area = this._surfaceArea();
    const mass = (FABRIC_DENSITY * this.area) / nodeCount;
    this.nodeMass = mass;
    this.invMass.fill(1 / mass);

    this._applyRestLengths();
    this._repin(0);
  }

  /** Rest lengths shrink as tension rises (1/tension, ported), modulated per corner by trim. */
  _applyRestLengths() {
    for (const s of this.springs) {
      let wsum = 0, tsum = 0;
      for (let k = 0; k < 4; k++) {
        const w = this._cornerWeight[k][s.a] + this._cornerWeight[k][s.b];
        wsum += w;
        tsum += w * this.corners[k].trim;
      }
      const trim = wsum > 1e-9 ? tsum / wsum : 1;
      s.rest = s.restBase / (this.tension * trim);
    }
  }

  /** Pinned corners are infinite-mass so springs stretch honestly against them. */
  _repin(t) {
    for (let k = 0; k < 4; k++) {
      const c = this.corners[k];
      const ci = this.cornerIdx[k];
      if (c.broken) {
        this.invMass[ci] = 1 / this.nodeMass; // freed node — flogging falls out of this
        continue;
      }
      this.invMass[ci] = 0;
      const p = this._anchorPos(c.anchor, t);
      this.pos[ci * 3] = p.x; this.pos[ci * 3 + 1] = p.y; this.pos[ci * 3 + 2] = p.z;
      this.prev[ci * 3] = p.x; this.prev[ci * 3 + 1] = p.y; this.prev[ci * 3 + 2] = p.z;
    }
  }

  _anchorPos(a, t) {
    const p = a.pos;
    if (!a.sway) return p;
    const s = a.sway(t);
    return { x: p.x + s.x, y: p.y + s.y, z: p.z + s.z };
  }

  _surfaceArea() {
    let total = 0;
    for (let i = 0; i < this.tris.length; i += 3) {
      const a = this.tris[i] * 3, b = this.tris[i + 1] * 3, c = this.tris[i + 2] * 3;
      const e1x = this.pos[b] - this.pos[a], e1y = this.pos[b + 1] - this.pos[a + 1], e1z = this.pos[b + 2] - this.pos[a + 2];
      const e2x = this.pos[c] - this.pos[a], e2y = this.pos[c + 1] - this.pos[a + 1], e2z = this.pos[c + 2] - this.pos[a + 2];
      const nx = e1y * e2z - e1z * e2y, ny = e1z * e2x - e1x * e2z, nz = e1x * e2y - e1y * e2x;
      total += Math.hypot(nx, ny, nz) * 0.5;
    }
    return total;
  }

  /**
   * Advance the sim. Accumulates real time and burns it in fixed SIM_DT chunks,
   * so a variable-rate render loop and a fast-forwarded selftest produce
   * identical traces. Never reads a clock.
   *
   * @param {number} dt    seconds elapsed since last call
   * @param {object} wind  { sample(pos, t) -> {x,y,z} }
   * @param {number} t     world time, seconds
   */
  step(dt, wind, t) {
    if (!this.rigged) return;
    this._acc += dt;
    let n = 0;
    while (this._acc >= SIM_DT && n < MAX_SUBSTEPS) {
      this._substep(SIM_DT, wind, this.t);
      this._acc -= SIM_DT;
      this.t += SIM_DT;
      n++;
    }
    if (n === MAX_SUBSTEPS) this._acc = 0; // dropped frames: don't try to catch up
  }

  _substep(dt, wind, t) {
    this._accumulateWind(wind, t, dt);
    this._integrate(dt);
    this.lambda.fill(0); // XPBD multipliers are per-substep
    for (let i = 0; i < RELAX_ITERS; i++) this._relax(dt * dt);
    this._pinCorners(t);
    this._measureLoads(dt);
    this._checkFailure(dt);
  }

  /** Wind force per FACE — the hypar mechanic lives here. */
  _accumulateWind(wind, t, dt) {
    const pos = this.pos, prev = this.prev, F = this.force;
    F.fill(0);
    const coeff = PRESSURE_COEFF * (1 - this.porosity);
    const tanCoeff = TANGENT_COEFF * (1 - this.porosity);
    const invDt = 1 / dt;
    const probe = this._probe;

    for (let i = 0; i < this.tris.length; i += 3) {
      const ia = this.tris[i] * 3, ib = this.tris[i + 1] * 3, ic = this.tris[i + 2] * 3;

      const e1x = pos[ib] - pos[ia], e1y = pos[ib + 1] - pos[ia + 1], e1z = pos[ib + 2] - pos[ia + 2];
      const e2x = pos[ic] - pos[ia], e2y = pos[ic + 1] - pos[ia + 1], e2z = pos[ic + 2] - pos[ia + 2];
      // |cross| is twice the area and its direction is the face normal
      let nx = e1y * e2z - e1z * e2y, ny = e1z * e2x - e1x * e2z, nz = e1x * e2y - e1y * e2x;
      const len = Math.hypot(nx, ny, nz);
      if (len < 1e-9) continue;
      const area = len * 0.5;
      nx /= len; ny /= len; nz /= len;

      probe.x = (pos[ia] + pos[ib] + pos[ic]) / 3;
      probe.y = (pos[ia + 1] + pos[ib + 1] + pos[ic + 1]) / 3;
      probe.z = (pos[ia + 2] + pos[ib + 2] + pos[ic + 2]) / 3;
      const w = wind.sample(probe, t);

      // Relative wind, not absolute: as the cloth accelerates downwind the load
      // bleeds off by itself. This is what stops flogging from exploding.
      const vx = (pos[ia] - prev[ia] + pos[ib] - prev[ib] + pos[ic] - prev[ic]) / 3 * invDt;
      const vy = (pos[ia + 1] - prev[ia + 1] + pos[ib + 1] - prev[ib + 1] + pos[ic + 1] - prev[ic + 1]) / 3 * invDt;
      const vz = (pos[ia + 2] - prev[ia + 2] + pos[ib + 2] - prev[ib + 2] + pos[ic + 2] - prev[ic + 2]) / 3 * invDt;
      const rx = w.x - vx, ry = w.y - vy, rz = w.z - vz;

      const d = clamp(rx * nx + ry * ny + rz * nz, -MAX_NORMAL_SPEED, MAX_NORMAL_SPEED);
      // d*|d| rather than d^2: keeps the v^2 magnitude but points the force the
      // way the wind is actually blowing. A sail is double-sided.
      const p = coeff * area * d * Math.abs(d);

      const tx = (rx - nx * d) * tanCoeff * area;
      const ty = (ry - ny * d) * tanCoeff * area;
      const tz = (rz - nz * d) * tanCoeff * area;

      const fx = (nx * p + tx) / 3, fy = (ny * p + ty) / 3, fz = (nz * p + tz) / 3;
      F[ia] += fx; F[ia + 1] += fy; F[ia + 2] += fz;
      F[ib] += fx; F[ib + 1] += fy; F[ib + 2] += fz;
      F[ic] += fx; F[ic + 1] += fy; F[ic + 2] += fz;
    }
  }

  _integrate(dt) {
    const pos = this.pos, prev = this.prev, F = this.force, im = this.invMass;
    const dt2 = dt * dt;
    for (let n = 0; n < im.length; n++) {
      if (im[n] === 0) continue; // pinned
      const i = n * 3;
      for (let k = 0; k < 3; k++) {
        const a = F[i + k] * im[n] + (k === 1 ? GRAVITY : 0);
        const x = pos[i + k];
        const nx = x + (x - prev[i + k]) * VEL_DAMP + a * dt2;
        prev[i + k] = x;
        pos[i + k] = nx;
      }
    }
  }

  /**
   * XPBD constraint solve. The plain-PBD version of this is simpler, but its
   * position corrections carry no force information — the leftover stretch
   * after a fixed iteration count is solver error, not fabric strain, so
   * reading load off it measures the solver. XPBD's Lagrange multiplier lambda
   * is the real constraint impulse, so lambda/dt^2 is a genuine newton value
   * and the corner reactions balance the applied wind by construction.
   */
  _relax(dt2) {
    const pos = this.pos, im = this.invMass, lam = this.lambda;
    for (let si = 0; si < this.springs.length; si++) {
      const s = this.springs[si];
      const wa = im[s.a], wb = im[s.b];
      const w = wa + wb;
      if (w === 0) continue; // both ends pinned
      const ia = s.a * 3, ib = s.b * 3;
      const dx = pos[ib] - pos[ia], dy = pos[ib + 1] - pos[ia + 1], dz = pos[ib + 2] - pos[ia + 2];
      const d = Math.hypot(dx, dy, dz);
      if (d < 1e-9) continue;
      const C = d - s.rest;
      const compliance = s.kind === 'bend' ? COMP_BEND : C > 0 ? COMP_STRETCH : COMP_COMPRESS;
      const at = compliance / dt2;
      const dLambda = (-C - at * lam[si]) / (w + at);
      lam[si] += dLambda;
      // grad C is -n for node a and +n for node b, with n = (b - a)/d
      const nx = dx / d, ny = dy / d, nz = dz / d;
      pos[ia] -= nx * dLambda * wa; pos[ia + 1] -= ny * dLambda * wa; pos[ia + 2] -= nz * dLambda * wa;
      pos[ib] += nx * dLambda * wb; pos[ib + 1] += ny * dLambda * wb; pos[ib + 2] += nz * dLambda * wb;
    }
  }

  _pinCorners(t) {
    for (let k = 0; k < 4; k++) {
      const c = this.corners[k];
      if (c.broken) continue;
      const ci = this.cornerIdx[k] * 3;
      const p = this._anchorPos(c.anchor, t);
      this.pos[ci] = p.x; this.pos[ci + 1] = p.y; this.pos[ci + 2] = p.z;
    }
  }

  /**
   * Corner load = magnitude of the VECTOR sum of the tensions in the springs
   * meeting that corner, each read off its XPBD multiplier as |lambda|/dt^2.
   *
   * Reading tension as FABRIC_K * leftover-stretch instead looks equivalent and
   * is not: after a fixed 5 iterations the leftover stretch is solver error, so
   * that number measures the solver rather than the fabric and comes out ~50x
   * hot. The multiplier is the actual constraint impulse, which is why the
   * statics assert balances.
   *
   * The vector sum (rather than a scalar total) is what
   * makes DESIGN.md's anchor-angle mechanic fall out for free: edges pulling in
   * nearly the same direction add up, edges pulling apart partly cancel — so a
   * pinched corner really does multiply its own load.
   */
  _measureLoads(dt) {
    const pos = this.pos, lam = this.lambda;
    const invDt2 = 1 / (dt * dt);
    const alpha = 1 - Math.exp(-dt / LOAD_TAU);
    for (let k = 0; k < 4; k++) {
      const c = this.corners[k];
      if (c.broken) { c.load = 0; c.loadVec.x = c.loadVec.y = c.loadVec.z = 0; continue; }
      const ci = this.cornerIdx[k], cix = ci * 3;
      let sx = 0, sy = 0, sz = 0;
      for (const { s, si } of this.cornerSprings[k]) {
        if (lam[si] >= 0) continue; // slack or compressed fabric pulls on nothing
        const tension = -lam[si] * invDt2; // the multiplier IS the impulse; /dt^2 makes it newtons
        const o = (s.a === ci ? s.b : s.a) * 3;
        const dx = pos[o] - pos[cix], dy = pos[o + 1] - pos[cix + 1], dz = pos[o + 2] - pos[cix + 2];
        const d = Math.hypot(dx, dy, dz);
        if (d < 1e-9) continue;
        sx += (dx / d) * tension; sy += (dy / d) * tension; sz += (dz / d) * tension;
      }
      c.loadVec.x += (sx - c.loadVec.x) * alpha;
      c.loadVec.y += (sy - c.loadVec.y) * alpha;
      c.loadVec.z += (sz - c.loadVec.z) * alpha;
      const raw = Math.hypot(sx, sy, sz);
      c.load += (raw - c.load) * alpha;
      if (c.load > c.peakLoad) c.peakLoad = c.load;
    }
  }

  /** Ported from the prototype: 0.4 s sustained over the rating and it lets go. */
  _checkFailure(dt) {
    for (const c of this.corners) {
      if (c.broken) continue;
      if (c.load > c.hw.rating) c.overload += dt;
      else c.overload = Math.max(0, c.overload - dt * OVERLOAD_RECOVER);
      if (c.overload > OVERLOAD_SECS) {
        c.broken = true;
        c.overload = 0;
        c.load = 0;
        this.events.push({ type: 'break', corner: c, anchorId: c.anchorId, hw: c.hw.name, t: this.t });
      }
    }
    if (this._dirtyRest) { this._applyRestLengths(); this._dirtyRest = false; }
  }

  /** Re-rig a blown corner with fresh hardware. Lane D's hold-E repair calls this. */
  repairCorner(index, hw = HARDWARE[1]) {
    const c = this.corners[index];
    if (!c || !c.broken) return false;
    c.broken = false;
    c.hw = hw;
    c.load = 0;
    c.overload = 0;
    this._repin(this.t);
    this.events.push({ type: 'repair', corner: c, anchorId: c.anchorId, hw: hw.name, t: this.t });
    return true;
  }

  /** Turnbuckle trim at ONE corner (Lane D, 1.2 s hold). Tightens/eases locally. */
  trimCorner(index, delta) {
    const c = this.corners[index];
    if (!c) return false;
    c.trim = clamp(c.trim + delta, TRIM_MIN, TRIM_MAX);
    this._dirtyRest = true;
    return true;
  }

  setTension(tension) {
    this.tension = clamp(tension, TENSION_MIN, TENSION_MAX);
    if (this.rigged) this._applyRestLengths();
  }

  drainEvents() {
    const out = this.events;
    this.events = [];
    return out;
  }

  /**
   * Ground-projected shade over a rect: fraction of sample points on the rect
   * that the sail blocks from the sun. This IS the shade mechanic, so it
   * raycasts toward the actual sun rather than projecting straight down —
   * which is what lets DESIGN.md's moving/seasonal sun change the answer.
   *
   * @param {object} rect    {x, z, w, d} on the ground, metres, world XZ
   * @param {object} sunDir  direction TO the sun; defaults to straight overhead
   */
  coverageOver(rect, sunDir = { x: 0, y: 1, z: 0 }) {
    if (!this.rigged) return 0;
    const len = Math.hypot(sunDir.x, sunDir.y, sunDir.z) || 1;
    const dx = sunDir.x / len, dy = sunDir.y / len, dz = sunDir.z / len;
    if (dy <= 0.01) return 0; // sun at or below the horizon casts no useful shade

    const COLS = 6, ROWS = 4; // prototype sampled 6x4 over the garden
    let hit = 0;
    for (let i = 0; i < COLS; i++) {
      for (let j = 0; j < ROWS; j++) {
        const ox = rect.x + ((i + 0.5) / COLS) * rect.w;
        const oz = rect.z + ((j + 0.5) / ROWS) * rect.d;
        if (this._rayHitsSail(ox, 0, oz, dx, dy, dz)) hit++;
      }
    }
    return hit / (COLS * ROWS);
  }

  /** Moller-Trumbore against every face; 162 tris, cheap enough to not bother accelerating. */
  _rayHitsSail(ox, oy, oz, dx, dy, dz) {
    const pos = this.pos;
    for (let i = 0; i < this.tris.length; i += 3) {
      const a = this.tris[i] * 3, b = this.tris[i + 1] * 3, c = this.tris[i + 2] * 3;
      const e1x = pos[b] - pos[a], e1y = pos[b + 1] - pos[a + 1], e1z = pos[b + 2] - pos[a + 2];
      const e2x = pos[c] - pos[a], e2y = pos[c + 1] - pos[a + 1], e2z = pos[c + 2] - pos[a + 2];
      const px = dy * e2z - dz * e2y, py = dz * e2x - dx * e2z, pz = dx * e2y - dy * e2x;
      const det = e1x * px + e1y * py + e1z * pz;
      if (Math.abs(det) < 1e-9) continue; // ray parallel to the face
      const inv = 1 / det;
      const tx = ox - pos[a], ty = oy - pos[a + 1], tz = oz - pos[a + 2];
      const u = (tx * px + ty * py + tz * pz) * inv;
      if (u < 0 || u > 1) continue;
      const qx = ty * e1z - tz * e1y, qy = tz * e1x - tx * e1z, qz = tx * e1y - ty * e1x;
      const v = (dx * qx + dy * qy + dz * qz) * inv;
      if (v < 0 || u + v > 1) continue;
      const hitT = (e2x * qx + e2y * qy + e2z * qz) * inv;
      if (hitT > 1e-6) return true;
    }
    return false;
  }

  /** Sum of the aerodynamic + weight force on the whole sail, N. Used by the statics assert. */
  netAppliedForce(wind, t) {
    this._accumulateWind(wind, t, SIM_DT);
    let fx = 0, fy = 0, fz = 0;
    for (let n = 0; n < this.invMass.length; n++) {
      fx += this.force[n * 3];
      fy += this.force[n * 3 + 1] + GRAVITY * this.nodeMass;
      fz += this.force[n * 3 + 2];
    }
    return { x: fx, y: fy, z: fz };
  }

  maxLoad() {
    let m = 0;
    for (const c of this.corners) if (c.load > m) m = c.load;
    return m;
  }
}

/**
 * three.js view over a rig. Imported lazily so the sim core above stays
 * headless-runnable; call this only from the browser, after Lane A's vendor/
 * exists. Returns a THREE.Group to add to the scene, with update() per frame.
 */
export async function createSailView(rig, { color = 0xd8c48a } = {}) {
  const THREE = await import('../vendor/three.module.js');

  const geo = new THREE.BufferGeometry();
  const verts = new Float32Array(rig.pos.length);
  geo.setAttribute('position', new THREE.BufferAttribute(verts, 3));
  geo.setIndex(new THREE.BufferAttribute(new Uint16Array(rig.tris), 1));

  const mat = new THREE.MeshStandardMaterial({
    color, side: THREE.DoubleSide, roughness: 0.92, metalness: 0.0,
  });
  const mesh = new THREE.Mesh(geo, mat);
  mesh.castShadow = true;    // the shadow IS the product
  mesh.receiveShadow = true;
  mesh.frustumCulled = false; // it flogs well outside its initial bounds

  const group = new THREE.Group();
  group.add(mesh);
  group.update = () => {
    verts.set(rig.pos);
    geo.attributes.position.needsUpdate = true;
    geo.computeVertexNormals();
    geo.computeBoundingSphere();
  };
  group.update();
  return group;
}