Release, fire, explosion and toxic screening

/** Liquid through a hole. dP gauge Pa over the outside, h liquid head above the hole (m). → { mdot kg/s, v m/s } */
export function liquidRelease({ A, d, Cd = 0.61, rho, dP = 0, h = 0 }) {
  const a = A ?? (d ? Math.PI * d * d / 4 : 0); if (!(a > 0) || !(rho > 0)) throw new Error('Liquid release needs a hole area (or diameter) and a density.');
  const v = Math.sqrt(Math.max(0, 2 * dP / rho + 2 * G * h)); return { mdot: Cd * a * rho * v, v, A: a };
}

/** Gas through a hole, P0 and Pa absolute (Pa). Cd defaults to 1 (the conservative value; a person may lower it). */
export function gasRelease({ A, d, Cd = 1, P0, Pa = P_ATM, T, M, k = 1.4, Z = 1 }) {
  const a = A ?? (d ? Math.PI * d * d / 4 : 0); if (!(a > 0) || !(P0 > 0) || !(T > 0) || !(M > 0) || !(k > 1)) throw new Error('Gas release needs area (or diameter), P0, T, M and k > 1.');
  if (P0 <= Pa) return { mdot: 0, choked: false, ratio: Pa / P0 };
  const crit = Math.pow((k + 1) / 2, k / (k - 1)), choked = P0 / Pa >= crit;
  if (choked) return { mdot: Cd * a * P0 * Math.sqrt(k * M / (Z * R * T) * Math.pow(2 / (k + 1), (k + 1) / (k - 1))), choked, ratio: Pa / P0 };
  const r = Pa / P0; return { mdot: Cd * a * P0 * Math.sqrt(2 * M / (Z * R * T) * k / (k - 1) * (Math.pow(r, 2 / k) - Math.pow(r, (k + 1) / k))), choked, ratio: r };
}

/** Fraction of a pressurised liquid that flashes to vapour on release to atmosphere (0 if T <= Tb). cp J/kg K, hfg J/kg. */
export function flashFraction({ cp, T, Tb, hfg }) { if (!(cp > 0 && hfg > 0)) throw new Error('Flash needs cp and the heat of vaporisation.'); return T <= Tb ? 0 : 1 - Math.exp(-cp * (T - Tb) / hfg); }

/** Evaporation from a pool, kg/s. U wind m/s at about 10 m, Pv Pa, Sc the vapour's Schmidt number in air (the person's). */
export function poolEvaporation({ area, U, Pv, M, T, Sc }) {
  if (!(area > 0 && U >= 0 && Pv >= 0 && M > 0 && T > 0 && Sc > 0)) throw new Error('Evaporation needs area, wind, vapour pressure, M, T and Sc.');
  const d = Math.sqrt(4 * area / Math.PI), kg = 0.004786 * Math.pow(Math.max(U, 0.1), 0.78) * Math.pow(d, -0.11) * Math.pow(Sc, -0.67);
  return { mdot: kg * area * M * Pv / (R * T), kg, d };
}

/** Radiation distance (m) from the point source to the flux level q (W/m²): r = sqrt(tau chi Q/(4 pi q)). chi (radiative fraction) and tau (transmissivity) are the person's. */
export function radiationDistance({ Q, chi, tau = 1, q }) { if (!(Q > 0 && chi > 0 && q > 0)) throw new Error('Radiation distance needs Q, the radiative fraction chi and the flux level q.'); return Math.sqrt(tau * chi * Q / (4 * Math.PI * q)); }

/** Pool fire: D pool diameter (m), mPrime burning flux kg/(m² s) (the person's, from the fuel), flame height by Thomas, point source at mid-height; ground distances to each level. */
export function poolFire({ D, mPrime, dHc, chi, tau = 1, levels, rhoAir = 1.2 }) {
  if (!(D > 0 && mPrime > 0)) throw new Error('A pool fire needs a diameter and a burning flux.');
  const H = 42 * D * Math.pow(mPrime / (rhoAir * Math.sqrt(G * D)), 0.61), area = Math.PI * D * D / 4, Q = mPrime * area * dHc;
  return { H, Q, mdot: mPrime * area, distances: levels.map((q) => { const r = radiationDistance({ Q, chi, tau, q }), h = H / 2; return { q, r, ground: r > h ? Math.sqrt(r * r - h * h) : 0 }; }) };
}

/** TNT-equivalent mass (kg). eta is the yield fraction the person justifies; no default is given. */
export function tntEquivalent({ m, dHc, eta, eTNT = 4.184e6 }) { if (!(m > 0 && dHc > 0 && eta > 0 && eta <= 1)) throw new Error('TNT equivalence needs mass, dHc and an explicit yield eta between 0 and 1.'); return eta * m * dHc / eTNT; }

/** Side-on overpressure (Pa) at R metres from W kg TNT (Kinney-Graham). */
export function blastOverpressure(W, Rm, Pa = P_ATM) { if (!(W > 0 && Rm > 0)) return 0; const Z = Rm / Math.cbrt(W); return Pa * 808 * (1 + (Z / 4.5) ** 2) / Math.sqrt((1 + (Z / 0.048) ** 2) * (1 + (Z / 0.32) ** 2) * (1 + (Z / 1.35) ** 2)); }

/** Distance (m) at which the overpressure falls to P (Pa); monotone, found by bisection. */
export function blastDistance(W, P, Pa = P_ATM) { if (!(W > 0 && P > 0)) return 0; let lo = 0.01, hi = 1e5; if (blastOverpressure(W, lo, Pa) < P) return 0; for (let i = 0; i < 80; i++) { const m = Math.sqrt(lo * hi); if (blastOverpressure(W, m, Pa) >= P) lo = m; else hi = m; } return lo; }

// ---- toxic --------------------------------------------------------------------------------------------------------------------------
export const probit = ({ a, b, n, C, t }) => a + b * Math.log(Math.pow(C, n) * t);

const erf = (x) => { const t = 1 / (1 + 0.3275911 * Math.abs(x)), y = 1 - (((((1.061405429 * t - 1.453152027) * t) + 1.421413741) * t - 0.284496736) * t + 0.254829592) * t * Math.exp(-x * x); return x >= 0 ? y : -y; };

export const probitToProbability = (pr) => 0.5 * (1 + erf((pr - 5) / Math.SQRT2));

function interpLog(pts, R, col) {
  if (R <= pts[0][0]) return pts[0][col]; if (R >= pts[pts.length - 1][0]) return pts[pts.length - 1][col];
  for (let i = 0; i + 1 < pts.length; i++) if (R >= pts[i][0] && R <= pts[i + 1][0]) return Math.exp(lerp(Math.log(R), Math.log(pts[i][0]), Math.log(pts[i + 1][0]), Math.log(pts[i][col]), Math.log(pts[i + 1][col])));
  return NaN;
}

export const sachsLength = (E, P0 = 101325) => Math.cbrt(E / P0);

export function overpressure(pack, strength, E, R, P0 = 101325) {
  const c = cleanPack(pack); if (!c.ok) { const e = new Error(NEEDS_SOURCE + ' (' + c.problems.join('; ') + ')'); e.needsSource = true; throw e; }
  const cv = c.pack.curves[strength]; if (!cv) throw new Error('the pack has no curve for strength ' + strength);
  const Rb = R / sachsLength(E, P0), Pb = interpLog(cv.points, Rb, 1), td = cv.points[0].length >= 3 ? interpLog(cv.points.filter((q) => q.length >= 3), Rb, 2) : null;
  return { dP: Pb * P0, Rbar: Rb, Pbar: Pb, duration: td != null ? td * sachsLength(E, P0) / 340 : null, reviewed: c.pack.review.state === 'confirmed' };
}

/** Distance at which the overpressure falls to dP (Pa); the curve must be monotone decreasing in the range. */
export function distanceFor(pack, strength, E, dP, P0 = 101325) {
  let lo = 1e-3 * sachsLength(E, P0), hi = 1e3 * sachsLength(E, P0);
  if (overpressure(pack, strength, E, lo, P0).dP < dP) return lo; if (overpressure(pack, strength, E, hi, P0).dP > dP) return hi;
  for (let i = 0; i < 100; i++) { const m = Math.sqrt(lo * hi); if (overpressure(pack, strength, E, m, P0).dP > dP) lo = m; else hi = m; }
  return Math.sqrt(lo * hi);
}

/** sources [{ id, x, y, volume (m3 of flammable mixture), energyDensity (J/m3, the person's), strength 1..10 }]; at {x, y}. The method treats separate sources separately: the largest governs. */
export function cloud(pack, sources, at, P0 = 101325) {
  const rows = sources.map((s) => { const E = s.volume * s.energyDensity, R = Math.max(Math.hypot(at.x - s.x, at.y - s.y), 0.01); return { id: s.id, E, R, ...overpressure(pack, s.strength, E, R, P0) }; });
  const g = rows.reduce((a, b) => (b.dP > a.dP ? b : a), rows[0]);
  return { rows, governing: g.id, dP: g.dP, notes: ['sources are treated separately (no interaction); this is a screening estimate, not an engineering approval'] };
}

/** Contour circles (plan shapes in site coordinates) of the governing source at given overpressures (Pa), for CivilCraft's site plan. */
export function contours(pack, source, levels, P0 = 101325) {
  const E = source.volume * source.energyDensity;
  return levels.map((l) => ({ id: `${source.id}:blast:${l}`, type: 'circle', cx: source.x, cy: source.y, r: distanceFor(pack, source.strength, E, l, P0), label: `${(l / 1000).toFixed(1)} kPa`, level: l, unit: 'Pa' }));
}

export function instantaneous(spec) {
  const { mass, rhoGas, rhoAir = 1.2, Cg, aTop, aEdge, ustar = 0, uWind = 0, R0, tEnd, cEnd = null } = spec;
  for (const [k, v] of Object.entries({ mass, rhoGas, Cg, R0, tEnd })) if (!(v > 0)) throw new Error('heavy gas: give ' + k);
  if (!(aTop >= 0 && aEdge >= 0)) throw new Error('heavy gas: give the entrainment constants aTop and aEdge (zero is allowed)');
  if (!(rhoGas > rhoAir)) throw new Error('the gas is not heavier than air at this temperature: use the plume models instead');
  const Vg = mass / rhoGas; let V = Vg / 1, R = R0, t = 0, x = 0; V = Math.max(V, 1e-9); let H = V / (Math.PI * R * R);
  const out = { t: [0], R: [R], H: [H], c: [Vg / V], x: [0], reach: null }, N = 40000; let dt = tEnd / N;
  const state = (R_, V_) => { const H_ = V_ / (Math.PI * R_ * R_), c = Vg / V_, rho = rhoAir + c * (rhoGas - rhoAir), gp = G * (rho - rhoAir) / rhoAir; return { H_, c, gp }; };
  for (let s = 0; s < N && t < tEnd; s++) {
    const f = (R_, V_) => { const { H_, gp } = state(R_, V_), dR = Cg * Math.sqrt(Math.max(gp, 0) * H_), dV = Math.PI * R_ * R_ * aTop * ustar + 2 * Math.PI * R_ * H_ * aEdge * dR; return [dR, dV]; };
    const [a1, b1] = f(R, V), [a2, b2] = f(R + dt / 2 * a1, V + dt / 2 * b1), [a3, b3] = f(R + dt / 2 * a2, V + dt / 2 * b2), [a4, b4] = f(R + dt * a3, V + dt * b3);
    R += dt / 6 * (a1 + 2 * a2 + 2 * a3 + a4); V += dt / 6 * (b1 + 2 * b2 + 2 * b3 + b4); t += dt; x += uWind * dt;
    const st = state(R, V);
    if (s % 40 === 0) { out.t.push(t); out.R.push(R); out.H.push(st.H_); out.c.push(st.c); out.x.push(x); }
    if (cEnd != null && st.c <= cEnd && !out.reach) { out.reach = { time: t, R, x }; break; }
  }
  const st = state(R, V); out.final = { R, H: st.H_, c: st.c, V, gprime: st.gp }; out.massConserved = out.final.c * V * rhoGas; out.notes = ['isothermal, flat ground, no obstacles, no heat transfer from the ground; once g\' is small the cloud is passive and the plume models take over'];
  return out;
}