"use strict";

// Transmission line physics — pure functions, no DOM or canvas dependencies.
//
// Spatial convention: z ∈ [0, 1] maps to [0, ℓ].
// Time convention:    tNorm is dimensionless time in units of τ_d (one-way delay).
//
// Single-segment model (N=1):
//   model.segments = [{Z0: 50}] — one segment, the whole line.
//   Degenerates exactly to the original two-terminal bounce series.
//
// Multi-segment model (N>1):
//   model.segments = [{Z0: Z1}, {Z0: Z2}, …]  — N equal-length segments.
//   Each internal boundary produces a reflected and a transmitted wave (lattice diagram).
//   buildBounceSeries uses a BFS priority queue ordered by tDie.
//
// Wave packet struct (new):
//   { n, A, dir, zStart, zEnd, tBorn, tDie, segIdx }
//   n       — sequential index (for display)
//   A       — voltage amplitude
//   dir     — +1 rightward, -1 leftward
//   zStart  — z where this packet's front begins its journey
//   zEnd    — z where this packet's front ends (next boundary, or 0/1)
//   tBorn   — tNorm when the front was at zStart
//   tDie    — tNorm when the front reaches zEnd  (= tBorn + |zEnd − zStart|)
//   segIdx  — index into model.segments for the segment this packet lives in
//
// After computeDynamicState, each wave also carries:
//   u       — tNorm − tBorn  (age of the packet)
//   front   — current position of the leading edge  ∈ [zStart, zEnd] (or reversed for leftward)
//
// Extension points:
//   • Rise time (causal shifted erf or linear ramp): model.riseTimeTr > 0  — see riseShape/waveVoltageAt.
//   • Non-equal segment lengths: store lengths[] summing to 1; the BFS solver already works
//     by position, so changing zBnds is all that's needed.

const TLPhysics = (() => {
  const BOUNCE_EPS  = 1e-6;
  const MAX_BOUNCES = 5000;

  // ---- reflection coefficients and initial voltage step ----
  // Uses the first segment's Z0 for the source divider and source-end Γ,
  // and the last segment's Z0 for the load-end Γ.
  function computeWaveParams(model) {
    const { Vg, Rg, RL } = model;
    const Z0_src  = model.segments[0].Z0;
    const Z0_load = model.segments[model.segments.length - 1].Z0;
    const V1 = Vg * Z0_src / (Rg + Z0_src);
    const gL = isFinite(RL) ? (RL - Z0_load) / (RL + Z0_load) : 1;   // RL=∞ → open circuit
    const gS = (Rg - Z0_src)  / (Rg + Z0_src);
    return { V1, gL, gS };
  }

  // ---- BFS lattice-diagram solver ----
  // Builds the full set of wave packets produced by successive reflections and
  // transmissions at every impedance boundary (source, internal, load).
  //
  // Returns {
  //   series:     [{n, A, dir, zStart, zEnd, tBorn, tDie, segIdx}]
  //   srcEvents:  [{t, dV}]   — cumulative voltage steps at z=0  (VS)
  //   loadEvents: [{t, dV}]   — cumulative voltage steps at z=ℓ  (VL)
  //   gSEffective, tolPct, ampTol, tEnd
  // }
  function buildBounceSeries(model, waves) {
    const segs   = model.segments;
    const N      = segs.length;
    const segLen = 1 / N;
    // z positions of all N+1 boundaries: [0, 1/N, 2/N, ..., 1]
    const zBnds  = Array.from({ length: N + 1 }, (_, i) => i * segLen);

    const gS     = waves.gS;
    const gL     = waves.gL;
    const tolPct = Number.isFinite(model.reflectTol) ? model.reflectTol : 1;
    const ampTol = Math.max(BOUNCE_EPS, Math.abs(waves.V1) * (tolPct / 100));

    const series     = [];
    const srcEvents  = [{ t: 0, dV: waves.V1 }]; // initial step seen at source at t=0
    const loadEvents = [];

    // Queue of pending packets; each entry is a wave packet description.
    // We process in order of tDie (earliest-completing first).
    const queue = [];

    function enqueue(pkt) {
      if (Math.abs(pkt.A) <= ampTol) return;
      if (series.length + queue.length >= MAX_BOUNCES) return;
      queue.push(pkt);
    }

    // Seed: first incident wave in segment 0, rightward from z=0 to z=zBnds[1].
    enqueue({ A: waves.V1, dir: +1, zStart: 0, zEnd: zBnds[1], tBorn: 0, segIdx: 0 });

    while (queue.length > 0) {
      // Pop earliest-completing packet (smallest tDie = tBorn + |zEnd − zStart|).
      queue.sort((a, b) =>
        (a.tBorn + Math.abs(a.zEnd - a.zStart)) -
        (b.tBorn + Math.abs(b.zEnd - b.zStart))
      );
      const pkt = queue.shift();
      const { A, dir, zStart, zEnd, tBorn, segIdx } = pkt;
      const tDie = tBorn + Math.abs(zEnd - zStart);
      const n    = series.length + 1;
      series.push({ n, A, dir, zStart, zEnd, tBorn, tDie, segIdx });

      const atSource = zEnd < 1e-9;
      const atLoad   = zEnd > 1 - 1e-9;

      if (atLoad) {
        // Load terminal: add load event, spawn reflected wave back in last segment.
        loadEvents.push({ t: tDie, dV: (1 + gL) * A });
        enqueue({
          A: gL * A, dir: -1,
          zStart: 1, zEnd: zBnds[N - 1],
          tBorn: tDie, segIdx: N - 1,
        });
      } else if (atSource) {
        // Source terminal: add source event, spawn reflected wave into first segment.
        srcEvents.push({ t: tDie, dV: (1 + gS) * A });
        enqueue({
          A: gS * A, dir: +1,
          zStart: 0, zEnd: zBnds[1],
          tBorn: tDie, segIdx: 0,
        });
      } else {
        // Internal boundary between two segments.
        // For a rightward wave in segIdx: left=segIdx, right=segIdx+1.
        // For a leftward wave in segIdx: left=segIdx-1, right=segIdx.
        const leftSeg  = dir > 0 ? segIdx     : segIdx - 1;
        const rightSeg = dir > 0 ? segIdx + 1 : segIdx;
        const Zl = segs[leftSeg].Z0;
        const Zr = segs[rightSeg].Z0;
        // Γ for a wave arriving from the left side of this boundary:
        const gBound = (Zr - Zl) / (Zr + Zl);

        if (dir > 0) {
          // Rightward hitting boundary from left:
          //   reflected  → Γ·A leftward, back into leftSeg
          //   transmitted → (1+Γ)·A rightward, into rightSeg
          enqueue({
            A: gBound * A, dir: -1,
            zStart: zEnd, zEnd: zBnds[leftSeg],
            tBorn: tDie, segIdx: leftSeg,
          });
          enqueue({
            A: (1 + gBound) * A, dir: +1,
            zStart: zEnd, zEnd: zBnds[rightSeg + 1],
            tBorn: tDie, segIdx: rightSeg,
          });
        } else {
          // Leftward hitting boundary from right:
          //   reflected  → −Γ·A rightward, back into rightSeg
          //   transmitted → (1−Γ)·A leftward, into leftSeg
          enqueue({
            A: -gBound * A, dir: +1,
            zStart: zEnd, zEnd: zBnds[rightSeg + 1],
            tBorn: tDie, segIdx: rightSeg,
          });
          enqueue({
            A: (1 - gBound) * A, dir: -1,
            zStart: zEnd, zEnd: zBnds[leftSeg],
            tBorn: tDie, segIdx: leftSeg,
          });
        }
      }
    }

    const tEnd = series.reduce((m, w) => Math.max(m, w.tDie), 2);
    return { series, srcEvents, loadEvents, gSEffective: gS, tolPct, ampTol, tEnd };
  }

  // ---- time query: node voltage (step model) ----
  function sumEventsAtTime(events, tn) {
    const eps = 1e-9;
    let v = 0;
    for (const e of events) {
      if (tn >= e.t - eps) v += e.dV;
    }
    return v;
  }

  // ---- rise-time wave shapes ----
  // Fraction of final amplitude reached after time dt since the wave front passed.

  // Abramowitz & Stegun 7.1.26 rational approximation (max error < 1.5e-7).
  function erf(x) {
    const sign = x < 0 ? -1 : 1;
    x = Math.abs(x);
    const t = 1 / (1 + 0.3275911 * x);
    const y = 1 - (((((1.061405429 * t - 1.453152027) * t) + 1.421413741) * t - 0.284496736) * t + 0.254829592) * t * Math.exp(-x * x);
    return sign * y;
  }

  // Causal erf S-curve: tr = 0 → hard step. tr > 0 → smooth S from 0 to 1 over [0, tr].
  // Uses shifted erf: 0.5·(1 + erf(k·(2·dt/tr − 1))) with k chosen so that
  // riseShape(0, tr) ≈ 0 and riseShape(tr, tr) ≈ 1.  k = 1.8 gives < 0.1% undershoot.
  // Causal: strictly 0 for dt ≤ 0, strictly 1 for dt ≥ tr.
  const ERF_K = 1.8;
  function riseShape(dt, tr) {
    if (dt <= 0) return 0;
    if (tr <= 0 || dt >= tr) return 1;
    return 0.5 * (1 + erf(ERF_K * (2 * dt / tr - 1)));
  }

  // Linear ramp (trapezoidal / SPICE PULSE): tr = 0 → hard step. tr > 0 → dt/tr clamped to [0,1].
  function riseShapeLinear(dt, tr) {
    if (dt <= 0) return 0;
    if (tr <= 0 || dt >= tr) return 1;
    return dt / tr;
  }

  // Voltage contribution of a single wave packet wf at position z.
  // wf must carry the fields added by computeDynamicState: u, front.
  // Returns 0 if z is outside this packet's z-range or ahead of the front.
  // shape: "linear" → riseShapeLinear, "erf" → riseShape (causal shifted erf).
  // tr: rise time in τ_d units (used by both shapes).
  function waveVoltageAt(wf, z, tr = 0, shape = "erf") {
    const { dir, zStart, A, u, front } = wf;
    const eps = 1e-9;
    const rise = shape === "linear" ? riseShapeLinear : riseShape;
    if (dir > 0) {
      // Rightward: nonzero in [zStart, front].
      if (z < zStart - eps || z > front + eps) return 0;
      const dt = u - (z - zStart);
      return A * rise(dt, tr);
    } else {
      // Leftward: nonzero in [front, zStart].
      if (z > zStart + eps || z < front - eps) return 0;
      const dt = u - (zStart - z);
      return A * rise(dt, tr);
    }
  }

  // Map z ∈ [0,1] to the index of the segment that owns it.
  // Boundaries are assigned to the right-hand (higher-index) segment; z=1 goes to the last.
  // The +1e-9 nudge handles floating-point imprecision when z is exactly i/N.
  function segmentForZ(z, N) {
    if (z >= 1 - 1e-9) return N - 1;
    return Math.min(N - 1, Math.floor(z * N + 1e-9));
  }

  // Sum of all wave contributions at position z.
  // For multi-segment models (N > 1), only waves whose segIdx matches the segment
  // that owns z are summed.  This prevents double-counting at segment boundaries
  // where adjacent waves share a z value (a plotting artifact, not a physics error).
  function totalVoltageAt(z, launchedWaves, tr = 0, shape = "erf", N = 1) {
    const segIdx = N > 1 ? segmentForZ(z, N) : 0;
    let V = 0;
    for (const wf of launchedWaves) {
      if (N > 1 && wf.segIdx !== segIdx) continue;
      V += waveVoltageAt(wf, z, tr, shape);
    }
    return V;
  }

  // Node voltage as smooth sum of causal erf-rise events.
  function sumEventsWithRise(events, tn, tr) {
    let v = 0;
    for (const e of events) {
      const dt = tn - e.t;
      if (dt > 0) v += e.dV * riseShape(dt, tr);
    }
    return v;
  }

  // Node voltage as smooth sum of linear-ramp events (matches SPICE PULSE with finite TR).
  function sumEventsWithLinearRamp(events, tn, tr) {
    let v = 0;
    for (const e of events) {
      const dt = tn - e.t;
      if (dt > 0) v += e.dV * riseShapeLinear(dt, tr);
    }
    return v;
  }

  // ---- dynamic state at a given tNorm ----
  // Annotates each launched wave with:
  //   u:     tNorm − tBorn
  //   front: current leading-edge position
  // A wave is "active" while u < |zEnd − zStart| (front still in transit).
  // shape: "step" | "linear" | "erf"
  function computeDynamicState(tn, bounce, riseTimeTr = 0, riseShape_ = "erf") {
    const { clamp } = TLUtils;
    const launchedWaves = [];
    const activeWaves   = [];

    for (const w of bounce.series) {
      const u = tn - w.tBorn;
      if (u < 0) continue;
      let front;
      if (w.dir > 0) {
        front = clamp(w.zStart + u, w.zStart, w.zEnd);
      } else {
        front = clamp(w.zStart - u, w.zEnd, w.zStart);
      }
      const ww = { ...w, u, front };
      launchedWaves.push(ww);
      if (u < Math.abs(w.zEnd - w.zStart) - 1e-9) activeWaves.push(ww);
    }

    const sumNode = (events) => {
      if (riseTimeTr <= 0) return sumEventsAtTime(events, tn);
      if (riseShape_ === "linear") return sumEventsWithLinearRamp(events, tn, riseTimeTr);
      return sumEventsWithRise(events, tn, riseTimeTr);
    };

    return {
      launchedWaves,
      activeWaves,
      VS: sumNode(bounce.srcEvents),
      VL: sumNode(bounce.loadEvents),
    };
  }

  // ---- Method of Characteristics time-domain simulator ----
  //
  // model.blocks = [
  //   { type: 'tl', Z0, tau },   // T-line segment; tau = one-way delay in τ_d units
  //   { type: 'R',  value },      // shunt resistor to ground (Ω)
  //   { type: 'C',  value },      // shunt capacitor (F, SI) — Step 3
  //   { type: 'L',  value },      // shunt inductor  (H, SI) — Step 4
  //   { type: 'short' },          // shunt short to ground
  // ]
  // model.terminal = { type: 'R'|'open'|'short'|'C'|'L', value }
  // model.Rg, model.Vg, model.riseTimeTr, model.riseShape: unchanged
  //
  // Consecutive shunt blocks (no 'tl' between them) share one electrical node;
  // their admittances accumulate in the node equation.
  // Consecutive 'tl' blocks share a pure impedance-mismatch boundary (no shunt).
  //
  // opts: { tEnd, oversample }
  //   tEnd      — simulation end time in τ_d units (default 10)
  //   oversample — delay steps per shortest segment (default 8)
  //
  // Returns { dt, nSteps, tGrid, segs, nDelay, nodeV, vPlusHist, vMinusHist }
  //   dt          — time step in τ_d units
  //   tGrid       — Float64Array of length nSteps: tGrid[k] = k * dt
  //   segs        — [{Z0, tau}] after parsing (T-line segments only)
  //   nDelay      — Int32Array: nDelay[i] = delay of segment i in time steps
  //   nodeV       — (N+1) × nSteps: nodeV[j][k] = voltage at node j at step k
  //   vPlusHist   — N × nSteps: V⁺ written into left end of segment i at step k
  //   vMinusHist  — N × nSteps: V⁻ written into right end of segment i at step k
  //
  // Spatial query — see voltageAt().
  function simulateTimeDomain(model, opts) {
    const { Vg: Vg_final, Rg } = model;
    const tEnd_      = (opts && opts.tEnd)      != null ? opts.tEnd      : 10;
    const oversample = (opts && opts.oversample) != null ? opts.oversample : 8;

    // --- parse model.blocks into T-line segments + per-node shunt lists ---
    // nodeShunts[j] are the shunt elements sitting at node j.
    // Node j sits between segment j-1 (on the left) and segment j (on the right).
    // Node 0 = source side; node N = load side.
    const segs       = [];   // {Z0, tau}
    const nodeShunts = [];   // nodeShunts[j] = [{type, value?}]
    let   pending    = [];   // shunts accumulating until the next 'tl' block
    for (const blk of model.blocks) {
      if (blk.type === 'tl') {
        nodeShunts.push(pending);
        pending = [];
        segs.push({ Z0: blk.Z0, tau: blk.tau });
      } else {
        pending.push(blk);
      }
    }
    nodeShunts.push(pending);   // shunts at load-side node (right of last T-line)
    const N = segs.length;
    if (N === 0) throw new Error('simulateTimeDomain: no tl blocks in model.blocks');

    // --- time step and per-segment delay lengths ---
    // dt chosen so the shortest segment spans exactly `oversample` steps.
    const minTau   = segs.reduce((m, s) => Math.min(m, s.tau), Infinity);
    const dt       = minTau / oversample;
    const nSteps   = Math.ceil(tEnd_ / dt) + 1;
    // dt_SI: real-time duration of one time step (seconds).
    // model.tau_d must be provided (SI seconds per tau unit) when any shunt C or L
    // is present; for purely resistive circuits it is unused (default 1 is harmless).
    const tau_d_SI = model.tau_d != null ? model.tau_d : 1;
    const dt_SI    = dt * tau_d_SI;
    // Each segment rounded to the nearest integer number of steps (≥ 1).
    const nDelay = new Int32Array(N);
    for (let i = 0; i < N; i++) nDelay[i] = Math.max(1, Math.round(segs[i].tau / dt));

    // --- ring buffers ---
    // delayFwd[i]: V⁺ entering left end of segment i.
    //   Written at step k, read (= arrives at right end) at step k + nDelay[i].
    // delayBwd[i]: V⁻ entering right end of segment i.
    //   Written at step k, arrives at left end at step k + nDelay[i].
    // Both use the same ring-buffer idiom: fwdIdx[i] / bwdIdx[i] is the current
    // write position; that same position holds the oldest (= arriving) value.
    const delayFwd = Array.from({ length: N }, (_, i) => new Float64Array(nDelay[i]));
    const delayBwd = Array.from({ length: N }, (_, i) => new Float64Array(nDelay[i]));
    const fwdIdx   = new Int32Array(N);
    const bwdIdx   = new Int32Array(N);

    // --- reactive state per node (C: vCap; L: iL) ---
    // Each element of nodeShuntState[j] is a mutable object parallel to nodeShunts[j].
    const nodeShuntState = nodeShunts.map(shunts =>
      shunts.map(sh => ({ vCap: 0, iL: 0 }))
    );
    const termState = { vCap: 0, iL: 0 };

    // --- output storage ---
    const nodeV      = Array.from({ length: N + 1 }, () => new Float64Array(nSteps));
    const vPlusHist  = Array.from({ length: N },     () => new Float64Array(nSteps));
    const vMinusHist = Array.from({ length: N },     () => new Float64Array(nSteps));

    // --- source waveform ---
    const riseTimeTr_ = model.riseTimeTr || 0;
    const riseMode_   = model.riseShape  || 'step';
    function vgAt(tn) {
      if (riseTimeTr_ <= 0) return tn > -1e-12 ? Vg_final : 0;
      if (riseMode_ === 'linear') return Vg_final * riseShapeLinear(tn, riseTimeTr_);
      return Vg_final * riseShape(tn, riseTimeTr_);
    }

    // --- helper: accumulate shunt contributions into node equation ---
    // Each shunt contributes conductance G and Norton current I_src such that:
    //   V * (G_base + ΣG) = I_base + ΣI_src
    // where G_base / I_base come from the Thevenin equivalent(s) of adjacent T-lines.
    // Returns true if the node is hard-shorted to ground.
    function accumulateShunts(shunts, states, Vnode_prev_fn) {
      let G = 0, I = 0, shorted = false;
      for (let s = 0; s < shunts.length; s++) {
        const sh = shunts[s], st = states[s];
        if (sh.type === 'R') {
          G += 1 / sh.value;
        } else if (sh.type === 'C') {
          const Gc = sh.value / dt_SI;
          G += Gc;
          I += Gc * st.vCap;
        } else if (sh.type === 'L') {
          const Gl = dt_SI / sh.value;
          G += Gl;
          I -= st.iL;            // see derivation: inductor history current subtracts
        } else if (sh.type === 'short') {
          shorted = true;
        }
      }
      return { G, I, shorted };
    }

    function updateShuntStates(shunts, states, V) {
      for (let s = 0; s < shunts.length; s++) {
        const sh = shunts[s], st = states[s];
        if (sh.type === 'C') st.vCap = V;
        if (sh.type === 'L') st.iL  += (dt_SI / sh.value) * V;
      }
    }

    // --- main simulation loop ---
    const vArr = {
      right: new Float64Array(N),  // V⁺ arriving at right end of each segment
      left:  new Float64Array(N),  // V⁻ arriving at left end of each segment
    };
    const Vnodes = new Float64Array(N + 1);

    for (let step = 0; step < nSteps; step++) {
      const tn = step * dt;

      // 1. Read arriving waves from delay lines.
      for (let i = 0; i < N; i++) {
        vArr.right[i] = delayFwd[i][fwdIdx[i]];
        vArr.left[i]  = delayBwd[i][bwdIdx[i]];
      }

      // 2. Solve node voltages.

      // Source node (j=0): Thevenin from {Vg, Rg} on left; segment 0 on right.
      // KCL: (vg - V)/Rg + (2·aR - V)/Z0 = 0  → V = (vg·Z0 + 2·aR·Rg)/(Rg + Z0)
      // (Source-node shunts, if any, are rare but handled.)
      {
        const vg = vgAt(tn);
        const aR = vArr.left[0];      // V⁻ arriving at left end of seg 0
        const Z0 = segs[0].Z0;
        let G = 1 / Rg + 1 / Z0;
        let I = vg / Rg + 2 * aR / Z0;
        const { G: Gs, I: Is, shorted } = accumulateShunts(
          nodeShunts[0], nodeShuntState[0]);
        if (shorted) {
          Vnodes[0] = 0;
        } else {
          Vnodes[0] = (I + Is) / (G + Gs);
        }
        updateShuntStates(nodeShunts[0], nodeShuntState[0], Vnodes[0]);
      }

      // Internal nodes (j = 1 … N-1).
      for (let j = 1; j < N; j++) {
        const ZL = segs[j - 1].Z0;
        const ZR = segs[j].Z0;
        const aL = vArr.right[j - 1];  // V⁺ arriving at right end of seg j-1
        const aR = vArr.left[j];        // V⁻ arriving at left end of seg j
        let G = 1 / ZL + 1 / ZR;
        let I = 2 * aL / ZL + 2 * aR / ZR;
        const { G: Gs, I: Is, shorted } = accumulateShunts(
          nodeShunts[j], nodeShuntState[j]);
        if (shorted) {
          Vnodes[j] = 0;
        } else {
          Vnodes[j] = (I + Is) / (G + Gs);
        }
        updateShuntStates(nodeShunts[j], nodeShuntState[j], Vnodes[j]);
      }

      // Load node (j=N): segment N-1 on left; terminal on right.
      {
        const aL   = vArr.right[N - 1];   // V⁺ arriving at right end of last seg
        const Z0   = segs[N - 1].Z0;
        const term = model.terminal;
        // First apply any mid-node shunts at the load-side node (nodeShunts[N]).
        let G = 1 / Z0;
        let I = 2 * aL / Z0;
        const { G: Gs, I: Is, shorted: sh0 } = accumulateShunts(
          nodeShunts[N], nodeShuntState[N]);
        if (sh0) {
          Vnodes[N] = 0;
        } else {
          // Then add the terminal element.
          let Gt = 0, It = 0, shorted = false;
          if (term.type === 'open') {
            // no terminal conductance
          } else if (term.type === 'short') {
            shorted = true;
          } else if (term.type === 'R') {
            Gt = 1 / term.value;
          } else if (term.type === 'C') {
            Gt = term.value / dt_SI;
            It = Gt * termState.vCap;
          } else if (term.type === 'L') {
            Gt = dt_SI / term.value;
            It = -termState.iL;
          }
          Vnodes[N] = shorted ? 0 : (I + Is + It) / (G + Gs + Gt);
        }
        updateShuntStates(nodeShunts[N], nodeShuntState[N], Vnodes[N]);
        if (term.type === 'C') termState.vCap = Vnodes[N];
        if (term.type === 'L') termState.iL  += (dt_SI / term.value) * Vnodes[N];
      }

      // 3. Compute outgoing waves and write into delay lines.
      // At node j: V = V⁺_out + V⁻_in  →  V⁺_out = V - V⁻_in  (outgoing into seg j rightward)
      // At node j: V = V⁺_in + V⁻_out  →  V⁻_out = V - V⁺_in  (outgoing into seg j-1 leftward)
      for (let i = 0; i < N; i++) {
        const vOutFwd = Vnodes[i]     - vArr.left[i];   // V⁺ leaving node i  into seg i
        const vOutBwd = Vnodes[i + 1] - vArr.right[i];  // V⁻ leaving node i+1 into seg i
        vPlusHist[i][step]  = vOutFwd;
        vMinusHist[i][step] = vOutBwd;
        delayFwd[i][fwdIdx[i]] = vOutFwd;
        delayBwd[i][bwdIdx[i]] = vOutBwd;
        fwdIdx[i] = (fwdIdx[i] + 1) % nDelay[i];
        bwdIdx[i] = (bwdIdx[i] + 1) % nDelay[i];
      }

      // 4. Record node voltages.
      for (let j = 0; j <= N; j++) nodeV[j][step] = Vnodes[j];
    }

    const tGrid = Float64Array.from({ length: nSteps }, (_, k) => k * dt);
    return { dt, nSteps, tGrid, segs, nDelay, nodeV, vPlusHist, vMinusHist };
  }

  // ---- Spatial voltage query for simulateTimeDomain output ----
  //
  // Returns V(z, tIdx) by reading from vPlusHist / vMinusHist at the appropriate
  // delay offsets.  z ∈ [0, 1] is normalized by total T-line delay (sum of taus).
  //
  // Within segment i, at fractional position f ∈ [0, 1]:
  //   V⁺ that is at position f at step tIdx was written into the segment
  //   kFwd = round(f * D) steps ago  → look up vPlusHist[i][tIdx - kFwd]
  //   V⁻ that is at position f at step tIdx was written into the segment
  //   kBwd = round((1-f) * D) steps ago → look up vMinusHist[i][tIdx - kBwd]
  function voltageAt(z, tIdx, sim) {
    const { segs, nDelay, vPlusHist, vMinusHist, nodeV } = sim;
    const N = segs.length;
    const totalTau = segs.reduce((s, seg) => s + seg.tau, 0);
    let cumZ = 0;
    for (let i = 0; i < N; i++) {
      const segZ = segs[i].tau / totalTau;
      if (z <= cumZ + segZ + 1e-9) {
        const f    = Math.min(1, Math.max(0, (z - cumZ) / segZ));
        const D    = nDelay[i];
        const kFwd = Math.round(f * D);
        const kBwd = D - kFwd;
        const tF   = tIdx - kFwd;
        const tB   = tIdx - kBwd;
        const vP   = tF >= 0 ? vPlusHist[i][tF]  : 0;
        const vM   = tB >= 0 ? vMinusHist[i][tB] : 0;
        return vP + vM;
      }
      cumZ += segZ;
    }
    // z at or past the load end: return load node voltage
    return nodeV[N][tIdx];
  }

  // ---- Convenience: convert old-style model (model.segments, model.RL) to blocks ----
  // This lets existing callers construct a blocks-based model without rewriting app.js.
  function segmentsToBlocks(model) {
    const N = model.segments.length;
    const tau = 1 / N;          // each segment = 1/N of total normalized delay
    const blocks = model.segments.map(seg => ({ type: 'tl', Z0: seg.Z0, tau }));
    const RL = model.RL;
    const terminal = isFinite(RL)
      ? { type: 'R', value: RL }
      : { type: 'open' };
    return { ...model, blocks, terminal };
  }

  return {
    BOUNCE_EPS,
    MAX_BOUNCES,
    computeWaveParams,
    buildBounceSeries,
    sumEventsAtTime,
    sumEventsWithRise,
    sumEventsWithLinearRamp,
    riseShape,
    riseShapeLinear,
    waveVoltageAt,
    totalVoltageAt,
    segmentForZ,
    computeDynamicState,
    simulateTimeDomain,
    voltageAt,
    segmentsToBlocks,
  };
})();