Storm hydrographs routed through the pipe network

/** One catchment's runoff hydrograph at times 0, dt, 2dt…: C·A converted to flow units. Flow per (area × intensity): ft³/s per ac·in/h (≈1), m³/s per ha·mm/h (1/360). */
export function catchmentHydrograph({ CA, tc, i, duration, dt, steps, units }) {
  const k = units === 'm' ? 1 / 360 : 1;
  const q = new Array(steps).fill(0), w = Math.max(tc, dt);
  for (let n = 0; n < steps; n++) {
    const t1 = n * dt;
    q[n] = k * CA * (rainDepth(i, duration, t1) - rainDepth(i, duration, t1 - w)) / w;
  }
  return q;
}

/** Muskingum coefficients for K, X, dt (minutes). `lag` = the step is too coarse to attenuate: C2 would go negative. */
export function muskingum(K, X, dt) {
  const d = K * (1 - X) + 0.5 * dt;
  const c0 = (-K * X + 0.5 * dt) / d, c1 = (K * X + 0.5 * dt) / d, c2 = (K * (1 - X) - 0.5 * dt) / d;
  return { c0, c1, c2, lag: c0 < 0 || c2 < 0 };
}

/** Route one series through a reach. Muskingum needs 2·K·X ≤ step ≤ 2·K·(1 − X) for its three coefficients to stay positive, so the reach is run on a grid of
 *  max(dt, 2·K·X) minutes (the inflow read off by linear interpolation, the outflow interpolated back to dt); if even that is longer than 2·K·(1 − X) the reach is too short
 *  to attenuate on this grid and the flow passes after a plain lag of K minutes (linear interpolation, peak and volume kept). `step` = the grid used. */
export function routeReach(inflow, K, X, dt) {
  const n = inflow.length, out = new Array(n).fill(0);
  const d = Math.max(dt, 2 * K * X), co = muskingum(K, X, d);
  const at = (a, u) => { if (u <= 0) return a[0]; const j = Math.floor(u), f = u - j; return j >= a.length - 1 ? a[a.length - 1] : a[j] * (1 - f) + a[j + 1] * f; };
  if (co.lag) {
    const s = K / dt;
    for (let j = 0; j < n; j++) { const f = j - s, i = Math.floor(f), b = f - i; out[j] = f < 0 ? 0 : (inflow[i] ?? 0) * (1 - b) + (inflow[i + 1] ?? inflow[i] ?? 0) * b; }
    return { out, lag: true, step: dt };
  }
  const m = Math.ceil(((n - 1) * dt) / d) + 1, g = new Array(m);
  let qi = at(inflow, 0), prevIn = qi, prevOut = qi;
  g[0] = qi;
  for (let k = 1; k < m; k++) {
    qi = at(inflow, (k * d) / dt);
    prevOut = Math.max(0, co.c0 * qi + co.c1 * prevIn + co.c2 * prevOut); prevIn = qi; g[k] = prevOut;
  }
  for (let j = 0; j < n; j++) out[j] = at(g, (j * dt) / d);
  return { out, lag: false, step: d };
}

const len = Math.hypot(b.x - a.x, b.y - a.y);
if (!(f.S > 0) || !(len > 0)) { problems.push({ id: p.id, kind: 'pipe', msg: t('rt.flat', { name: p.name }) }); continue; }
const full = manningFull(D, f.S, h.n, units), yn = normalDepth(f.Q, D, f.S, h.n, units);
const A = yn === null ? full.A : partFull(yn, D).A, V = A > 0 ? f.Q / A : full.V, Tw = yn === null ? D : Math.max(partFull(yn, D).T, 0.1 * D);
const c = CELERITY_RATIO * Math.max(V, 0.1), K = len / c / 60;      // minutes
const X = Math.min(0.5, Math.max(0, 0.5 * (1 - f.Q / (Tw * f.S * c * len))));
reach.set(p.id, { len, S: f.S, cap: full.Q, K, X, c, V });

const pipes = new Map();
for (const p of ord) {
  const r = reach.get(p.id), inflow = into.get(p.from);
  if (!r) continue;
  const { out, lag } = routeReach(inflow, r.K, r.X, dt);
  const a = peakOf(inflow, dt), b = peakOf(out, dt), dn = into.get(p.to);
  for (let j = 0; j < steps; j++) dn[j] += out[j];
  const over = b.peak > r.cap * (1 + 1e-9);
  pipes.set(p.id, { peakIn: a.peak, peakOut: b.peak, tPeakOut: b.tPeak, lag, att: a.peak > 0 ? Math.max(0, 1 - b.peak / a.peak) : 0, cap: r.cap, ratio: r.cap > 0 ? b.peak / r.cap : 0, K: r.K, X: r.X, c: r.c, V: r.V, over });
  if (over) problems.push({ id: p.id, kind: 'pipe', msg: t('rt.over', { name: p.name }) });
}
const structs = new Map(); for (const [id, q] of into) structs.set(id, { ...peakOf(q, dt), q });

/** The outlet's hydrograph (or any structure's) as CSV text: minutes, flow. */
export function hydrographCsv(r, id, units = 'ft') {
  const s = r.structs.get(id);
  if (!s) return '';
  return ['minutes,flow ' + (units === 'm' ? 'm3/s' : 'cfs'), ...r.t.map((m, j) => `${+m.toFixed(3)},${+s.q[j].toFixed(6)}`)].join('\r\n') + '\r\n';
}