Add 3D sail cloth sim with per-face wind and XPBD corner loads

Verlet cloth on a bilinear patch between 4 anchors, N=10 grid,
structural/shear/bend constraints, 5 relaxation iterations at a fixed
1/60 substep. Wind is applied per FACE so hypar twist genuinely sheds
load rather than being cosmetic.

Two deviations from PLAN3D worth flagging:

- Load is read from each constraint's XPBD Lagrange multiplier, not from
  FABRIC_K * leftover-stretch. After a fixed iteration count the leftover
  stretch is solver error, not fabric strain, so the naive reading came
  out ~50x hot (60 kN peaks on a 5x5 m sail). The multiplier is the real
  constraint impulse, which the statics assert confirms by balancing the
  corner reactions against the applied wind to 8%.

- Wind uses a signed square (d*|d|) rather than clamp(d)^2, so the
  leeward face is pushed too. A sail is double-sided.

The sim core deliberately does not import three.js: it runs headless
under node today, stays allocation-free in the hot loop, and replays
bit-for-bit. createSailView() pulls three in lazily for rendering.

Loads land in real newtons (~1-4 kN on a 5x5 m sail in a 34 m/s storm),
so hardware ratings are real working load limits.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
m3ultra 2026-07-16 21:28:59 +10:00
parent 8d76340f49
commit 06ec4cbea2
2 changed files with 969 additions and 0 deletions

617
web/world/js/sail.js Normal file
View File

@ -0,0 +1,617 @@
/**
* 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;
}

View File

@ -0,0 +1,352 @@
/**
* sail.selftest.js assert suite for the sail sim. [Lane B]
*
* DOM-free and three-free on purpose: runs today under `node sail.selftest.js`
* before Lane A's shell exists, and Lane A's selftest.html can import
* runSailSelftest() unchanged once it lands. Drives the sim with fixed-dt loops
* only never rAF, never a clock.
*
* The headline assert is `hypar sheds load vs flat`: it is the game's thesis
* stated as a test. If it ever goes red, the sail has stopped being a sail.
*/
import { SailRig, HARDWARE } from './sail.js';
const SIM_DT = 1 / 60;
// ---------- deterministic stub wind (Lane C owns the real weather.js) ----------
function mulberry32(seed) {
return function () {
seed |= 0; seed = (seed + 0x6d2b79f5) | 0;
let t = Math.imul(seed ^ (seed >>> 15), 1 | seed);
t = (t + Math.imul(t ^ (t >>> 7), 61 | t)) ^ t;
return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
};
}
/**
* Ports the prototype's gust scheduler shape: telegraph 1.5 s, ramp 0.8 s,
* hold ~1.7 s, fade 1 s. The schedule is precomputed from the seed so
* sample(pos, t) stays a pure function of t which is what makes the
* determinism assert meaningful and lets the selftest fast-forward.
*/
function makeStubWind({ seed = 7, stormLen = 90, dir = { x: 0, y: 0, z: 1 }, calm = false } = {}) {
const rng = mulberry32(seed);
const gusts = [];
for (let t = 3; t < stormLen; t += 5 + rng() * 7) {
gusts.push({ start: t, pow: 12 + rng() * 16 + 10 * (t / stormLen) });
}
const len = Math.hypot(dir.x, dir.y, dir.z) || 1;
const dx = dir.x / len, dy = dir.y / len, dz = dir.z / len;
const out = { x: 0, y: 0, z: 0 };
return {
speedAt(t) {
if (calm) return 0;
const p = Math.min(1, t / stormLen);
let speed = 8 + 26 * Math.min(1, p * 1.6);
for (const g of gusts) {
const gt = t - g.start;
if (gt < 0 || gt >= 5) continue;
if (gt < 1.5) continue; // telegraph: you see it, you don't feel it
else if (gt < 2.3) speed += g.pow * (gt - 1.5) / 0.8; // ramp
else if (gt < 4.0) speed += g.pow; // hold
else speed += g.pow * (5.0 - gt); // fade
}
return speed;
},
sample(pos, t) {
const s = this.speedAt(t);
out.x = dx * s; out.y = dy * s; out.z = dz * s;
return out;
},
};
}
const constantWind = (v) => ({ sample: () => v, speedAt: () => Math.hypot(v.x, v.y, v.z) });
// ---------- test rigs ----------
// Same 5x5 m footprint, same multiset of corner heights {4.0, 4.0, 2.5, 2.5}.
// Only the ARRANGEMENT differs: coplanar (flat, pitched) vs permuted (twisted
// hypar). Any load difference is therefore purely geometry, nothing else.
const FOOT = [
{ x: -2.5, z: -2.5 }, { x: 2.5, z: -2.5 }, { x: 2.5, z: 2.5 }, { x: -2.5, z: 2.5 },
];
const HEIGHTS_FLAT = [4.0, 4.0, 2.5, 2.5]; // y is linear in z -> one plane
const HEIGHTS_HYPAR = [4.0, 2.5, 4.0, 2.5]; // opposite corners up/down -> saddle
const makeAnchors = (heights) =>
FOOT.map((f, i) => ({
id: `a${i}`,
type: 'post',
pos: { x: f.x, y: heights[i], z: f.z },
}));
const ALL_IDS = ['a0', 'a1', 'a2', 'a3'];
const UNBREAKABLE = { name: 'test rig', cost: 0, rating: Infinity };
function rig(heights, { hw = UNBREAKABLE, tension = 1.0, porosity = 0 } = {}) {
return new SailRig({ anchors: makeAnchors(heights), gridN: 10, porosity })
.attach(ALL_IDS, [hw, hw, hw, hw], tension);
}
/** Fixed-dt fast-forward. Returns peak corner load seen over the whole run, N. */
function runStorm(r, wind, secs, onStep) {
const steps = Math.round(secs / SIM_DT);
let peak = 0;
for (let i = 0; i < steps; i++) {
r.step(SIM_DT, wind, i * SIM_DT);
const m = r.maxLoad();
if (m > peak) peak = m;
if (onStep) onStep(r, i);
}
return peak;
}
// ---------- tiny assert harness ----------
const results = [];
function test(name, fn) {
try {
const detail = fn();
results.push({ name, pass: true, detail: detail || '' });
} catch (e) {
results.push({ name, pass: false, detail: e.message });
}
}
function assert(cond, msg) {
if (!cond) throw new Error(msg);
}
const kN = (n) => `${(n / 1000).toFixed(2)} kN`;
// ---------- the suite ----------
export function runSailSelftest() {
results.length = 0;
test('sim stays finite through a full storm', () => {
const r = rig(HEIGHTS_HYPAR);
runStorm(r, makeStubWind({ stormLen: 90 }), 90);
for (const v of r.pos) assert(Number.isFinite(v), 'node position went NaN/Infinity');
for (const c of r.corners) assert(Number.isFinite(c.load), 'corner load went NaN');
return `peak ${kN(r.corners.reduce((m, c) => Math.max(m, c.peakLoad), 0))}`;
});
test('sail sags under gravity when calm', () => {
const r = rig(HEIGHTS_FLAT);
runStorm(r, makeStubWind({ calm: true }), 6);
const N = r.N, mid = (Math.floor(N / 2) * N + Math.floor(N / 2)) * 3;
const midY = r.pos[mid + 1];
const cornerMeanY = HEIGHTS_FLAT.reduce((a, b) => a + b) / 4;
assert(midY < cornerMeanY, `belly (${midY.toFixed(2)}m) should hang below corner mean (${cornerMeanY}m)`);
return `belly sags ${(cornerMeanY - midY).toFixed(2)} m below corner plane`;
});
// Newton's third law. This is what calibrates FABRIC_K into real newtons:
// if the corner reactions don't sum to the actual aerodynamic + weight force
// on the fabric, the load meter is lying and every kN rating is meaningless.
test('statics: corner reactions balance the applied force', () => {
const w = constantWind({ x: 0, y: 0, z: 18 });
const r = rig(HEIGHTS_FLAT);
runStorm(r, w, 12); // settle
// A membrane in steady wind never fully stops moving, so this compares the
// TIME-AVERAGED reaction against the time-averaged applied force. That is
// the momentum balance that has to hold; instant by instant it need not.
let n = 0, ax = 0, ay = 0, az = 0, rx = 0, ry = 0, rz = 0;
const steps = Math.round(4 / SIM_DT);
for (let i = 0; i < steps; i++) {
const t = 12 + i * SIM_DT;
r.step(SIM_DT, w, t);
const f = r.netAppliedForce(w, t);
ax += f.x; ay += f.y; az += f.z;
for (const c of r.corners) { rx += c.loadVec.x; ry += c.loadVec.y; rz += c.loadVec.z; }
n++;
}
ax /= n; ay /= n; az /= n; rx /= n; ry /= n; rz /= n;
const appliedMag = Math.hypot(ax, ay, az);
const err = Math.hypot(rx - ax, ry - ay, rz - az) / appliedMag;
assert(err < 0.2, `reactions ${kN(Math.hypot(rx, ry, rz))} vs applied ${kN(appliedMag)}${(err * 100).toFixed(0)}% out of balance`);
return `applied ${kN(appliedMag)}, reactions ${kN(Math.hypot(rx, ry, rz))}, residual ${(err * 100).toFixed(1)}%`;
});
// THE THESIS. A twisted sail resists bellying into one coherent pocket, so
// its worst moment is gentler than a flat sail's worst moment.
//
// Scored on WORST CASE over wind direction, not per-direction. Lane C's
// storms veer, so the player never gets to choose the wind, and worst-case is
// what the hardware actually has to survive. Per-direction would be a false
// assert: a flat sail that happens to sit edge-on to the wind genuinely does
// have low drag, and from that one angle it beats the hypar. Demanding
// otherwise would mean tuning the sim into a lie.
test('hypar sheds load vs flat, worst case over wind direction (the thesis)', () => {
const DIRS = [
{ name: 'N', x: 0, z: 1 }, { name: 'NE', x: 0.707, z: 0.707 },
{ name: 'E', x: 1, z: 0 }, { name: 'SE', x: 0.707, z: -0.707 },
{ name: 'S', x: 0, z: -1 }, { name: 'SW', x: -0.707, z: -0.707 },
{ name: 'W', x: -1, z: 0 }, { name: 'NW', x: -0.707, z: 0.707 },
];
const sweep = (heights) => {
let worst = 0, at = '';
for (const d of DIRS) {
const storm = makeStubWind({ seed: 7, stormLen: 45, dir: { x: d.x, y: 0, z: d.z } });
const p = runStorm(rig(heights), storm, 45);
if (p > worst) { worst = p; at = d.name; }
}
return { worst, at };
};
const flat = sweep(HEIGHTS_FLAT);
const hypar = sweep(HEIGHTS_HYPAR);
assert(
hypar.worst < flat.worst * 0.8,
`hypar worst ${kN(hypar.worst)} (${hypar.at}) should be well under flat worst ${kN(flat.worst)} (${flat.at})`
);
const shed = (1 - hypar.worst / flat.worst) * 100;
return `flat worst ${kN(flat.worst)} from ${flat.at} -> hypar worst ${kN(hypar.worst)} from ${hypar.at} (sheds ${shed.toFixed(0)}%)`;
});
test('cascade: losing a corner spikes its neighbours', () => {
const w = constantWind({ x: 0, y: 0, z: 22 });
const r = rig(HEIGHTS_HYPAR);
runStorm(r, w, 6); // settle
const before = Math.max(r.corners[1].load, r.corners[3].load);
r.corners[0].broken = true;
r._repin(r.t);
runStorm(r, w, 2.5); // let the load redistribute
const after = Math.max(r.corners[1].load, r.corners[3].load);
assert(after >= before * 2, `neighbour went ${kN(before)} -> ${kN(after)}, wanted >= 2x`);
return `neighbour ${kN(before)} -> ${kN(after)} (${(after / before).toFixed(1)}x)`;
});
test('determinism: identical inputs give byte-equal load traces', () => {
const trace = () => {
const r = rig(HEIGHTS_HYPAR);
const w = makeStubWind({ seed: 3, stormLen: 30 });
const out = [];
runStorm(r, w, 30, (rr) => { for (const c of rr.corners) out.push(c.load); });
return out;
};
const a = trace(), b = trace();
assert(a.length === b.length, 'traces differ in length');
for (let i = 0; i < a.length; i++) {
assert(a[i] === b[i], `sample ${i} diverged: ${a[i]} vs ${b[i]}`);
}
return `${a.length} load samples identical`;
});
test('determinism: variable frame dt matches fixed dt', () => {
// Lane A's render loop delivers ragged dt. The internal accumulator has to
// absorb that, or nothing the selftest proves applies to the real game.
const w1 = makeStubWind({ seed: 5, stormLen: 20 });
const fixed = rig(HEIGHTS_HYPAR);
for (let i = 0; i < Math.round(20 / SIM_DT); i++) fixed.step(SIM_DT, w1, i * SIM_DT);
const w2 = makeStubWind({ seed: 5, stormLen: 20 });
const ragged = rig(HEIGHTS_HYPAR);
const rng = mulberry32(99);
let acc = 0;
while (acc < 20) {
const dt = 0.004 + rng() * 0.02; // 4-24 ms frames
ragged.step(dt, w2, acc);
acc += dt;
}
for (let k = 0; k < 4; k++) {
const d = Math.abs(fixed.corners[k].load - ragged.corners[k].load);
assert(d < 1e-6, `corner ${k} drifted ${d.toFixed(6)} N between fixed and ragged dt`);
}
return 'ragged frame times converge on the fixed-dt trace';
});
test('tension dial changes load (drum tight shock-loads)', () => {
const w = constantWind({ x: 0, y: 0, z: 20 });
const loose = rig(HEIGHTS_HYPAR, { tension: 0.7 });
const tight = rig(HEIGHTS_HYPAR, { tension: 1.35 });
const loosePeak = runStorm(loose, w, 8);
const tightPeak = runStorm(tight, w, 8);
assert(tightPeak > loosePeak, `tight ${kN(tightPeak)} should exceed loose ${kN(loosePeak)}`);
return `loose ${kN(loosePeak)} vs tight ${kN(tightPeak)}`;
});
test('porous shade cloth carries less load than solid membrane', () => {
const w = constantWind({ x: 0, y: 0, z: 20 });
const solid = runStorm(rig(HEIGHTS_HYPAR, { porosity: 0 }), w, 8);
const porous = runStorm(rig(HEIGHTS_HYPAR, { porosity: 0.35 }), w, 8);
assert(porous < solid, `porous ${kN(porous)} should be under solid ${kN(solid)}`);
return `solid ${kN(solid)} vs porous ${kN(porous)}`;
});
test('coverage: sail shades the ground under it, not beside it', () => {
const r = rig(HEIGHTS_FLAT);
runStorm(r, makeStubWind({ calm: true }), 4);
const under = r.coverageOver({ x: -2, z: -2, w: 4, d: 4 });
const beside = r.coverageOver({ x: 12, z: 12, w: 4, d: 4 });
assert(under > 0.9, `ground under the sail only ${(under * 100).toFixed(0)}% shaded`);
assert(beside === 0, `ground 12 m away reported ${(beside * 100).toFixed(0)}% shaded`);
return `under sail ${(under * 100).toFixed(0)}%, off to the side ${(beside * 100).toFixed(0)}%`;
});
test('coverage tracks a low sun off to the side', () => {
const r = rig(HEIGHTS_FLAT);
runStorm(r, makeStubWind({ calm: true }), 4);
const noon = r.coverageOver({ x: -2, z: -2, w: 4, d: 4 }, { x: 0, y: 1, z: 0 });
const lowSun = r.coverageOver({ x: -2, z: -2, w: 4, d: 4 }, { x: 0.9, y: 0.25, z: 0 });
assert(noon > lowSun, `shadow should slide off the bed as the sun drops (noon ${noon}, low ${lowSun})`);
return `noon ${(noon * 100).toFixed(0)}% -> low sun ${(lowSun * 100).toFixed(0)}%`;
});
// PLAN3D §7 definition of done, in miniature.
test('cheap flat rig cascades; twisted mixed rig survives', () => {
const storm = () => makeStubWind({ seed: 11, stormLen: 90 });
const cheap = rig(HEIGHTS_FLAT, { hw: HARDWARE[0], tension: 1.35 });
runStorm(cheap, storm(), 90);
const cheapBroken = cheap.corners.filter((c) => c.broken).length;
const good = rig(HEIGHTS_HYPAR, { hw: HARDWARE[2], tension: 0.95 });
runStorm(good, storm(), 90);
const goodBroken = good.corners.filter((c) => c.broken).length;
assert(cheapBroken >= 2, `flat drum-tight carabiner rig only lost ${cheapBroken} corners — should cascade`);
assert(goodBroken === 0, `twisted rated-shackle rig lost ${goodBroken} corners — should survive`);
return `cheap flat lost ${cheapBroken}/4, good hypar lost ${goodBroken}/4`;
});
test('repair re-pins a blown corner and it carries load again', () => {
const w = constantWind({ x: 0, y: 0, z: 20 });
const r = rig(HEIGHTS_HYPAR);
runStorm(r, w, 4);
r.corners[0].broken = true;
r._repin(r.t);
runStorm(r, w, 1);
assert(r.corners[0].load === 0, 'broken corner should carry no load');
assert(r.repairCorner(0, HARDWARE[1]), 'repairCorner should report success');
runStorm(r, w, 3);
assert(r.corners[0].load > 100, `repaired corner only pulling ${kN(r.corners[0].load)}`);
const evs = r.drainEvents();
assert(evs.some((e) => e.type === 'repair'), 'no repair event emitted');
return `repaired corner back to ${kN(r.corners[0].load)}`;
});
const pass = results.every((r) => r.pass);
return { pass, results };
}
// ---------- entry points ----------
function report(out) {
const lines = out.results.map(
(r) => `${r.pass ? 'PASS' : 'FAIL'} ${r.name}${r.detail ? `\n ${r.detail}` : ''}`
);
return `${lines.join('\n')}\n\n${out.pass ? 'ALL GREEN' : 'FAILURES'}${out.results.filter((r) => r.pass).length}/${out.results.length}`;
}
// Run only when invoked directly: `node web/world/js/sail.selftest.js`.
// Importing the module must not run the suite. In the browser `process` is
// undefined, so Lane A's selftest.html just calls runSailSelftest() itself.
if (typeof process !== 'undefined' && process.versions?.node && import.meta.filename === process.argv[1]) {
const out = runSailSelftest();
console.log(report(out));
process.exit(out.pass ? 0 : 1);
}
export { report, makeStubWind, HEIGHTS_FLAT, HEIGHTS_HYPAR, makeAnchors };