Railway cant, cant ramps and track quantities

const G = 9.81;
const KMH_PER_MPH = 1.609344;
const M_PER_FT = 1200 / 3937;   // US survey foot, as units.js

export const kmh = (v, unit) => (fin(v) && v > 0 ? (unit === 'mph' ? v * KMH_PER_MPH : v) : null);
/** Project length → metres (US survey feet for ft projects). */
export const toM = (v, units) => (units === 'm' ? v : v * M_PER_FT);
/** A length in metres (a cant in mm ÷ a gradient in mm/m) → project units. */
const toProj = (m, units) => (units === 'm' ? m : m / M_PER_FT);
/** mm → project units. */
export const mmTo = (v, units) => (units === 'm' ? v / 1000 : v / 1000 / M_PER_FT);

export function equilibriumCant(vKmh, Rm, base = 1500) {
  if (!(vKmh > 0) || !(Rm > 0)) return 0;
  const v = vKmh / 3.6;
  return (base * v * v) / (G * Rm);
}

/** The cant ramp's shape g(u) (0 → 1) for a transition kind: the cant follows the curvature. */
export function rampShape(kind) {
  switch (kind) {
    case 'bloss': return (u) => u * u * (3 - 2 * u);
    case 'helmert': return (u) => (u <= 0.5 ? 2 * u * u : 1 - 2 * (1 - u) ** 2);
    case 'sinehalf': return (u) => (1 - Math.cos(Math.PI * u)) / 2;
    case 'sine': return (u) => u - Math.sin(2 * Math.PI * u) / (2 * Math.PI);
    default: return (u) => u;   // clothoid, cubic parabola, a ramp on the straight
  }
}

export const RAMP_PEAK = { linear: 1, bloss: 1.5, helmert: 2, sinehalf: Math.PI / 2, sine: 2 };

/** Each rail's height (mm) over the profile for a signed cant c (+ = the right rail higher) on a curve whose outer side is `high`. */
function heights(E, high, pivot) {
  const c = high === 'R' ? E : -E;
  if (pivot === 'centre') return { L: -c / 2, R: c / 2 };
  if (pivot === 'inner') return high === 'R' ? { L: 0, R: E } : { L: E, R: 0 };   // the inner rail stays; a negative cant lowers the outer
  return { L: Math.max(0, -c), R: Math.max(0, c) };   // 'low': the lower rail stays on the profile
}

const ramp = (end, on, from, to) => {
  const len = to - from, shape = on ? SHAPE_OF[kind] || 'linear' : 'linear';
  const gradient = len > 1e-9 ? (Math.abs(E) / toM(len, units)) * RAMP_PEAK[shape] : Infinity;
  return { end, from, to, len, shape, onTransition: on, gradient, rate: v ? gradient * v : null };
};
rec.rampIn = c.ls1 > 0 ? ramp('in', true, c.ts, c.sc) : ramp('in', false, c.ts - toProj(Math.abs(E) / L.gradient, units), c.ts);
rec.rampOut = c.ls2 > 0 ? ramp('out', true, c.cs, c.st) : ramp('out', false, c.st, c.st + toProj(Math.abs(E) / L.gradient, units));

if (Math.abs(E) > L.cant + 1e-9) warn('rail.cant.max', 'warning', part, { pi: c.piIndex + 1, cant: n1(E), max: n1(L.cant) }, { limit: L.cant, value: Math.abs(E), unit: 'mm' });
if (rec.deficiency !== null && rec.deficiency > L.deficiency + 1e-9) warn('rail.cant.deficiency', 'warning', part, { pi: c.piIndex + 1, d: n1(rec.deficiency), max: n1(L.deficiency), eq: n1(eq), cant: n1(E) }, { limit: L.deficiency, value: rec.deficiency, unit: 'mm' });
if (rec.excess !== null && rec.excess > L.excess + 1e-9) warn('rail.cant.excess', 'warning', part, { pi: c.piIndex + 1, x: n1(rec.excess), max: n1(L.excess), eq: n1(eqMin), cant: n1(E) }, { limit: L.excess, value: rec.excess, unit: 'mm' });

const seg = (from, to, type, e0, e1, high) => {
  const a = Math.max(from, al.start), b = Math.min(to, al.end);
  if (!(b - a > 1e-9)) return;
  const h0 = heights(e0, high, plan.pivot), h1 = heights(e1, high, plan.pivot);
  segments.push({ start: a, length: b - a, type, startCant: e0, endCant: e1, startLeft: h0.L, startRight: h0.R, endLeft: h1.L, endRight: h1.R, high });
};
for (const p of pieces) {
  const from = Math.max(p.from, at);
  if (from > at + 1e-9) seg(at, from, 'CONSTANTCANT', 0, 0, p.c.high);
  if (p.to - from > 1e-9) {
    // a ramp cut short by the one before starts where that one ended, at the cant this ramp has there
    const g = p.shape ? rampShape(p.shape) : null, val = (s) => (g ? p.e0 + (p.e1 - p.e0) * g(Math.min(1, Math.max(0, (s - p.from) / (p.to - p.from)))) : p.e0);
    seg(from, p.to, p.shape ? CANT_TYPES[p.shape] : 'CONSTANTCANT', val(from), val(p.to), p.c.high);
  }
  at = Math.max(at, p.to);
}
if (al.end > at + 1e-9) seg(at, al.end, 'CONSTANTCANT', 0, 0, curves[curves.length - 1]?.high || 'L');

/** The applied cant and each rail's height at a station. */
export function cantAt(plan, sta) {
  for (const s of plan?.segments || []) {
    if (sta < s.start - 1e-9 || sta > s.start + s.length + 1e-9) continue;
    if (s.type === 'CONSTANTCANT') return { E: s.startCant, L: s.startLeft, R: s.startRight };
    const shape = Object.keys(CANT_TYPES).find((k) => CANT_TYPES[k] === s.type), g = rampShape(shape);
    // the segment may be a ramp cut short: map back to the ramp's own u through its curve record
    const c = plan.curves.find((x) => x.rampIn && sta >= x.rampIn.from - 1e-9 && sta <= x.rampIn.to + 1e-9 && CANT_TYPES[x.rampIn.shape] === s.type) || plan.curves.find((x) => x.rampOut && sta >= x.rampOut.from - 1e-9 && sta <= x.rampOut.to + 1e-9 && CANT_TYPES[x.rampOut.shape] === s.type);
    let E;
    if (c) {
      const inn = c.rampIn && sta >= c.rampIn.from - 1e-9 && sta <= c.rampIn.to + 1e-9, r = inn ? c.rampIn : c.rampOut, u = (sta - r.from) / Math.max(r.len, 1e-12);
      E = c.cant * (inn ? g(u) : 1 - g(u));
    } else E = s.startCant + ((s.endCant - s.startCant) * (sta - s.start)) / Math.max(s.length, 1e-12);
    const h = heights(E, s.high, plan.pivot);
    return { E, L: h.L, R: h.R };
  }
  return { E: 0, L: 0, R: 0 };
}

/**
 * The track's cross section square to the alignment (offsets o right +, z from the profile's elevation; project units). The rails'
 * heads sit at hL, hR (mm) over the profile, ±base/2 from each track's centre; the sleeper lies parallel to them, the rail + pad
 * below; the ballast bed is filled to `crib` under the sleepers' tops, its shoulder level out to `shoulder` past the sleeper end, then the
 * ballast slope down to the sub-ballast, whose top falls away from the centre at `fall` and lies `ballast.depth` under the sleeper
 * at the lowest low rail; the sub-ballast runs `cess` past the ballast toe, `sub.depth` thick. The formation (its bottom) is not
 * tilted by the cant; the track is.
 */
export function trackShape(rail, { hL = 0, hR = 0 } = {}, units = 'm') {
  const mm = (v) => mmTo(v, units), P = RAIL_PROFILES[rail.rail] || RAIL_PROFILES['60E1'];
  const base = mm(baseOf(rail)), n = rail.tracks || 1, S = rail.sleeper, B = rail.ballast, SB = rail.sub;
  const below = mm(P.h + (rail.pad || 0));   // rail head to the sleeper's top
  const slope = (mm(hR) - mm(hL)) / base;    // dz / do across the track
  const tracks = [];
  for (let k = 0; k < n; k++) {
    const c = (k - (n - 1) / 2) * rail.centres;
    const zl = mm(hL), plane = (o) => zl + (o - (c - base / 2)) * slope - below;   // the sleeper's top
    const a = c - S.length / 2, b = c + S.length / 2;
    const lowO = hL <= hR ? c - base / 2 : c + base / 2;
    tracks.push({ c, slope, below, rails: [{ o: c - base / 2, z: mm(hL) }, { o: c + base / 2, z: mm(hR) }], sleeper: [{ o: a, z: plane(a) }, { o: b, z: plane(b) }, { o: b, z: plane(b) - S.depth }, { o: a, z: plane(a) - S.depth }], plane, lowO, under: plane(lowO) - S.depth });
  }
  // the sub-ballast's top: z0 − |o|·fall, ballast.depth under the sleeper at each track's low rail (the deepest one decides)
  const z0 = Math.min(...tracks.map((k) => k.under - B.depth + Math.abs(k.lowO) * SB.fall));
  const sbTop = (o) => z0 - Math.abs(o) * SB.fall;
  // the ballast top from the centre out, each side
  const crib = B.crib || 0;   // the ballast's top lies this far under the sleepers' tops (they stand proud of it)
  const ends = tracks.flatMap((k) => [{ o: k.sleeper[0].o, z: k.sleeper[0].z - crib }, { o: k.sleeper[1].o, z: k.sleeper[1].z - crib }]);
  const topAt0 = () => {
    const on = tracks.find((k) => k.sleeper[0].o <= 1e-9 && k.sleeper[1].o >= -1e-9);
    if (on) return on.plane(0) - crib;
    const l = ends.filter((e) => e.o < 0).sort((p, q) => q.o - p.o)[0], r = ends.filter((e) => e.o > 0).sort((p, q) => p.o - q.o)[0];
    return l && r ? l.z + ((r.z - l.z) * (0 - l.o)) / (r.o - l.o) : (l || r).z;
  };
  const top = {}, tz = topAt0();
  for (const [side, sx] of [['L', -1], ['R', 1]]) {
    const out = [{ o: 0, z: tz, code: 'CL' }];
    const mine = ends.filter((e) => e.o * sx > 1e-9).sort((p, q) => Math.abs(p.o) - Math.abs(q.o));
    mine.forEach((e, i) => out.push({ o: e.o, z: e.z, code: i ? 'SLP' + (i + 1) : 'SLP' }));
    const last = mine.length ? mine[mine.length - 1] : { o: 0, z: tz };
    const bsh = { o: last.o + sx * B.shoulder, z: last.z, code: 'BSH' };
    // down the ballast slope to the sub-ballast top: bsh.z − d/slope = z0 − (|bsh.o| + d)·fall
    const d = Math.max(0, (bsh.z - z0 + Math.abs(bsh.o) * SB.fall) / Math.max(1 / B.slope - SB.fall, 1e-6));
    const toe = { o: bsh.o + sx * d, z: sbTop(bsh.o + sx * d), code: 'BTOE' };
    const cess = { o: toe.o + sx * SB.cess, z: sbTop(toe.o + sx * SB.cess), code: 'CESS' };
    out.push(bsh, toe, cess);
    top[side] = out;
  }
  // the ballast: its top (toe, shoulder, sleeper ends …, shoulder, toe), its bottom along the sub-ballast top (with the crown)
  const Lr = top.L.slice().reverse(), Rr = top.R;
  const tL = Lr.filter((q) => q.code !== 'CESS'), tR = Rr.filter((q) => q.code !== 'CESS' && q.code !== 'CL');
  const ballast = [...tL.map(({ o, z }) => ({ o, z })), ...tR.map(({ o, z }) => ({ o, z }))];
  const toeL = top.L.find((q) => q.code === 'BTOE'), toeR = top.R.find((q) => q.code === 'BTOE');
  ballast.push({ o: 0, z: z0 });   // the crown of the sub-ballast under it (after the right toe, before closing at the left toe)
  const cL = top.L.find((q) => q.code === 'CESS'), cR = top.R.find((q) => q.code === 'CESS');
  const sub = [{ o: cL.o, z: cL.z }, { o: 0, z: z0 }, { o: cR.o, z: cR.z }, { o: cR.o, z: cR.z - SB.depth }, { o: 0, z: z0 - SB.depth }, { o: cL.o, z: cL.z - SB.depth }];
  void toeL; void toeR;
  return { tracks: tracks.map(({ plane, lowO, under, ...k }) => k), top, ballast, sub, areas: { ballast: Math.abs(area(ballast)), sub: Math.abs(area(sub)) }, z0 };
}

function area(P) { let s = 0; for (let i = 0; i < P.length; i++) { const a = P[i], b = P[(i + 1) % P.length]; s += a.o * b.z - b.o * a.z; } return s / 2; }

/** Quantities from the corridor's sections (each with `track` from corridor.js): average end areas between sections, track and rail
 *  lengths along each track's own line (shorter inside a curve), sleepers every `spacing` along each track. */
export function trackQuantities(sections, rail, units = 'm') {
  const secs = (sections || []).filter((s) => s.track);
  const n = rail?.tracks || 1, out = { length: 0, tracks: [], sleepers: 0, railLength: 0, sleeperVolume: 0, ballastGross: 0, ballast: 0, sub: 0 };
  if (secs.length < 2) return out;
  out.length = secs[secs.length - 1].sta - secs[0].sta;
  for (let i = 1; i < secs.length; i++) {
    const a = secs[i - 1].track, b = secs[i].track, ds = secs[i].sta - secs[i - 1].sta;
    out.ballastGross += ((a.shape.areas.ballast + b.shape.areas.ballast) / 2) * ds;
    out.sub += ((a.shape.areas.sub + b.shape.areas.sub) / 2) * ds;
  }
  const S = rail.sleeper, world = (sec, o) => ({ x: sec.track.p.x + sec.track.r.x * o, y: sec.track.p.y + sec.track.r.y * o });
  for (let k = 0; k < n; k++) {
    let len = 0, rl = 0;
    for (let i = 1; i < secs.length; i++) {
      const A = secs[i - 1].track.shape.tracks[k], B = secs[i].track.shape.tracks[k];
      const p = world(secs[i - 1], A.c), q = world(secs[i], B.c);
      len += Math.hypot(q.x - p.x, q.y - p.y);
      for (let j = 0; j < 2; j++) { const u = world(secs[i - 1], A.rails[j].o), w = world(secs[i], B.rails[j].o); rl += Math.hypot(w.x - u.x, w.y - u.y); }
    }
    const sleepers = Math.floor(len / S.spacing + 1e-9) + 1;
    out.tracks.push({ c: secs[0].track.shape.tracks[k].c, length: len, sleepers, rail: rl });
    out.sleepers += sleepers; out.railLength += rl;
  }
  out.sleeperVolume = out.sleepers * S.length * S.width * S.depth;
  out.ballast = Math.max(0, out.ballastGross - out.sleepers * S.length * S.width * Math.max(0, S.depth - (rail.ballast?.crib || 0)));   // the part of each sleeper in the bed
  void units;
  return out;
}