Layer of protection and SIL arithmetic

/** Mitigated frequency per year: sum over initiating events of f x prod(modifier p) x prod(IPL PFD). Missing inputs are listed, never assumed. */
export function scenarioResult(sc) {
  const missing = [];
  if (!sc.initiating.length) missing.push('an initiating event');
  for (const e of sc.initiating) if (!(e.f > 0)) missing.push('frequency of ' + (e.name || e.id));
  for (const m of sc.modifiers) if (m.p === undefined) missing.push('probability of ' + (m.name || 'a modifier'));
  for (const p of sc.ipls) if (!(p.pfd > 0)) missing.push('PFD of ' + (p.name || p.id));
  const fIE = sc.initiating.reduce((s, e) => s + (e.f > 0 ? e.f : 0), 0), pm = sc.modifiers.reduce((a, m) => a * (m.p ?? 1), 1), pi = sc.ipls.reduce((a, p) => a * (p.pfd > 0 ? p.pfd : 1), 1);
  const f = fIE * pm * pi, t = sc.target;
  const out = { id: sc.id, complete: !missing.length, missing, fIE, modifiers: pm, ipls: pi, fMitigated: f, target: t ?? null };
  if (t > 0) { out.meets = f <= t; out.gap = f / t; out.rrf = f > t ? f / t : 1; out.requiredPfd = f > t ? t / f : 1; }
  return out;
}

/** The PFD a new layer (a SIF) must reach: target / (mitigated frequency without it). */
export function requiredPfd(sc, withoutSifId) {
  const copy = clone(sc); copy.ipls = copy.ipls.filter((p) => p.id !== withoutSifId); const r = scenarioResult(copy);
  return { requiredPfd: r.requiredPfd ?? null, silTarget: r.requiredPfd ? silBand(r.requiredPfd) : null, fWithout: r.fMitigated };
}

/** SIL for a PFDavg: the band it falls in; 0 when above SIL1's upper edge (no SIL), null when below SIL4's lower edge ("better than SIL4": not claimed). */
export function silBand(pfd) {
  if (!(pfd > 0)) return null;
  if (pfd >= 1e-1) return 0;
  for (const b of SIL_LOW) if (pfd >= b.lo && pfd < b.hi) return b.sil;
  return pfd < 1e-5 ? 4 : 0;
}

/** Unavailability of ONE channel at time t since new: dangerous undetected failures that the proof test finds (rate PTC lDU) are cleared every TI; the rest (1-PTC) stay until the end of life LT. */
export function pfdChannel(t, { lambdaDU, TI, PTC = 1, beta = 0 }) {
  const lt = (1 - beta) * lambdaDU * PTC, ln = (1 - beta) * lambdaDU * (1 - PTC), tau = t % TI;
  return 1 - Math.exp(-lt * tau - ln * t);
}

/** PFDavg by the closed forms (needs PTC = 1). */
export function pfdAvgFormula({ arch, lambdaDU, TI, beta = 0 }) {
  const l = lambdaDU * TI, a = (1 - beta) * l, c = beta * l / 2;
  switch (arch) { case '1oo1': return l / 2; case '1oo2': return a * a / 3 + c; case '2oo2': return l; case '1oo3': return a ** 3 / 4 + c; case '2oo3': return a * a + c; default: throw new Error('Unknown architecture ' + arch); }
}

/** PFDavg by averaging the exact channel unavailability over the life (Simpson), with KooN voting and the beta part as a common-cause element added. Handles PTC < 1. */
export function pfdAvgNumeric({ arch, lambdaDU, TI, beta = 0, PTC = 1, LT }) {
  if (!ARCH[arch]) throw new Error('Unknown architecture ' + arch);
  const [K, N] = ARCH[arch]; if (K === N) beta = 0;   // beta has no meaning with one channel, and in a 2oo2 any single failure already defeats the function
  const life = LT && LT >= TI ? LT : TI, steps = Math.max(200, Math.ceil(life / TI) * 400) & ~1, h = life / steps;
  const f = (t) => {
    const u = pfdChannel(Math.min(t, life - 1e-9), { lambdaDU, TI, PTC, beta });
    let p = 0; for (let j = N - K + 1; j <= N; j++) p += choose(N, j) * u ** j * (1 - u) ** (N - j);   // at least N-K+1 channels failed dangerously
    const ccf = N > 1 ? beta * lambdaDU * (PTC * (t % TI) + (1 - PTC) * t) : 0;   // first order common-cause term
    return p + ccf;
  };
  let s = f(0) + f(life); for (let i = 1; i < steps; i++) s += f(i * h) * (i % 2 ? 4 : 2);
  const v = s * h / 3 / life;
  return v;
}

/** A SIF: PFDavg per subsystem and in total (sum), the SIL of the total and of the weakest subsystem; notes list what was not checked. */
export function verifySif(sif, { target = null } = {}) {
  const subs = sif.subsystems.map((u) => {
    const ok = u.lambdaDU >= 0 && u.TI > 0;
    if (!ok) return { name: u.name, arch: u.arch, pfd: null, missing: [u.lambdaDU >= 0 ? '' : 'dangerous undetected failure rate', u.TI > 0 ? '' : 'proof-test interval'].filter(Boolean) };
    const pfd = pfdAvgNumeric({ arch: u.arch, lambdaDU: u.lambdaDU, TI: u.TI, beta: u.beta || 0, PTC: u.PTC ?? 1, LT: u.LT });
    return { name: u.name, arch: u.arch, pfd, formula: (u.PTC ?? 1) === 1 ? pfdAvgFormula({ arch: u.arch, lambdaDU: u.lambdaDU, TI: u.TI, beta: u.beta || 0 }) : null, sil: silBand(pfd), HFT: ARCH[u.arch][1] - ARCH[u.arch][0] };
  });
  const complete = subs.length > 0 && subs.every((s) => s.pfd != null), total = complete ? subs.reduce((a, s) => a + s.pfd, 0) : null;
  const out = { id: sif.id, complete, subsystems: subs, pfdAvg: total, sil: total == null ? null : silBand(total), demand: sif.demand,
    notes: ['Architectural constraints (hardware fault tolerance, safe failure fraction, route 1H or 2H) and systematic capability are not checked.', 'Failure rates and test intervals are the person\'s inputs; none is assumed.', REVIEW.note] };
  if (sif.demand === 'high') out.notes.push('High-demand mode needs PFH, not PFDavg: not computed here.');
  if (target != null && total != null) { out.meetsTarget = total <= target; out.targetPfd = target; }
  return out;
}