Hazardous area extents and zone geometry

/** Dilution method: minimum ventilation flow (m³/s) and the hypothetical volume (m³) and the sphere radius that volume fills. */
export function dilution(s, pack = SCREENING_PACK, { Pa = 101325 } = {}) {
  const rr = releaseRate(s), v = s.ventilation; if (!rr || !(s.substance.lfl > 0) || !(s.substance.M > 0)) return null;
  const T = v.T || 293.15, lflm = s.substance.lfl * s.substance.M * Pa / (R * T), k = pack.k[s.grade] || 0.5, Qmin = rr.rate / (k * lflm) * (T / 293);
  const out = { rate: rr.rate, lflMass: lflm, k, Qmin };
  if (v.C > 0) { out.Vz = (v.f || 1) * Qmin / v.C; out.radius = Math.cbrt(3 * out.Vz / (4 * Math.PI)); }
  return out;
}

/** Classify every source. → { zones: [{ id, source, tag, zone, grade, extent (m), extentFrom, shape }], problems, notes } */
export function classify(sources, pack = SCREENING_PACK) {
  const zones = [], problems = [], notes = [DISCLAIMER, 'Classification pack: ' + pack.name + ' (review ' + pack.review.state + ').', REVIEW.note];
  for (const s of sources) {
    const about = [{ kind: 'release-source', id: s.id, label: s.tag }], z = zoneFor(pack, s);
    if (!z) { problems.push(problem('area.zone', 'warning', `Source ${s.id}: the pack gives no zone for ${s.grade} / ${s.ventilation.degree} / ${s.ventilation.availability}.`, about)); continue; }
    if (!(s.substance.lfl > 0)) problems.push(problem('area.lfl', 'warning', `Source ${s.id} has no lower flammable limit; the extent cannot be worked out.`, about));
    if (z === 'NE') { zones.push({ id: s.id, source: s.id, tag: s.tag, zone: 'NE', grade: s.grade, extent: 0, extentFrom: 'none', shape: null }); continue; }
    let r = 0, from = '';
    if (s.extent?.manual > 0) { r = s.extent.manual; from = 'given by the person'; }
    else { const d = dilution(s, pack); if (d && d.radius > 0) { r = d.radius; from = 'dilution volume (' + d.Vz.toPrecision(3) + ' m³)'; } else problems.push(problem('area.extent', 'warning', `Source ${s.id}: no extent. Give the release rate (or a hole), the lower flammable limit and the air changes per second, or enter the distance from the standard's chart.`, about)); }
    if (!r) continue;
    if (!s.ventilation.indoor && !s.extent?.manual) problems.push(problem('area.outdoor', 'note', `Source ${s.id} is outdoors: the dilution volume method is for ventilated spaces; take the distance from the standard's outdoor chart and enter it as the extent.`, about));
    if (s.substance.heavy) problems.push(problem('area.heavy', 'warning', `Source ${s.id} releases a vapour heavier than air: the zone extends low, into pits and trenches. The shape below is a ground-to-release-height bounding cylinder.`, about));
    const shape = s.substance.heavy ? { type: 'cylinder', cx: s.at.x, cy: s.at.y, r, z0: 0, z1: s.at.z + r } : { type: 'sphere', cx: s.at.x, cy: s.at.y, cz: s.at.z, r };
    zones.push({ id: s.id, source: s.id, tag: s.tag, zone: z, grade: s.grade, extent: r, extentFrom: from, shape, gasGroup: s.gasGroup, tempClass: s.tempClass });
  }
  return { zones, problems, notes };
}

/** Section through the vertical plane x = const or y = const (axis 'x' | 'y', at): the chord of each shape in (s, z) with s along the other axis. */
export function sectionShapes(zones, { axis = 'x', at = 0 } = {}) {
  const out = [];
  for (const z of zones) if (z.shape) {
    const c = axis === 'x' ? z.shape.cx : z.shape.cy, o = axis === 'x' ? z.shape.cy : z.shape.cx, d = Math.abs(at - c); if (d >= z.shape.r) continue;
    const half = Math.sqrt(z.shape.r ** 2 - d ** 2);
    if (z.shape.type === 'cylinder') out.push({ id: z.id, layer: LAYER[z.zone], label: z.zone, type: 'rect', s0: o - half, s1: o + half, z0: z.shape.z0, z1: z.shape.z1 });
    else out.push({ id: z.id, layer: LAYER[z.zone], label: z.zone, type: 'circle-chord', s: o, z: z.shape.cz, r: half });
  }
  return out;
}

const inside = (sh, p) => (sh.type === 'sphere' ? Math.hypot(p.x - sh.cx, p.y - sh.cy, p.z - sh.cz) <= sh.r : Math.hypot(p.x - sh.cx, p.y - sh.cy) <= sh.r && p.z >= sh.z0 && p.z <= sh.z1);

const ring = (cx, cy, r, n = 48) => Array.from({ length: n + 1 }, (_, k) => { const a = 2 * Math.PI * (k % n) / n; return [f6(cx + r * Math.cos(a)), f6(cy + r * Math.sin(a))]; });

export function siteGeoJson(zones, { origin = [0, 0], name = 'consequence zones' } = {}) {
  const features = [];
  for (const z of zones || []) {
    const r = footprint(z); if (!r) continue;
    features.push({ type: 'Feature', properties: { id: String(z.id), label: String(z.label || ''), kind: z.type === 'circle' ? 'radius' : 'plume', level: z.level ?? null, radius_m: f6(z.r || 0), review: 'not engineer-reviewed' },
      geometry: { type: 'Polygon', coordinates: [r.map(([x, y]) => [f6(x + origin[0]), f6(y + origin[1])])] } });
  }
  return { type: 'FeatureCollection', name, features };
}

export function siteDxf(zones, { origin = [0, 0] } = {}) {
  const lv = [...new Set((zones || []).map((z) => String(z.level)))], out = ['0', 'SECTION', '2', 'ENTITIES'];
  for (const z of zones || []) {
    const r = footprint(z); if (!r) continue;
    const layer = 'CONSEQ-' + (lv.indexOf(String(z.level)) + 1);
    out.push('0', 'POLYLINE', '8', layer, '66', '1', '70', '1');
    for (const [x, y] of r.slice(0, -1)) out.push('0', 'VERTEX', '8', layer, '10', String(f6(x + origin[0])), '20', String(f6(y + origin[1])), '30', '0');
    out.push('0', 'SEQEND', '8', layer);
  }
  out.push('0', 'ENDSEC', '0', 'EOF'); return out.join('\n') + '\n';
}

/** A closed outward-facing mesh of a zone volume. volume: { type 'sphere'|'hemisphere'|'cylinder'|'polygon', cx, cy, cz | z, r, z0, z1, polygon[[x,y]], height }. Metres. */
export function volumeMesh(v, { sides = 24 } = {}) {
  const verts = [], tris = [], push = (p) => verts.push(p) - 1;
  const lathe = (profile) => {                                                 // [r, z] from the bottom pole up; about (cx, cy)
    const rings = profile.map(([r, z]) => (r < 1e-9 ? { pole: push([v.cx, v.cy, z]) } : { ring: Array.from({ length: sides }, (_, i) => push([v.cx + r * Math.cos(2 * Math.PI * i / sides), v.cy + r * Math.sin(2 * Math.PI * i / sides), z])) }));
    for (let k = 0; k + 1 < rings.length; k++) for (let i = 0; i < sides; i++) {
      const a = rings[k], b = rings[k + 1], j = (i + 1) % sides;
      if (a.pole !== undefined && b.ring) tris.push([a.pole, b.ring[j], b.ring[i]]);
      else if (a.ring && b.pole !== undefined) tris.push([a.ring[i], a.ring[j], b.pole]);
      else if (a.ring && b.ring) { tris.push([a.ring[i], a.ring[j], b.ring[j]]); tris.push([a.ring[i], b.ring[j], b.ring[i]]); }
    }
  };
  if (v.type === 'sphere' || v.shape === 'sphere') { const n = 12; lathe(Array.from({ length: n + 1 }, (_, i) => { const a = -Math.PI / 2 + Math.PI * i / n; return [i === 0 || i === n ? 0 : v.r * Math.cos(a), v.cz + v.r * Math.sin(a)]; })); }
  else if (v.type === 'hemisphere' || v.shape === 'hemisphere') { const n = 8, z0 = v.z ?? 0; lathe([[0, z0]].concat(Array.from({ length: n + 1 }, (_, i) => { const a = Math.PI / 2 * i / n; return [i === n ? 0 : v.r * Math.cos(a), z0 + v.r * Math.sin(a)]; }).slice(0))); }
  else if (v.type === 'cylinder' || v.shape === 'cylinder') lathe([[0, v.z0], [v.r, v.z0], [v.r, v.z1], [0, v.z1]]);
  else if ((v.type === 'polygon' || v.shape === 'extrusion') && v.polygon && v.polygon.length > 2) {
    const z0 = v.z0 ?? 0, z1 = z0 + (v.height ?? v.zMax ?? 1), P = v.polygon, area = P.reduce((s, p, i) => { const q = P[(i + 1) % P.length]; return s + (p[0] * q[1] - q[0] * p[1]); }, 0);
    const poly = area < 0 ? P.slice().reverse() : P, lo = poly.map(([x, y]) => push([x, y, z0])), hi = poly.map(([x, y]) => push([x, y, z1])), n = poly.length;
    for (let i = 0; i < n; i++) { const j = (i + 1) % n; tris.push([lo[i], lo[j], hi[j]], [lo[i], hi[j], hi[i]]); }
    const cLo = push([poly.reduce((s, p) => s + p[0], 0) / n, poly.reduce((s, p) => s + p[1], 0) / n, z0]), cHi = push([verts[cLo][0], verts[cLo][1], z1]);   // fan caps (star-shaped footprints)
    for (let i = 0; i < n; i++) { const j = (i + 1) % n; tris.push([cLo, lo[j], lo[i]], [cHi, hi[i], hi[j]]); }
  } else throw new Error('volume: sphere | hemisphere | cylinder | polygon');
  return { verts, tris };
}

export function meshVolume(m) { let s = 0; for (const [a, b, c] of m.tris) { const p = m.verts[a], q = m.verts[b], r = m.verts[c]; s += (p[0] * (q[1] * r[2] - q[2] * r[1]) - p[1] * (q[0] * r[2] - q[2] * r[0]) + p[2] * (q[0] * r[1] - q[1] * r[0])) / 6; } return s; }

/** opts { substances: { <service>: { name, lfl, M, heavy? } }, ventilation: { degree, availability, C?, f?, indoor? }, hole (m, default 0.00025), unitMm (drawing unit in mm → metres via x/1000), z (release height m, default 1), only: [syms] } */
export function sourcesFromPid(items, opts = {}) {
  const sources = [], problems = [], assumptions = ['grade of release by what the item is (' + 'seals, flanges and valve stems secondary; vents and sampling points primary): a screening default, not the standard\'s table', 'hole diameter ' + (opts.hole ?? 0.00025) * 1000 + ' mm unless the item gives a rate'], skipped = new Set();
  const vent = { degree: 'medium', availability: 'fair', indoor: false, ...(opts.ventilation || {}) }, k = opts.toMetres ?? 1;
  let n = 0;
  for (const it of items || []) {
    const d = DEFAULT_GRADE[it.sym]; if (!d || (opts.only && !opts.only.includes(it.sym))) continue;
    const svc = (it.svc || it.attrs?.svc || '').trim(), sub = svc && opts.substances?.[svc];
    if (!sub) { skipped.add(svc || '(no service)'); continue; }
    const vented = it.attrs?.vent === 'atmosphere' || it.attrs?.vent === 'yes', P = num(it.attrs?.dp), grade = it.sym === 'relief' && vented ? 'primary' : d[1];
    const src = { id: 'PS' + (++n), tag: it.tag || '', kind: d[0], grade, substance: { name: sub.name, lfl: sub.lfl, M: sub.M, heavy: !!sub.heavy }, at: { x: f6((it.x ?? 0) * k), y: f6((it.y ?? 0) * k), z: opts.z ?? 1 }, ventilation: { ...vent }, pid: it.tag || '' };
    if (it.attrs?.rate && num(it.attrs.rate) > 0) src.rate = num(it.attrs.rate);
    else if (P && (opts.hole ?? 0.00025) > 0) src.hole = { d: opts.hole ?? 0.00025, P0: 101325 + P * 1e5, T: (num(it.attrs?.dt) != null ? num(it.attrs.dt) + 273.15 : 293.15), k: sub.k || 1.3 };
    else { problems.push({ code: 'area.nosize', severity: 'info', text: `${it.tag || it.sym}: no design pressure (dp) or release rate (rate) on the symbol, so no extent can be worked out; the source is listed without one.` }); }
    sources.push(src);
  }
  for (const s of skipped) problems.push({ code: 'area.nosubstance', severity: 'warning', text: `Service "${s}" has no substance in the list given (name, lower flammable limit, molar mass): its items are not sources yet.` });
  return { sources, problems, assumptions };
}