/* mna-engine.jsx — a small TRANSIENT circuit simulator for the free-build
   oscillator bench. Unlike breadboard.jsx's static topology recognizer, this
   actually steps through time: capacitors integrate, the 555 flip-flop holds
   state, and an astable genuinely oscillates.

   Formulation: nodal analysis with the Norton/companion method — every element
   is reduced to (conductance G between two nodes) + (current source I into a
   node). Voltage sources are stamped as stiff Thevenin (small series R), so the
   only unknowns are node voltages and each step is a single dense linear solve
   G·v = i, wrapped in a few fixed-point passes for the nonlinear LED/diode/BJT.

   No React, no DOM — pure functions, exported on window for the bench + tests.
   Names mna* prefixed. */

const MNA_GND = 0;          // net id 0 is always ground (reference, V=0)
const MNA_VT = 0.05;        // diode thermal-ish scale (tuned, not physical)

/* ── dense linear solve (Gaussian elimination w/ partial pivot) ──────────── */
function mnaSolve(A, b) {
  const n = b.length;
  // augment
  const M = A.map((row, i) => row.concat(b[i]));
  for (let col = 0; col < n; col++) {
    // pivot
    let piv = col;
    for (let r = col + 1; r < n; r++) if (Math.abs(M[r][col]) > Math.abs(M[piv][col])) piv = r;
    if (Math.abs(M[piv][col]) < 1e-12) { M[piv][col] += 1e-9; }  // nudge singular
    if (piv !== col) { const t = M[piv]; M[piv] = M[col]; M[col] = t; }
    const d = M[col][col];
    for (let r = 0; r < n; r++) {
      if (r === col) continue;
      const f = M[r][col] / d;
      if (f === 0) continue;
      for (let c = col; c <= n; c++) M[r][c] -= f * M[col][c];
    }
  }
  const x = new Array(n);
  for (let i = 0; i < n; i++) x[i] = M[i][n] / M[i][i];
  return x;
}

/* ── the simulator ───────────────────────────────────────────────────────
   netlist: {
     nNodes,                       // total nets incl. ground (ids 0..nNodes-1)
     elements: [ ... ]             // see element shapes below
   }
   Element shapes (nodes are net ids; ground = 0):
     { k:"V",   a, b, volts, r }                  ideal-ish source a→b (Thevenin r)
     { k:"R",   a, b, ohms }
     { k:"C",   a, b, farads }                    // state Vc carried internally
     { k:"D",   a, b, vf, ron }                   // diode a(anode)→b(cathode)
     { k:"LED", a, b, vf, ron }                   // like D, + brightness readout
     { k:"NPN", b, c, e }                          // base/collector/emitter
     { k:"IC555", vcc, gnd, out, dis, thr, trig, rst, ctrl }
*/
function mnaCreate(netlist) {
  const els = netlist.elements;
  const n = netlist.nNodes;

  // persistent state
  const capV = els.map(e => (e.k === "C" ? 0 : null));         // last cap voltage
  const capI = els.map(e => (e.k === "C" ? 0 : null));         // last cap (displacement) current
  const indI = els.map(e => (e.k === "L" ? 0 : null));         // last inductor current
  const indVp = els.map(e => (e.k === "L" ? 0 : null));        // last inductor voltage (trapezoidal companion)
  const ff   = { q: false };                                    // 555 flip-flop (per IC; supports one — fine here)
  const ffs  = els.map(e => (e.k === "IC555" ? { q: false } : null));
  let icPrevV = null;   // last settled solution — clocks the 555 latch (see step)
  // nonlinear element operating voltages (for warm-start)
  const nlV  = els.map(() => 0);

  function step(dt) {
    // ── fixed-point passes for nonlinear devices ──
    let v = new Array(n).fill(0);
    let lastSettled = false;

    // Clock every 555 latch ONCE per timestep, from the PREVIOUS settled
    // solution. (Updating the comparators inside the fixed-point passes below —
    // on half-solved node voltages — made the latch flip every pass, which froze
    // the timing cap at ⅓·Vcc and chattered OUT at the step rate. dt ≪ R·C, so a
    // one-step evaluation delay is invisible.)
    if (icPrevV) {
      els.forEach((e, idx) => {
        if (e.k !== "IC555") return;
        const vcc = (e.vcc === MNA_GND ? 0 : icPrevV[e.vcc]) || 0;
        const supply = Math.max(0.1, vcc);
        const hi = 2 / 3 * supply, lo = 1 / 3 * supply;
        const vthr = (e.thr === MNA_GND ? 0 : icPrevV[e.thr]) || 0;
        const vtrig = (e.trig === MNA_GND ? 0 : icPrevV[e.trig]) || 0;
        const vrst = (e.rst == null || e.rst === MNA_GND) ? 0 : (icPrevV[e.rst] || 0);
        const st = ffs[idx];
        if (e.rst != null && vrst < 0.4 * supply) st.q = false;   // reset dominates
        else {
          if (vthr > hi) st.q = false;        // threshold → reset (out low)
          if (vtrig < lo) st.q = true;         // trigger → set  (out high)
        }
      });
    }
    for (let pass = 0; pass < 24; pass++) {
      const G = Array.from({ length: n }, () => new Array(n).fill(0));
      const I = new Array(n).fill(0);
      const stampG = (a, b, g) => {
        if (a !== MNA_GND) G[a][a] += g;
        if (b !== MNA_GND) G[b][b] += g;
        if (a !== MNA_GND && b !== MNA_GND) { G[a][b] -= g; G[b][a] -= g; }
      };
      const stampI = (a, b, cur) => {          // current source from a→b (push into b, out of a)
        if (a !== MNA_GND) I[a] -= cur;
        if (b !== MNA_GND) I[b] += cur;
      };

      els.forEach((e, idx) => {
        if (e.k === "R") {
          stampG(e.a, e.b, 1 / Math.max(1e-3, e.ohms));
        } else if (e.k === "V") {
          const g = 1 / Math.max(1e-3, e.r || 0.05);
          stampG(e.a, e.b, g);
          stampI(e.b, e.a, e.volts * g);        // Thevenin → Norton
        } else if (e.k === "C") {
          // trapezoidal companion (energy-conserving — backward Euler would add
          // numerical damping that kills LC ringing). Geq=2C/dt, Ieq=Geq·Vprev+Iprev.
          const geq = 2 * e.farads / dt;
          stampG(e.a, e.b, geq);
          stampI(e.b, e.a, geq * capV[idx] + capI[idx]);
        } else if (e.k === "L") {
          // trapezoidal companion (dual of the capacitor): geq=dt/(2L) in parallel
          // with a source carrying the prior current plus the prior-voltage term.
          const geq = dt / (2 * Math.max(1e-9, e.henries));
          stampG(e.a, e.b, geq);
          stampI(e.a, e.b, indI[idx] + geq * indVp[idx]);
        } else if (e.k === "D" || e.k === "LED") {
          // PWL diode: conduct (Ron, offset Vf) when forward, else tiny leak.
          const va = (e.a === MNA_GND ? 0 : v[e.a]);
          const vb = (e.b === MNA_GND ? 0 : v[e.b]);
          const vd = va - vb;
          // turn-on at exactly Vf (not a soft Vf−0.15 knee): the soft knee let a
          // reverse-biased diode leak a small NEGATIVE current near the threshold,
          // which showed up as phantom reverse flow (e.g. the inverter's "off" LED).
          // A real diode blocks reverse, so conduct only when forward past Vf.
          const on = vd > e.vf;
          if (on) {
            const g = 1 / Math.max(1, e.ron);
            stampG(e.a, e.b, g);
            // companion offset source: inject +g·Vf INTO the anode (e.a), out of
            // the cathode (e.b), so branch current = g·(Va−Vb) − g·Vf = (Vd−Vf)/Ron.
            // (stampI(x,y,c) does I[x]−=c, I[y]+=c, so x=cathode, y=anode.)
            stampI(e.b, e.a, e.vf * g);
          } else {
            stampG(e.a, e.b, 1e-9);
          }
        } else if (e.k === "NPN") {
          const vb = (e.b === MNA_GND ? 0 : v[e.b]);
          const vc = (e.c === MNA_GND ? 0 : v[e.c]);
          const ve = (e.e === MNA_GND ? 0 : v[e.e]);
          const vbe = vb - ve;
          const beta = e.beta || 120;
          if (vbe > 0.5) {
            // B-E as a (PWL) diode — its on-resistance + the external base resistor
            // set the base current Ib. Offset source injects INTO the base (anode).
            const gbe = 1 / 200;
            stampG(e.b, e.e, gbe);
            stampI(e.e, e.b, 0.65 * gbe);
            const ib = Math.max(0, vbe - 0.65) * gbe;
            const icA = beta * ib;                 // active-region collector current = β·Ib
            const vce = vc - ve;
            if (vce > icA * 40) {
              // active: a β-controlled current source C→E, so the base current
              // genuinely controls the collector current (change Rb → LED dims).
              stampI(e.c, e.e, icA);
              stampG(e.c, e.e, 1e-5);              // tiny output conductance (numerical stability)
            } else {
              stampG(e.c, e.e, 1 / 40);            // saturation: Vce collapses, Rce small
            }
          } else {
            stampG(e.b, e.e, 1e-9);
            stampG(e.c, e.e, 1e-9);
          }
        } else if (e.k === "IC555") {
          const vcc = (e.vcc === MNA_GND ? 0 : v[e.vcc]);
          const supply = Math.max(0.1, vcc);
          const st = ffs[idx];   // latch already clocked for this timestep (above)
          // OUT pin: push-pull Thevenin to supply or 0
          const gout = 1 / 12;
          if (e.out != null) {
            stampG(e.out, MNA_GND, gout);
            stampI(MNA_GND, e.out, (st.q ? supply : 0) * gout);
          }
          // DISCHARGE pin: closed to ground when out low, else open
          if (e.dis != null) {
            if (!st.q) stampG(e.dis, MNA_GND, 1 / 12);
            else stampG(e.dis, MNA_GND, 1e-9);
          }
        } else if (e.k === "OPAMP") {
          // Op-amp as a voltage-controlled voltage source: Vout = A·(v+ − v−),
          // clamped to the rails. Stamped LINEARLY into the matrix as a
          // transconductance (Gout behind gain A) so any feedback network solves
          // in one matrix solve — no relaxation needed. The rail clamp is the only
          // nonlinearity and is handled by the fixed-point loop: once the output
          // saturates we switch to a stiff source pinned at the rail.
          const vp = (e.np === MNA_GND ? 0 : v[e.np]);
          const vn = (e.nn === MNA_GND ? 0 : v[e.nn]);
          const vhi = (e.vpos == null) ? 12 : (typeof e.vpos === "number" ? e.vpos : (e.vpos === MNA_GND ? 0 : v[e.vpos]));
          const vlo = (e.vneg == null) ? 0  : (typeof e.vneg === "number" ? e.vneg : (e.vneg === MNA_GND ? 0 : v[e.vneg]));
          const A = e.gain || 1e5;
          const gout = 1 / (e.rout || 50);
          const open = (vp - vn) * A;                 // ideal open-loop output
          if (e.out != null && e.out !== MNA_GND) {
            G[e.out][e.out] += gout;                  // output conductance to ref
            if (open > vhi) {
              I[e.out] += vhi * gout;                 // saturated HIGH: stiff source at +rail
            } else if (open < vlo) {
              I[e.out] += vlo * gout;                 // saturated LOW: stiff source at −rail
            } else {
              // linear region: Vout = A·(v+ − v−) ⇒ controlled-source entries
              if (e.np !== MNA_GND) G[e.out][e.np] -= gout * A;
              if (e.nn !== MNA_GND) G[e.out][e.nn] += gout * A;
            }
          }
          // inputs are high-Z — weak tie so a floating input stays defined
          if (e.np != null && e.np !== MNA_GND) G[e.np][e.np] += 1e-12;
          if (e.nn != null && e.nn !== MNA_GND) G[e.nn][e.nn] += 1e-12;
        }
      });

      // ground row/col → identity (v[0] = 0)
      // build reduced system over nodes 1..n-1
      const m = n - 1;
      if (m <= 0) { v = [0]; break; }
      const A = Array.from({ length: m }, (_, i) => {
        const row = new Array(m).fill(0);
        for (let j = 0; j < m; j++) row[j] = G[i + 1][j + 1];
        return row;
      });
      const rhs = new Array(m);
      for (let i = 0; i < m; i++) rhs[i] = I[i + 1];
      const sol = mnaSolve(A, rhs);
      const vNew = [0, ...sol];

      // convergence check
      let maxd = 0;
      for (let i = 0; i < n; i++) maxd = Math.max(maxd, Math.abs((vNew[i] || 0) - (v[i] || 0)));
      v = vNew;
      if (maxd < 1e-4) { if (lastSettled) break; lastSettled = true; } else lastSettled = false;
    }

    // ── commit capacitor & inductor states (trapezoidal) ──
    els.forEach((e, idx) => {
      if (e.k === "C") {
        const va = (e.a === MNA_GND ? 0 : v[e.a]);
        const vb = (e.b === MNA_GND ? 0 : v[e.b]);
        const vnew = va - vb;
        const geq = 2 * e.farads / dt;
        capI[idx] = geq * (vnew - capV[idx]) - capI[idx];   // i = Geq·ΔV − Iprev
        capV[idx] = vnew;
      } else if (e.k === "L") {
        const va = (e.a === MNA_GND ? 0 : v[e.a]);
        const vb = (e.b === MNA_GND ? 0 : v[e.b]);
        const vnew = va - vb;
        const geq = dt / (2 * Math.max(1e-9, e.henries));
        indI[idx] = indI[idx] + geq * (vnew + indVp[idx]);  // I += dt/(2L)·(Vnew+Vprev)
        indVp[idx] = vnew;
      }
    });

    icPrevV = v;   // remember settled solution to clock the 555 latch next step
    return v;
  }

  // element-current readout (post-step), for the scope/LED brightness
  function elementCurrent(idx, v) {
    const e = els[idx];
    const va = (e.a === MNA_GND ? 0 : v[e.a]);
    const vb = (e.b === MNA_GND ? 0 : v[e.b]);
    if (e.k === "R") return (va - vb) / Math.max(1e-3, e.ohms);
    if (e.k === "V") return (e.volts - (va - vb)) / Math.max(1e-3, e.r || 0.05);  // Thevenin branch current
    if (e.k === "L") return indI[idx];
    if (e.k === "D" || e.k === "LED") {
      const vd = va - vb;
      return vd > e.vf ? (vd - e.vf) / Math.max(1, e.ron) : 0;
    }
    if (e.k === "C") return capI[idx] || 0;
    if (e.k === "NPN") {
      const vbb = (e.b === MNA_GND ? 0 : v[e.b]);
      const vcc = (e.c === MNA_GND ? 0 : v[e.c]);
      const vee = (e.e === MNA_GND ? 0 : v[e.e]);
      const vbe = vbb - vee, vce = vcc - vee;
      if (vbe > 0.5 && vce > 0.02) {
        const beta = e.beta || 120;
        const ib = Math.max(0, vbe - 0.65) / 200;
        return Math.min(beta * ib, vce / 40);    // active (β·Ib) or saturation-limited
      }
      return 0;
    }
    return 0;
  }

  function ledState(idx) { return ffs[idx]; }
  function getFF(idx) { return ffs[idx] ? ffs[idx].q : null; }
  function setCapV(idx, val) { capV[idx] = val; }
  function getCapV(idx) { return capV[idx]; }
  function setIndI(idx, val) { indI[idx] = val; }
  function getIndI(idx) { return indI[idx]; }
  function reset() { capV.forEach((_, i) => { if (capV[i] != null) capV[i] = 0; }); capI.forEach((_, i) => { if (capI[i] != null) capI[i] = 0; }); indI.forEach((_, i) => { if (indI[i] != null) indI[i] = 0; }); indVp.forEach((_, i) => { if (indVp[i] != null) indVp[i] = 0; }); ffs.forEach(s => { if (s) s.q = false; }); icPrevV = null; }

  return { step, elementCurrent, getFF, setCapV, getCapV, setIndI, getIndI, reset, els, nNodes: n };
}

Object.assign(window, { mnaCreate, mnaSolve, MNA_GND });