Stack plume rise and air dispersion screening

/** Wind speed (m/s) at height z from the 10 m value (power law). */
export function windAt(u10, z, stab, terrain = 'rural') { return u10 * Math.pow(Math.max(z, 1) / 10, (P_EXP[terrain] || P_EXP.rural)[STABILITY.indexOf(stab)]); }

/** Effective release height. stack { hs, d, vs, Ts }; opts { u (at stack top), stab, Ta, x (m, for the transitional rise; omit for the final rise), dthetadz }.
 *  → { H, dh, hs, Fb, downwash, kind } */
export function plumeRise(stack, { u, stab = 'D', Ta = 293.15, x = null, dthetadz = null } = {}) {
  const { hs, d, vs, Ts } = stack;
  if (!(hs >= 0 && d > 0 && vs >= 0 && Ts > 0 && u > 0)) throw new Error('A stack needs height, diameter, exit velocity and temperature, and the wind must be above zero.');
  let h = hs, downwash = 0;
  if (vs < 1.5 * u) { downwash = 2 * d * (1.5 - vs / u); h = hs - downwash; }          // stack-tip downwash: h' = hs + 2d(vs/u - 1.5)
  const Fb = Ts > Ta ? G * vs * d * d * (Ts - Ta) / (4 * Ts) : 0;
  let dh, kind;
  if (Fb > 0) {
    if (stab === 'E' || stab === 'F') { const s = G / Ta * (dthetadz ?? (stab === 'E' ? 0.02 : 0.035)); dh = 2.6 * Math.cbrt(Fb / (u * s)); kind = 'buoyant stable'; }
    else { dh = Fb < 55 ? 21.425 * Math.pow(Fb, 0.75) / u : 38.71 * Math.pow(Fb, 0.6) / u; kind = 'buoyant'; }
    if (x != null && x > 0) dh = Math.min(dh, 1.6 * Math.cbrt(Fb) * Math.pow(x, 2 / 3) / u);
  } else { dh = 3 * d * vs / u; kind = 'momentum'; }
  return { H: Math.max(h, 0) + dh, dh, hs, Fb, downwash, kind };
}

/** A stack from a flowsheet stream (CH.2): `th` the Thermo of the flowsheet, `stream` {T K, P Pa, F mol/s, z[]}; species = ids to report (g/s each).
 *  Gas volume flow by the ideal gas law (the stack is near atmospheric pressure); exit velocity from the diameter. */
export function stackFromStream(th, stream, { hs, d, species = [], id = 'stack', x = 0, y = 0 } = {}) {
  if (!(hs >= 0 && d > 0)) throw new Error('Give the stack height and inner diameter.');
  const vol = stream.F * R * stream.T / stream.P, vs = vol / (Math.PI * d * d / 4), Q = {};
  for (const sp of species) { const i = th.ids.indexOf(sp); if (i < 0) throw new Error('Component ' + sp + ' is not in the flowsheet.'); Q[sp] = stream.F * stream.z[i] * th.c[i].MW; }   // MW g/mol -> g/s
  return { id, x, y, hs, d, vs, Ts: stream.T, vol, Q, from: 'stream' };
}

function atDistance(source, species, ctx, x, y = 0, z = 0) {
  const { u10, stab, terrain, Ta, zi } = ctx;
  const u0 = windAt(u10, source.hs, stab, terrain), r = plumeRise(source, { u: u0, stab, Ta, x });
  const u = windAt(u10, Math.max(r.H, 1), stab, terrain);
  return { c: plumeConc({ Q: (source.Q[species] ?? 0), u, H: r.H, stab, terrain, zi }, x, y, z) * 1e6, H: r.H, u };   // g/m3 -> ug/m3
}

/** Screening: worst 1-hour ground-level concentration (ug/m3) over the class x wind grid and distances. opts { species, terrain, Ta, zi, xMin, xMax, averaging: [3, 8, 24, 'annual'] }.
 *  → { max: { c, x, stab, u10, H }, byClass: { A: { c, x, u10 } … }, averaged: { 3: c … }, notes } */
export function screen(source, { species, terrain = 'rural', Ta = 293.15, zi = 0, xMin = 20, xMax = 30000, steps = 160, averaging = [3, 8, 24, 'annual'], wind = GRID } = {}) {
  if (!species || source.Q?.[species] == null) throw new Error('Name the species and give its rate (g/s) on the source.');
  const best = { c: 0 }, byClass = {};
  for (const stab of STABILITY) for (const u10 of wind[stab] || []) {
    for (let i = 0; i <= steps; i++) {
      const x = xMin * Math.pow(xMax / xMin, i / steps), r = atDistance(source, species, { u10, stab, terrain, Ta, zi }, x);
      if (!byClass[stab] || r.c > byClass[stab].c) byClass[stab] = { c: r.c, x, u10, H: r.H };
      if (r.c > best.c) Object.assign(best, { c: r.c, x, stab, u10, H: r.H });
    }
  }
  const averaged = {}; for (const a of averaging) averaged[a] = best.c * (AVERAGING[a] ?? NaN);
  const notes = ['Stability-class screening on flat ground, one source, no buildings. Ta ' + Ta + ' K' + (Ta === 293.15 ? ' (assumed: give the real ambient temperature)' : '') + '.',
    'Averaging multipliers are screening guidance values typed without the document open; the 1-hour maximum is the model\'s number.'];
  return { max: best, byClass, averaged, notes, species };
}

/** One source's contribution (ug/m3) at a receptor for one met hour. met { u10, dirFrom (deg clockwise from north, where the wind blows FROM), stab, zi?, Ta? } */
export function receptorHour(source, species, rec, met, terrain = 'rural') {
  const th = toRad((met.dirFrom + 180) % 360), dx = rec.x - source.x, dy = rec.y - source.y;
  const xd = dx * Math.sin(th) + dy * Math.cos(th), yc = dx * Math.cos(th) - dy * Math.sin(th);
  if (xd <= 1) return 0;
  return atDistance(source, species, { u10: met.u10, stab: met.stab, terrain, Ta: met.Ta ?? 293.15, zi: met.zi ?? 0 }, xd, yc, rec.z ?? 0).c;
}

/** Hour-by-hour run. sources [{ id, x, y, hs, d, vs, Ts, Q }], receptors [{ id, x, y, z? }], hours [met…]. opts { species, terrain, calm: 1 (m/s at 10 m, hours below it are skipped) }.
 *  → { receptors: { id: { max1h, second1h (high second high), mean, hourOfMax, rolling24 (max 24-h mean) } }, hours, calms, notes } */
export function hourly(sources, receptors, hours, { species, terrain = 'rural', calm = 1 } = {}) {
  if (!species) throw new Error('Name the species.');
  const out = {}; let calms = 0;
  const series = Object.fromEntries(receptors.map((r) => [r.id, new Array(hours.length).fill(0)]));
  hours.forEach((m, h) => {
    if (!(m.u10 >= calm)) { calms++; return; }
    for (const r of receptors) { let c = 0; for (const s of sources) c += receptorHour(s, species, r, m, terrain); series[r.id][h] = c; }
  });
  for (const r of receptors) {
    const s = series[r.id], sorted = s.map((c, i) => [c, i]).sort((a, b) => b[0] - a[0]);
    let roll = 0; if (s.length >= 24) { let acc = s.slice(0, 24).reduce((a, b) => a + b, 0); roll = acc; for (let i = 24; i < s.length; i++) { acc += s[i] - s[i - 24]; if (acc > roll) roll = acc; } roll /= 24; }
    out[r.id] = { max1h: sorted[0]?.[0] ?? 0, second1h: sorted[1]?.[0] ?? 0, hourOfMax: sorted[0]?.[1] ?? -1, mean: s.reduce((a, b) => a + b, 0) / Math.max(1, s.length), rolling24: roll };
  }
  return { receptors: out, hours: hours.length, calms, series, notes: ['Hour-by-hour stability-class Gaussian model; ' + calms + ' calm hour(s) (10 m wind below ' + calm + ' m/s) are skipped and count as zero in the means.', 'Not AERMOD: no boundary-layer scaling, no terrain, no building downwash.'] };
}

/** The Cend contour(s) of one source as plan shapes in site coordinates, for the wind blowing FROM dirFrom (deg clockwise from north). levels in ug/m3. */
export function toPlan(source, { species, levels, u10, stab, dirFrom, terrain = 'rural', Ta = 293.15, zi = 0 }) {
  const u0 = windAt(u10, source.hs, stab, terrain), r = plumeRise(source, { u: u0, stab, Ta }), u = windAt(u10, Math.max(r.H, 1), stab, terrain);
  const toAz = (dirFrom + 180) % 360, toMath = 90 - toAz, shapes = [];
  for (const lv of levels) {
    const f = plumeFootprint({ Q: source.Q[species], u, H: r.H, stab, terrain, zi }, lv * 1e-6, { cx: source.x, cy: source.y, toDeg: toMath });
    if (f.reach) shapes.push({ id: source.id + ':' + species + ':' + lv, type: 'plume', level: lv, unit: 'ug/m3', label: lv + ' ug/m3', polygon: f.polygon, reach: f.reach, H: r.H });
  }
  return { shapes, H: r.H, notes: ['One met condition, ground-level contour, flat ground.'] };
}

/** Compare results { pollutant, avg, c (ug/m3), where? }[] with a pack. Problems carry the pack and clause; an unconfirmed pack says so (rule 29). */
export function checkLimits(results, pack) {
  const pk = cleanPack(pack), ps = [];
  for (const r of results) {
    const lim = pk.limits.find((l) => l.pollutant === r.pollutant && l.avg === r.avg);
    if (!lim) { ps.push(problem('nolimit', 'note', `No limit for ${r.pollutant} (${r.avg}) in pack "${pk.name}": nothing was compared.`, [{ kind: 'pollutant', id: r.pollutant }])); continue; }
    const frac = r.c / lim.value, tag = pk.review.state === 'confirmed' ? '' : ' The pack is NOT confirmed by a person.';
    if (frac > 1) ps.push(problem('exceeds', 'error', `${r.pollutant} ${r.avg}: ${r.c.toPrecision(3)} ug/m3 is above ${lim.value} (${pk.name}, ${lim.clause || 'no clause'}).${tag}`, [{ kind: 'pollutant', id: r.pollutant }], { req: { op: '<=', value: lim.value, unit: 'ug/m3' }, found: r.c, pack: pk.name, clause: lim.clause }));
    else if (frac > 0.8) ps.push(problem('near', 'warning', `${r.pollutant} ${r.avg}: ${(frac * 100).toFixed(0)} % of the limit (${pk.name}).${tag}`, [{ kind: 'pollutant', id: r.pollutant }], { found: r.c, pack: pk.name, clause: lim.clause }));
  }
  return ps;
}

/**
 * plumedraw.js — the plume contours as drawing entities (0.8 Track CH, CH.19; DOM-free). A host (CadCraft's P&ID add-on, CivilCraft's site plan) draws them.
 *   plumeDrawing({ shapes (dispersion.toPlan(…).shapes), stack: { x, y, label }, windFrom (deg clockwise from north), receptors?, scale = 1, layer = 'AIR' })
 *   → { entities: [{ type: 'polyline' | 'line' | 'circle' | 'text', layer, … }], legend: [{ layer, label }], extent }
 * Units are the site coordinates of the shapes times `scale` (use 1000 for a drawing in mm). It adds the stack marker, a wind arrow that points the way the wind
 * blows, a label at the far end of each contour (level, reach) and a note that this is a stability-class screening, not AERMOD.
 */
export function plumeDrawing({ shapes, stack, windFrom, receptors = [], scale = 1, layer = 'AIR' }) {
  const E = [], legend = [], toAz = ((windFrom + 180) % 360) * Math.PI / 180, dx = Math.sin(toAz), dy = Math.cos(toAz);
  let minx = Infinity, miny = Infinity, maxx = -Infinity, maxy = -Infinity; const ext = (x, y) => { minx = Math.min(minx, x); miny = Math.min(miny, y); maxx = Math.max(maxx, x); maxy = Math.max(maxy, y); };
  for (const s of shapes) {
    const L = `${layer}-PLUME-${String(s.level).replace(/[^0-9A-Za-z.]/g, '_')}`; legend.push({ layer: L, label: s.label });
    E.push({ type: 'polyline', layer: L, closed: true, points: s.polygon.map((p) => [p[0] * scale, p[1] * scale]), id: s.id });
    for (const p of s.polygon) ext(p[0] * scale, p[1] * scale);
    const far = s.polygon.reduce((a, p) => ((p[0] - stack.x) * dx + (p[1] - stack.y) * dy > (a[0] - stack.x) * dx + (a[1] - stack.y) * dy ? p : a), s.polygon[0]);
    E.push({ type: 'text', layer: layer + '-TEXT', x: far[0] * scale, y: far[1] * scale, text: `${s.label}, reach ${s.reach.toFixed(0)} m`, h: 3 * scale });
  }
  const sx = stack.x * scale, sy = stack.y * scale, len = 30 * scale; ext(sx, sy);
  E.push({ type: 'circle', layer: layer + '-STACK', cx: sx, cy: sy, r: 1.5 * scale }, { type: 'text', layer: layer + '-TEXT', x: sx + 2 * scale, y: sy - 4 * scale, text: stack.label || 'Stack', h: 3 * scale });
  const ax = sx + dx * len, ay = sy + dy * len, hx = 0.25 * len;
  E.push({ type: 'line', layer: layer + '-WIND', x1: sx, y1: sy, x2: ax, y2: ay }, { type: 'line', layer: layer + '-WIND', x1: ax, y1: ay, x2: ax - hx * (dx * 0.87 - dy * 0.5), y2: ay - hx * (dy * 0.87 + dx * 0.5) }, { type: 'line', layer: layer + '-WIND', x1: ax, y1: ay, x2: ax - hx * (dx * 0.87 + dy * 0.5), y2: ay - hx * (dy * 0.87 - dx * 0.5) });
  E.push({ type: 'text', layer: layer + '-TEXT', x: ax, y: ay, text: `Wind from ${windFrom.toFixed(0)} deg`, h: 3 * scale });
  for (const r of receptors) { E.push({ type: 'circle', layer: layer + '-RECEPTOR', cx: r.x * scale, cy: r.y * scale, r: 1 * scale }, { type: 'text', layer: layer + '-TEXT', x: r.x * scale + 1.5 * scale, y: r.y * scale, text: `${r.id}${r.c != null ? ' ' + r.c.toPrecision(3) + ' ug/m3' : ''}`, h: 2.5 * scale }); ext(r.x * scale, r.y * scale); }
  E.push({ type: 'text', layer: layer + '-TEXT', x: minx, y: miny - 6 * scale, text: 'Stability-class screening, not AERMOD. Not an engineering approval.', h: 2.5 * scale });
  return { entities: E, legend, extent: { minx, miny, maxx, maxy } };
}