Water supply and duct network sizing

export function engineRuns(pipes) {
  const runs = [], attrs = {}, ends = [], joints = [];
  for (const { id, data: d } of pipes) {
    if (!d?.pts || d.pts.length < 2 || !d.segs?.length) continue;
    const m = mmDraw(d), z0 = (Number(d.z) || 0) * m;
    const cum = [0]; for (let i = 1; i < d.pts.length; i++) cum.push(cum[i - 1] + Math.hypot(d.pts[i][0] - d.pts[i - 1][0], d.pts[i][1] - d.pts[i - 1][1]));
    const slope = Number(d.slope) || 0;   // a ratio: the run falls slope × (length along it), as pipe.js zAt
    const zAt = (i) => z0 - slope * cum[i] * m;
    let g = 0, from = 0;
    for (let i = 1; i <= d.segs.length; i++) {
      if (i < d.segs.length && sameSize(d.segs[i], d.segs[from])) continue;
      const verts = d.pts.slice(from, i + 1).map((p, j) => [p[0] * m, p[1] * m, zAt(from + j)]);
      const runId = `${id}.${g++}`, size = d.segs[from].size;
      const sizes = sizesFor(d.mat, d.cls).map((s) => pipeSize(d.mat, d.cls, s)).filter(Boolean).map((p) => ({ key: p.size, idMm: p.id * MM_IN }));
      const cur = pipeSize(d.mat, d.cls, size);
      runs.push({ id: runId, verts, valves: [] });
      attrs[runId] = { mat: d.mat, idMm: cur ? cur.id * MM_IN : 0, sizes, label: size };
      if (from === 0) ends.push({ run: runId, p: verts[0], tag: d.tag || String(id) });
      if (i === d.segs.length) ends.push({ run: runId, p: verts[verts.length - 1], tag: d.tag || String(id) });
      from = i;
    }
  }
  return { runs, attrs, ends, joints };
}

let near = null;
for (const e of ends) { const d = Math.hypot(e.p[0] - source.x * m, e.p[1] - source.y * m); if (!near || d < near.d) near = { d, e }; }
if (!near || near.d > tolMm) return { ok: false, why: 'source' };

const rows = result.segs.map((s) => ({ run: s.run, lengthM: s.length, fu: s.fu, lps: s.flow * 1000, size: s.size, V: s.V, rate: s.rate, dpKpa: s.dp / 1000, status: s.status, beyond: s.beyond }));

export const csvRows = (a) => ['run,length_m,fixture_units,flow_lps,size,velocity_ms,rate_pa_m,loss_kpa,status', ...a.rows.map((r) => [r.run, r.lengthM.toFixed(2), r.fu, r.lps.toFixed(3), r.size, r.V.toFixed(2), r.rate.toFixed(1), r.dpKpa.toFixed(2), r.status].join(','))].join('\n') + '\n';

export const G = 9.80665;
export const WATER = { rho: 998.2, nu: 1.004e-6 };            // 20 °C — a labelled default
export const GPM = 0.003785411784 / 60;                        // m³/s per US gpm
export const KPA_PSI = 6.894757;

/** Hazen–Williams friction loss, Pa/m, for Q (m³/s) in a pipe of inside diameter D (m). */
export function hazenWilliams(Q, D, C = 130) {
  if (!(Q > 0) || !(D > 0)) return 0;
  return (10.67 * Q ** 1.852 / (C ** 1.852 * D ** 4.871)) * WATER.rho * G;
}

/** Darcy–Weisbach friction loss for water, Pa/m (eps in m). */
export function darcyWater(Q, D, eps = 0.0000015) {
  if (!(Q > 0) || !(D > 0)) return 0;
  const V = Q / (Math.PI * D * D / 4), Re = V * D / WATER.nu;
  return (colebrook(Re, eps / D) / D) * WATER.rho * V * V / 2;
}

export const velocity = (Q, D) => (D > 0 ? Q / (Math.PI * D * D / 4) : 0);

/** Demand (m³/s) for `fu` fixture units: log–log interpolation; below the first point in proportion; above the last by the last slope (`beyond` = true). */
export function demandFor(fu, curve = DEFAULT_CURVE) {
  const c = cleanCurve(curve) || DEFAULT_CURVE, k = c.unit === 'lps' ? 0.001 : GPM, P = c.points;
  if (!(fu > 0)) return { q: 0, beyond: false };
  if (fu <= P[0][0]) return { q: P[0][1] * k * (fu / P[0][0]), beyond: false };
  let i = P.findIndex((p) => p[0] >= fu);
  const beyond = i < 0;
  if (beyond) i = P.length - 1;
  const [x0, y0] = P[i - 1], [x1, y1] = P[i];
  const s = Math.log(y1 / y0) / Math.log(x1 / x0);
  return { q: y0 * (fu / x0) ** s * k, beyond };
}

/** Pressure (Pa) lost in the fittings of one segment: its bends, valves and the junction it leaves (a tee branch when it is not the biggest child). */
function fittingLoss(seg, V, rho, K, branchOf, round) {
  let k = 0;
  for (const a of seg.bends || []) if (a > 15) k += bendK(round ? (K.elbow90_round ?? K.elbow90) : K.elbow90, a);
  for (const v of seg.valves || []) k += K[v] ?? 0;
  if (branchOf) k += branchOf === 'branch' ? K.tee_branch : K.tee_run;
  return k * rho * V * V / 2;
}

const Q = flow.get(id), C = o.C?.[a.mat] ?? HW_C[a.mat] ?? HW_C.default, D = (a.idMm || 0) / 1000;
const V = velocity(Q, D), rate = hazenWilliams(Q, D, C);
let need = null;
for (const z of a.sizes || []) { const d = z.idMm / 1000; if (velocity(Q, d) <= o.vmax + 1e-9 && (!o.rateMax || hazenWilliams(Q, d, C) / 1000 <= o.rateMax + 1e-12)) { need = z; break; } }
const fit = D > 0 ? fittingLoss(s, V, WATER.rho, o.K, jn.get(id) || null, true) : 0;
const dp = rate * s.length + fit + WATER.rho * G * s.dz;     // friction + fittings + static (up is a cost)

const path = pathTo(tree, l.node);
if (!path.length && l.node !== source) { out.push({ ...l, reached: false }); continue; }
let dp = 0; for (const id of path) dp += segs.get(id).dp;
const left = o.availKpa != null ? o.availKpa - dp / 1000 : null;
out.push({ ...l, reached: true, path, dpKpa: dp / 1000, leftKpa: left, short: o.minKpa != null && left != null ? left < o.minKpa : false });

const w = (a.wMm || 0) / 1000, h = (a.shape === 'round' ? a.wMm : a.hMm || 0) / 1000;
let V = 0, rate = 0, de = 0, area = 0;
if (w > 0 && h > 0 && Q > 0) { const L = sectionLoss({ q: Q * 1000, shape: a.shape || 'rect', w: w * 1000, h: h * 1000, u: 'mm', mat }); V = L.V; rate = L.rate; de = L.de / 1000; area = section(a.shape || 'rect', w, h).area; }
const fit = V > 0 ? fittingLoss(s, V, AIR.rho, { ...o.C, elbow90: o.C.elbow90, tee_run: o.C.tee_run, tee_branch: o.C.tee_branch }, jn.get(id) || null, a.shape === 'round') : 0;
// needed size: the equal-friction round diameter, or the diameter at the velocity limit — in round equivalents
let needDe = null;
if (Q > 0) needDe = o.method === 'velocity' ? Math.sqrt(4 * Q / (Math.PI * o.vmax)) : diameterAtRate(Q, o.rate, ROUGHNESS[mat] ?? ROUGHNESS.galv);
const aspect = a.shape === 'round' || !(w && h) ? 1 : Math.max(w, h) / Math.min(w, h);
const dp = rate * s.length + fit + 0;     // ducts: static head cancels in air at building heights (stack effect is ignored, said in the help)

const crit = res.filter((r) => r.reached).reduce((m, r) => (!m || r.dp > m.dp ? r : m), null);
for (const r of res) if (r.reached && crit) r.damperPa = crit.dp - r.dp;    // pressure the balancing damper must take up
return { tree, segs: [...segs.values()], terminals: res, critical: crit, fanPa: crit ? crit.dp + (o.equipDp || 0) : 0, cycles: tree.cycles, unreached: tree.unreached, opts: o };

/** The round diameter (m) giving `rate` Pa/m at Q — sizing.js bisection. */
function diameterAtRate(Q, rate, eps) {
  let lo = Math.log(0.01), hi = Math.log(5);
  if (roundLoss(Q, 5, eps).rate > rate) return null;
  for (let i = 0; i < 100; i++) { const m = (lo + hi) / 2; if (roundLoss(Q, Math.exp(m), eps).rate > rate) lo = m; else hi = m; }
  return Math.exp(hi);
}

/** The foot of p on the polyline: { s (arc length mm), d (distance mm) }. */
function project(verts, p) {
  let best = { s: 0, d: Infinity }, acc = 0;
  for (let i = 0; i + 1 < verts.length; i++) {
    const a = verts[i], v = sub(verts[i + 1], a), L2 = dot(v, v), L = Math.sqrt(L2);
    let u = L2 > 0 ? dot(sub(p, a), v) / L2 : 0;
    u = Math.max(0, Math.min(1, u));
    const q = [a[0] + v[0] * u, a[1] + v[1] * u, (a[2] || 0) + v[2] * u], d = dist(p, q);
    if (d < best.d) best = { s: acc + L * u, d };
    acc += L;
  }
  return best;
}

/** The point at arc length s of the polyline. */
function at(verts, s) {
  let acc = 0;
  for (let i = 0; i + 1 < verts.length; i++) {
    const L = dist(verts[i], verts[i + 1]);
    if (s <= acc + L || i + 2 === verts.length) { const u = L > 0 ? Math.max(0, Math.min(1, (s - acc) / L)) : 0; const a = verts[i], b = verts[i + 1]; return [a[0] + (b[0] - a[0]) * u, a[1] + (b[1] - a[1]) * u, (a[2] || 0) + ((b[2] || 0) - (a[2] || 0)) * u]; }
    acc += L;
  }
  return verts[0].slice();
}

const turn = (a, b, c) => {
  const u = sub(b, a), v = sub(c, b), lu = Math.hypot(...u), lv = Math.hypot(...v);
  if (!lu || !lv) return 0;
  return Math.acos(Math.max(-1, Math.min(1, dot(u, v) / (lu * lv)))) * 180 / Math.PI;
};

const hit = (xi, p, tag) => { const pr = project(runInfo[xi].verts, p); if (pr.d <= tol) { cuts[xi].push({ s: pr.s, p: at(runInfo[xi].verts, pr.s), tag }); return true; } return false; };
// a run's end lying on another run's line cuts the other run (a tee)
runInfo.forEach((x, i) => { for (const e of [x.verts[0], x.verts[x.verts.length - 1]]) runInfo.forEach((y, j) => { if (i !== j) hit(j, e, null); }); });
for (const pt of points) runInfo.forEach((x, i) => hit(i, pt.p, pt.key));
// nodes: cluster every cut point within tol (a small union-find)
const all = [];
cuts.forEach((c, xi) => c.forEach((k) => all.push({ xi, k })));
const parent = all.map((_, i) => i), find = (i) => { while (parent[i] !== i) { parent[i] = parent[parent[i]]; i = parent[i]; } return i; };
for (let i = 0; i < all.length; i++) for (let j = i + 1; j < all.length; j++) if (dist(all[i].k.p, all[j].k.p) <= tol) parent[find(j)] = find(i);

for (let k = 0; k + 1 < uniq.length; k++) {
  const A = uniq[k], B = uniq[k + 1];
  if (A.node === B.node) continue;
  const verts = [A.k.p];
  const bends = [];
  for (const cr of corners) if (cr.s > A.k.s + tol * 0.5 && cr.s < B.k.s - tol * 0.5) { verts.push(cr.b); bends.push(turn(cr.a, cr.b, cr.c)); }
  verts.push(B.k.p);
  let len = 0; for (let i = 1; i < verts.length; i++) len += dist(verts[i - 1], verts[i]);
  const valves = (x.r.valves || []).filter((v) => v.s >= A.k.s && v.s < B.k.s).map((v) => v.kind);
  segs.push({ id: x.r.id + ':' + k, run: x.r.id, a: A.node, b: B.node, verts, length: len / 1000, dz: ((B.k.p[2] || 0) - (A.k.p[2] || 0)) / 1000, bends, valves });
}

export function accumulate(tree, loadAt) {
  const sum = new Map();
  for (let i = tree.order.length - 1; i >= 0; i--) {
    const s = tree.segOf.get(tree.order[i]);
    let v = loadAt.get(s.b) || 0;
    for (const c of tree.children.get(s.b) || []) v += sum.get(c) || 0;
    sum.set(s.id, v);
  }
  return sum;
}

export const AIR = { rho: 1.204, nu: 1.51e-5 };          // kg/m³, m²/s — ASHRAE standard air
export const ROUGHNESS = { smooth: 0.03e-3, galv: 0.09e-3, board: 0.9e-3, flex: 3.0e-3 };   // m
export const CFM = 0.028316846592 / 60;                  // m³/s per cfm (1 ft³ = 0.028316846592 m³)
export const LPS = 0.001;                                 // m³/s per L/s
export const IN_WG = 249.08891;                           // Pa per inch of water (conventional, 4 °C)
export const PER_100FT = IN_WG / 30.48;                   // Pa/m per in. w.g./100 ft  (≈ 8.1722)
export const FPM = 0.00508;                               // m/s per ft/min
export const IN = 0.0254, MM = 0.001;

/** Colebrook friction factor (iterated from the Swamee–Jain estimate; laminar 64/Re below Re 2300). */
export function colebrook(Re, rel) {
  if (!(Re > 0)) return NaN;
  if (Re < 2300) return 64 / Re;
  let f = 0.25 / Math.log10(rel / 3.7 + 5.74 / Re ** 0.9) ** 2;
  for (let i = 0; i < 60; i++) {
    const x = -2 * Math.log10(rel / 3.7 + 2.51 / (Re * Math.sqrt(f)));
    const g = 1 / (x * x);
    if (Math.abs(g - f) < 1e-14) { f = g; break; }
    f = g;
  }
  return f;
}

/** Altshul–Tsal (ASHRAE's explicit form) — a cross-check, not used for sizing. */
export function altshulTsal(Re, rel) {
  const f1 = 0.11 * (rel + 68 / Re) ** 0.25;
  return f1 >= 0.018 ? f1 : 0.85 * f1 + 0.0028;
}

/** Friction loss of a round duct: Pa/m for airflow Q (m³/s) in diameter D (m), roughness eps (m). Also V, Re, f. */
export function roundLoss(Q, D, eps = ROUGHNESS.galv) {
  const A = (Math.PI * D * D) / 4, V = Q / A, Re = reynolds(V, D), f = colebrook(Re, eps / D);
  return { rate: (f / D) * (AIR.rho * V * V) / 2, V, Re, f };
}

/** The round diameter (m) giving friction rate `rate` (Pa/m) at Q (m³/s): bisection on log D, 10 mm … 5 m. */
export function diameterFor(Q, rate, eps = ROUGHNESS.galv) {
  if (!(Q > 0) || !(rate > 0)) return NaN;
  let lo = Math.log(0.01), hi = Math.log(5);
  if (roundLoss(Q, 5, eps).rate > rate) return NaN;   // more than a 5 m duct carries at that rate
  for (let i = 0; i < 200; i++) {
    const mid = (lo + hi) / 2;
    if (roundLoss(Q, Math.exp(mid), eps).rate > rate) lo = mid; else hi = mid;
  }
  return Math.exp(hi);
}

/** Huebscher's circular equivalent of a rectangle a × b (any one length unit). */
export const deRect = (a, b) => (1.3 * (a * b) ** 0.625) / (a + b) ** 0.25;

export const ovalArea = (a, b) => (Math.PI * b * b) / 4 + b * (a - b);
export const ovalPerimeter = (a, b) => Math.PI * b + 2 * (a - b);
export const deOval = (a, b) => (a - b <= 1e-12 ? b : (1.55 * ovalArea(a, b) ** 0.625) / ovalPerimeter(a, b) ** 0.25);

/** Cross-section area and perimeter of a size (same length unit as the size): rect w×h, round d (= w), oval w×h. */
export function section(shape, w, h) {
  if (shape === 'round') return { area: (Math.PI * w * w) / 4, perimeter: Math.PI * w, de: w };
  if (shape === 'oval') return { area: ovalArea(w, h), perimeter: ovalPerimeter(w, h), de: deOval(w, h) };
  return { area: w * h, perimeter: 2 * (w + h), de: deRect(w, h) };
}

/**
 * Size a section by equal friction. q = airflow in cfm ('in') or L/s ('mm'); rate in in. w.g./100 ft or Pa/m; shape and,
 * for rect/oval, the fixed side `h` (height, or the oval's minor axis) in in or mm. → { shape, w, h, d (exact round
 * diameter), de, V (fpm or m/s), rate (actual, typed units), aspect, note } in the typed units; null when it cannot be sized.
 */
export function sizeSection({ q, rate, shape = 'round', h = null, u = 'in', mat = 'galv', maxAspect = 4 }) {
  const Q = q * (u === 'in' ? CFM : LPS), R = rate * (u === 'in' ? PER_100FT : 1), unit = u === 'in' ? IN : MM;
  const eps = ROUGHNESS[mat] ?? ROUGHNESS.galv;
  const D = diameterFor(Q, R, eps);
  if (!Number.isFinite(D)) return null;
  const dReq = D / unit;
  let w, hh = h, de;
  if (shape === 'round') {
    const list = u === 'in' ? ROUND_IN : ROUND_MM;
    w = list.find((x) => x >= dReq - 1e-9);
    if (!w) return null;
    hh = w; de = w;
  } else {
    const { step, min } = STEP[u];
    if (!(hh > 0)) return null;
    const deOf = shape === 'oval' ? (a) => deOval(a, hh) : (a) => deRect(a, hh);
    w = shape === 'oval' ? Math.max(hh, min) : min;
    while (deOf(w) < dReq - 1e-9) { w += step; if (w > 400 * step) return null; }
    de = deOf(w);
  }
  const s = section(shape, w, hh);
  const V = Q / (s.area * unit * unit);
  const actual = roundLoss(Q, de * unit, eps).rate / (u === 'in' ? PER_100FT : 1);
  const aspect = shape === 'round' ? 1 : Math.max(w, hh) / Math.min(w, hh);
  return { shape, w, h: hh, d: dReq, de, V: u === 'in' ? V / FPM : V, rate: actual, aspect, note: aspect > maxAspect ? 'aspect' : null };
}

/** Friction rate and velocity of a section as drawn (typed units): rate in in. w.g./100 ft or Pa/m, V in fpm or m/s. */
export function sectionLoss({ q, shape, w, h, u = 'in', mat = 'galv' }) {
  const Q = q * (u === 'in' ? CFM : LPS), unit = u === 'in' ? IN : MM;
  const s = section(shape, w, h);
  const r = roundLoss(Q, s.de * unit, ROUGHNESS[mat] ?? ROUGHNESS.galv);
  const V = Q / (s.area * unit * unit);
  return { rate: r.rate / (u === 'in' ? PER_100FT : 1), V: u === 'in' ? V / FPM : V, Re: r.Re, f: r.f, de: s.de };
}