// Smoke tests for TLPhysics pure functions.
// Requires Node 18+.  Run with:  node --test physics.test.js
"use strict";

const { test }       = require("node:test");
const assert         = require("node:assert/strict");
const { readFileSync, existsSync } = require("node:fs");
const path           = require("path");

// utils.js calls window.devicePixelRatio inside getDPR(); stub it so the IIFE
// doesn't throw during load.  Physics functions never call getDPR.
globalThis.window = { devicePixelRatio: 1 };

// eval() scopes `const` to the eval block, so we append an explicit globalThis
// assignment after each IIFE declaration to make the namespace visible.
eval(readFileSync(path.join(__dirname, "utils.js"),   "utf8") + "\nglobalThis.TLUtils = TLUtils;");
eval(readFileSync(path.join(__dirname, "physics.js"), "utf8") + "\nglobalThis.TLPhysics = TLPhysics;");

const {
  computeWaveParams,
  buildBounceSeries,
  sumEventsAtTime,
  sumEventsWithLinearRamp,
  riseShape,
  riseShapeLinear,
  totalVoltageAt,
  computeDynamicState,
} = globalThis.TLPhysics;

// Helper: assert two numbers are within `tol` of each other.
function near(a, b, msg, tol = 1e-9) {
  assert.ok(
    Math.abs(a - b) <= tol,
    `${msg}: expected ~${b}, got ${a} (diff ${(a - b).toExponential(2)})`
  );
}

// Build a minimal single-segment model and run buildBounceSeries.
function makeModel(Vg, Rg, Z0, RL, reflectTol = 0.001) {
  const model = { Vg, Rg, RL, segments: [{ Z0 }], reflectTol };
  const waves = computeWaveParams(model);
  const bounce = buildBounceSeries(model, waves);
  return { model, waves, bounce };
}

// ────────────────────────────────────────────────────────────────────────────
// computeWaveParams
// ────────────────────────────────────────────────────────────────────────────

test("computeWaveParams: matched line (Rg=Z0=RL=50)", () => {
  const m = { Vg: 1, Rg: 50, RL: 50, segments: [{ Z0: 50 }] };
  const { V1, gL, gS } = computeWaveParams(m);
  near(V1, 0.5,  "V1");
  near(gL, 0,    "gL");
  near(gS, 0,    "gS");
});

test("computeWaveParams: open-circuit load (RL=Infinity)", () => {
  const m = { Vg: 1, Rg: 50, RL: Infinity, segments: [{ Z0: 50 }] };
  const { V1, gL, gS } = computeWaveParams(m);
  near(V1, 0.5,  "V1");
  near(gL, 1,    "gL");
  near(gS, 0,    "gS");
});

test("computeWaveParams: short-circuit load (RL=0)", () => {
  const m = { Vg: 1, Rg: 50, RL: 0, segments: [{ Z0: 50 }] };
  const { V1, gL, gS } = computeWaveParams(m);
  near(V1, 0.5,  "V1");
  near(gL, -1,   "gL");
  near(gS, 0,    "gS");
});

test("computeWaveParams: source Γ with Rg mismatch", () => {
  // Rg=100, Z0=50 → gS = (100-50)/(100+50) = 50/150 = 1/3
  const m = { Vg: 1, Rg: 100, RL: 50, segments: [{ Z0: 50 }] };
  const { gS } = computeWaveParams(m);
  near(gS, 1 / 3, "gS", 1e-9);
});

// ────────────────────────────────────────────────────────────────────────────
// buildBounceSeries: wave count and first-wave correctness
// ────────────────────────────────────────────────────────────────────────────

test("buildBounceSeries: matched line produces exactly 1 wave", () => {
  const { bounce } = makeModel(1, 50, 50, 50);
  assert.equal(bounce.series.length, 1, "exactly one wave packet");
  assert.equal(bounce.series[0].dir, +1, "rightward");
  near(bounce.series[0].A, 0.5, "amplitude = V1 = 0.5");
});

test("buildBounceSeries: first wave geometry", () => {
  const { bounce } = makeModel(1, 50, 50, 50);
  const w = bounce.series[0];
  near(w.zStart,  0, "zStart = 0");
  near(w.zEnd,    1, "zEnd = 1");
  near(w.tBorn,   0, "tBorn = 0");
  near(w.tDie,    1, "tDie = 1");
});

test("buildBounceSeries: open-circuit — load event dV = 2·V1 at t=1", () => {
  // gL = 1 → (1+gL)·A = 1.0.  gS=0 → no further reflections → exactly 1 load event.
  const { bounce } = makeModel(1, 50, 50, Infinity);
  assert.equal(bounce.loadEvents.length, 1, "one load event");
  near(bounce.loadEvents[0].t,  1,   "event time = 1τ");
  near(bounce.loadEvents[0].dV, 1.0, "dV = 1.0 (= 2·0.5)");
});

test("buildBounceSeries: short-circuit — load event dV = 0 at t=1", () => {
  // gL = -1 → (1+gL)·A = 0
  const { bounce } = makeModel(1, 50, 50, 0);
  assert.equal(bounce.loadEvents.length, 1, "one load event");
  near(bounce.loadEvents[0].dV, 0, "dV = 0 at short circuit");
});

// ────────────────────────────────────────────────────────────────────────────
// DC steady state: VL(t→∞) = Vg · RL / (Rg + RL)
// ────────────────────────────────────────────────────────────────────────────

function checkDCSteadyState(label, Vg, Rg, Z0, RL, tol = 1e-3) {
  test(`DC steady state: ${label}`, () => {
    const { bounce } = makeModel(Vg, Rg, Z0, RL, 0.001);
    const VL_dc      = isFinite(RL) ? Vg * RL / (Rg + RL) : Vg;
    const VL_actual  = sumEventsAtTime(bounce.loadEvents, bounce.tEnd + 10);
    near(VL_actual, VL_dc, "VL", tol);
  });
}

// Matched source, various loads
checkDCSteadyState("matched (50/50/50)",           1, 50,  50,  50);
checkDCSteadyState("open circuit (50/50/∞)",        1, 50,  50,  Infinity);
checkDCSteadyState("short circuit (50/50/0)",      1, 50,  50,  0);
// Mismatched source
checkDCSteadyState("Rg=100 Z0=50 RL=150, Vg=5",   5, 100, 50,  150, 1e-3);
checkDCSteadyState("Rg=25  Z0=75 RL=200, Vg=3.3", 3.3, 25, 75, 200, 1e-3);

// ────────────────────────────────────────────────────────────────────────────
// Multi-segment: two identical segments should match single segment
// ────────────────────────────────────────────────────────────────────────────

test("multi-segment: two identical Z0 segments == single segment", () => {
  const single = makeModel(1, 50, 50, 100, 0.001);
  const model2 = { Vg: 1, Rg: 50, RL: 100, segments: [{ Z0: 50 }, { Z0: 50 }], reflectTol: 0.001 };
  const waves2 = computeWaveParams(model2);
  const bounce2 = buildBounceSeries(model2, waves2);
  const t = single.bounce.tEnd + 10;
  const VL_single = sumEventsAtTime(single.bounce.loadEvents, t);
  const VL_double = sumEventsAtTime(bounce2.loadEvents, t);
  near(VL_single, VL_double, "VL matches", 1e-3);
});

test("multi-segment: two different Z0 — DC still converges to Vg·RL/(Rg+RL)", () => {
  const model = { Vg: 1, Rg: 50, RL: 100, segments: [{ Z0: 50 }, { Z0: 75 }], reflectTol: 0.001 };
  const waves  = computeWaveParams(model);
  const bounce = buildBounceSeries(model, waves);
  const VL_dc  = 1 * 100 / (50 + 100);  // ≈ 0.6667
  const VL_actual = sumEventsAtTime(bounce.loadEvents, bounce.tEnd + 10);
  near(VL_actual, VL_dc, "VL", 1e-3);
});

// ────────────────────────────────────────────────────────────────────────────
// sumEventsAtTime
// ────────────────────────────────────────────────────────────────────────────

test("sumEventsAtTime: empty events → 0", () => {
  near(sumEventsAtTime([], 999), 0, "empty");
});

test("sumEventsAtTime: single event at t=0 seen immediately", () => {
  near(sumEventsAtTime([{ t: 0, dV: 0.5 }], 0), 0.5, "at t=0");
});

test("sumEventsAtTime: event at t=1 not seen before it", () => {
  near(sumEventsAtTime([{ t: 1, dV: 0.5 }], 0.5), 0, "before event");
  near(sumEventsAtTime([{ t: 1, dV: 0.5 }], 1.0), 0.5, "at event time");
  near(sumEventsAtTime([{ t: 1, dV: 0.5 }], 2.0), 0.5, "after event time");
});

test("sumEventsAtTime: accumulates multiple events", () => {
  const evts = [{ t: 0, dV: 0.5 }, { t: 1, dV: 0.25 }, { t: 2, dV: 0.1 }];
  near(sumEventsAtTime(evts, 1.5), 0.75, "after first two events");
  near(sumEventsAtTime(evts, 2.0), 0.85, "after all events");
});

// ────────────────────────────────────────────────────────────────────────────
// riseShape
// ────────────────────────────────────────────────────────────────────────────

test("riseShape: tr=0 (step) — hard step at dt=0", () => {
  near(riseShape(0,     0), 0, "dt=0 tr=0");
  near(riseShape(-1,    0), 0, "dt<0 tr=0");
  near(riseShape(0.001, 0), 1, "tiny dt, tr=0");
  near(riseShape(100,   0), 1, "large dt, tr=0");
});

test("riseShape: tr>0 — causal shifted erf over [0, tr]", () => {
  // Causal: dt <= 0 → exactly 0
  near(riseShape(0, 1), 0, "dt=0 → 0", 1e-3);
  near(riseShape(-1, 1), 0, "dt<0 → 0");
  // dt = tr → essentially 1
  near(riseShape(1, 1), 1, "dt=tr → 1", 1e-3);
  // Midpoint dt = tr/2 → 0.5 (erf(0) = 0)
  near(riseShape(0.5, 1), 0.5, "dt=tr/2 → 0.5", 1e-7);
  // Monotonically increasing
  const v1 = riseShape(0.2, 1);
  const v2 = riseShape(0.5, 1);
  const v3 = riseShape(0.8, 1);
  assert(v1 < v2, `monotonic: r(0.2)=${v1} < r(0.5)=${v2}`);
  assert(v2 < v3, `monotonic: r(0.5)=${v2} < r(0.8)=${v3}`);
  // Clamped to 1 for dt > tr
  near(riseShape(2, 1), 1, "dt>tr → 1");
});

// ────────────────────────────────────────────────────────────────────────────
// riseShapeLinear
// ────────────────────────────────────────────────────────────────────────────

test("riseShapeLinear: dt≤0 → 0", () => {
  near(riseShapeLinear(0,   1), 0, "dt=0");
  near(riseShapeLinear(-1,  1), 0, "dt<0");
});

test("riseShapeLinear: tr=0 (step) → 1 for dt>0", () => {
  near(riseShapeLinear(0.001, 0), 1, "tiny dt, tr=0");
  near(riseShapeLinear(100,   0), 1, "large dt, tr=0");
});

test("riseShapeLinear: linear ramp 0→1 over tr", () => {
  near(riseShapeLinear(0.1, 1), 0.1,  "10% of tr",  1e-12);
  near(riseShapeLinear(0.5, 1), 0.5,  "50% of tr",  1e-12);
  near(riseShapeLinear(0.9, 1), 0.9,  "90% of tr",  1e-12);
  near(riseShapeLinear(1.0, 1), 1.0,  "exactly tr", 1e-12);
  near(riseShapeLinear(2.0, 1), 1.0,  "past tr",    1e-12);
});

test("riseShapeLinear: 10–90% rise time is exactly 0.8·tr", () => {
  const tr = 0.3;
  near(riseShapeLinear(0.1 * tr, tr), 0.1, "10%", 1e-12);
  near(riseShapeLinear(0.9 * tr, tr), 0.9, "90%", 1e-12);
});

// ────────────────────────────────────────────────────────────────────────────
// totalVoltageAt: no spike at segment boundaries (multi-segment plotting bug)
//
// In smooth mode, drawSampledWave evaluates totalVoltageAt at exact boundary
// z values (e.g. z=0.5 for 2 segments).  Without segment filtering, the parent
// wave (dir=+1, zEnd=boundary) and the transmitted wave (dir=+1, zStart=boundary)
// both pass the range check, doubling the voltage at the boundary.
// ────────────────────────────────────────────────────────────────────────────

test("totalVoltageAt: no spike at segment boundary (2 segments, step mode)", () => {
  // Z0=40 → Z0=60, Rg=50, RL=∞.
  // Γ_bound = (60-40)/(60+40) = 0.2.  V1 = 1 * 40/(50+40) = 4/9.
  // At tNorm=0.6 the boundary crossing (tNorm=0.5) is complete:
  //   seg-0 has incident(4/9) + reflected(0.2·4/9), seg-1 has transmitted(1.2·4/9).
  // Correct V at z=0.5 = (1+0.2)·(4/9) = 4.8/9 = 8/15.
  // Buggy code (no segIdx filter) would give 2·(8/15) = 16/15 — a visible spike.
  const model  = { Vg: 1, Rg: 50, RL: Infinity,
                   segments: [{ Z0: 40 }, { Z0: 60 }], reflectTol: 0.001 };
  const waves  = computeWaveParams(model);
  const bounce = buildBounceSeries(model, waves);
  const dyn    = computeDynamicState(0.6, bounce);
  const N      = model.segments.length;
  const eps    = 1e-4;

  const vL = totalVoltageAt(0.5 - eps, dyn.launchedWaves, 0, "step", N);
  const vB = totalVoltageAt(0.5,       dyn.launchedWaves, 0, "step", N);
  const vR = totalVoltageAt(0.5 + eps, dyn.launchedWaves, 0, "step", N);

  // Correct value: (1 + 0.2) * (4/9) = 8/15
  near(vB, 8 / 15, "V at boundary = (1+Γ)·V1", 1e-6);
  near(vB, vL, "no spike: V(boundary) ≈ V(boundary-ε)", 1e-6);
  near(vB, vR, "no spike: V(boundary) ≈ V(boundary+ε)", 1e-6);
});

test("totalVoltageAt: no spike at boundaries with 4 segments", () => {
  // 4 equal segments: boundaries at 0.25, 0.5, 0.75.
  // tNorm=0.3 puts the front in seg 1; all three boundary z-values should be spike-free.
  const model  = { Vg: 1, Rg: 50, RL: 100,
                   segments: [{ Z0: 50 }, { Z0: 75 }, { Z0: 50 }, { Z0: 75 }],
                   reflectTol: 0.001 };
  const waves  = computeWaveParams(model);
  const bounce = buildBounceSeries(model, waves);
  const N      = model.segments.length;
  const eps    = 1e-4;

  for (const tNorm of [0.3, 0.6, 1.2, 2.0]) {
    const dyn = computeDynamicState(tNorm, bounce);
    for (const zB of [0.25, 0.5, 0.75]) {
      const vL = totalVoltageAt(zB - eps, dyn.launchedWaves, 0, "step", N);
      const vB = totalVoltageAt(zB,       dyn.launchedWaves, 0, "step", N);
      const vR = totalVoltageAt(zB + eps, dyn.launchedWaves, 0, "step", N);
      // No spike: boundary value must lie between its two neighbours (within tolerance).
      const lo = Math.min(vL, vR) - 1e-6;
      const hi = Math.max(vL, vR) + 1e-6;
      assert.ok(vB >= lo && vB <= hi,
        `spike at z=${zB}, tNorm=${tNorm}: V=${vB}, neighbours ${vL}..${vR}`);
    }
  }
});

// ────────────────────────────────────────────────────────────────────────────

// ────────────────────────────────────────────────────────────────────────────
// SPICE golden-reference comparisons
//
// Common circuit (all netlists in sim/):
//   Z0=50Ω  (LTRA l=50nH/m c=20pF/m len=1m → τ_d=1ns)
//   Rg=10Ω  (Rs1),  Vg=1V  PULSE(0 1 0 100p 100p 1u 2u) → TR=100ps=0.1·τ_d
//
// TSV columns: t[s]  v(a)  t[s]  v(b)
// Run simulations: cd sim && sh run.sh   (or ngspice <netlist>.sp individually)
// ────────────────────────────────────────────────────────────────────────────

// Compare every row of a SPICE TSV against the physics model.
// Skips automatically if the TSV has not been generated yet.
function spiceCompare(label, tsvFile, RL, TOL = 1e-4) {
  const TAU_D   = 1e-9;
  const TR_NORM = 0.1;
  const tsvPath = path.join(__dirname, "sim", tsvFile);
  const skip    = !existsSync(tsvPath);

  test(`SPICE comparison: ${label}`, { skip: skip ? "TSV not found — run: cd sim && ngspice" : false }, () => {
    const model  = { Vg: 1, Rg: 10, RL, segments: [{ Z0: 50 }], reflectTol: 0.001 };
    const waves  = computeWaveParams(model);
    const bounce = buildBounceSeries(model, waves);

    const rows = readFileSync(tsvPath, "utf8").trim().split("\n")
      .map((line) => line.trim().split(/\s+/).map(Number))
      .filter((cols) => cols.length >= 4);
    assert.ok(rows.length > 10, "TSV should have many rows");

    for (const [t_s, va_spice, , vb_spice] of rows) {
      const tn       = t_s / TAU_D;
      const va_model = sumEventsWithLinearRamp(bounce.srcEvents,  tn, TR_NORM);
      const vb_model = sumEventsWithLinearRamp(bounce.loadEvents, tn, TR_NORM);
      near(va_model, va_spice, `v(a) at t=${t_s.toExponential(3)}s`, TOL);
      near(vb_model, vb_spice, `v(b) at t=${t_s.toExponential(3)}s`, TOL);
    }
  });
}

// RL=∞  (open circuit)    ΓL=+1,   ΓS=−2/3
spiceCompare("open circuit  (tline-oc.sp)",       "results-tline-oc.tsv", Infinity);
// RL=0   (short circuit)   ΓL=−1,   ΓS=−2/3
spiceCompare("short circuit (tline-sc.sp)",       "results-tline-sc.tsv", 0);
// RL=100 (resistive load)  ΓL=+1/3, ΓS=−2/3
spiceCompare("resistive load 100Ω (tline-rl.sp)", "results-tline-rl.tsv", 100);

// ────────────────────────────────────────────────────────────────────────────
// simulateTimeDomain: validate against buildBounceSeries (resistive cases)
//
// Strategy: compare node voltages at the source (node 0) and load (node N) at
// several tNorm values.  We sample at half-integer multiples of τ_d so that
// we land between bounce events, where both models agree on a constant value.
// With integer delay steps the MoC simulation is exact for resistive circuits.
// ────────────────────────────────────────────────────────────────────────────

const { simulateTimeDomain, voltageAt, segmentsToBlocks } = globalThis.TLPhysics;

// Helper: build a blocks-format model from old-style params and run the sim.
function makeSim(Vg, Rg, Z0, RL, tEnd = 12, oversample = 16) {
  const oldModel = { Vg, Rg, RL, segments: [{ Z0 }], reflectTol: 0.001 };
  const newModel = segmentsToBlocks(oldModel);
  return { sim: simulateTimeDomain(newModel, { tEnd, oversample }), oldModel };
}

// Compare VS (node 0) and VL (node N) between MoC sim and bounce series.
function checkSimVsBounce(label, Vg, Rg, Z0, RL, tNorms, tol = 1e-9) {
  test(`simulateTimeDomain vs bounce: ${label}`, () => {
    const { sim, oldModel } = makeSim(Vg, Rg, Z0, RL);
    const waves  = computeWaveParams(oldModel);
    const bounce = buildBounceSeries(oldModel, waves);
    const { dt, nSteps, nodeV } = sim;
    for (const tn of tNorms) {
      const step = Math.round(tn / dt);
      if (step >= nSteps) continue;
      near(nodeV[0][step], sumEventsAtTime(bounce.srcEvents,  tn), `${label} VS @ tNorm=${tn}`, tol);
      near(nodeV[1][step], sumEventsAtTime(bounce.loadEvents, tn), `${label} VL @ tNorm=${tn}`, tol);
    }
  });
}

// Sample at half-integer τ values: between bounce events both models are constant.
const tSamples = [0.5, 1.5, 2.5, 3.5, 4.5, 5.5];

checkSimVsBounce("matched (50/50/50)",       1, 50,  50, 50,       tSamples);
checkSimVsBounce("open circuit (50/50/∞)",    1, 50,  50, Infinity, tSamples);
checkSimVsBounce("short circuit (50/50/0)",   1, 50,  50, 0,        tSamples);
checkSimVsBounce("mismatched Rg=100 RL=200",  1, 100, 50, 200,      tSamples);
checkSimVsBounce("Vg=5 Rg=20 Z0=50 RL=100",  5, 20,  50, 100,      tSamples);

test("simulateTimeDomain: 2 equal segments match bounce series VL", () => {
  const oldModel = { Vg: 1, Rg: 50, RL: 100, segments: [{ Z0: 50 }, { Z0: 50 }], reflectTol: 0.001 };
  const newModel = segmentsToBlocks(oldModel);
  const sim      = simulateTimeDomain(newModel, { tEnd: 12, oversample: 16 });
  const bounce   = buildBounceSeries(oldModel, computeWaveParams(oldModel));
  const { dt, nSteps, nodeV } = sim;
  for (const tn of tSamples) {
    const step = Math.round(tn / dt);
    if (step >= nSteps) continue;
    near(nodeV[2][step], sumEventsAtTime(bounce.loadEvents, tn), `VL @ tNorm=${tn}`, 1e-9);
  }
});

test("simulateTimeDomain: 2 segments Z0=50/75, DC = Vg·RL/(Rg+RL)", () => {
  const oldModel = { Vg: 1, Rg: 50, RL: 100, segments: [{ Z0: 50 }, { Z0: 75 }], reflectTol: 0.001 };
  const newModel = segmentsToBlocks(oldModel);
  const sim      = simulateTimeDomain(newModel, { tEnd: 20, oversample: 8 });
  const step     = Math.min(Math.round(18 / sim.dt), sim.nSteps - 1);
  near(sim.nodeV[2][step], 1 * 100 / (50 + 100), "VL DC steady state", 1e-3);
});

test("voltageAt: z=0 and z=1 match source/load node voltages", () => {
  const { sim } = makeSim(1, 50, 50, Infinity);
  const { dt, nSteps } = sim;
  for (const tn of [0.5, 1.5, 2.5]) {
    const step = Math.round(tn / dt);
    if (step >= nSteps) continue;
    near(voltageAt(0, step, sim), sim.nodeV[0][step], `z=0 @ tNorm=${tn}`, 1e-9);
    near(voltageAt(1, step, sim), sim.nodeV[1][step], `z=1 @ tNorm=${tn}`, 1e-9);
  }
});

test("voltageAt: midpoint of matched line = V1 once front has passed (tNorm=0.7)", () => {
  // Matched line Rg=Z0=RL=50: single rightward wave V1=0.5, no reflections.
  // At tNorm=0.7 the front is past z=0.5, so V(0.5) = V1 = 0.5.
  const { sim } = makeSim(1, 50, 50, 50, 12, 16);
  const step = Math.round(0.7 / sim.dt);
  near(voltageAt(0.5, step, sim), 0.5, "V(0.5) = 0.5", 1e-9);
});

test("voltageAt: midpoint = 0 before wave arrives (tNorm=0.3)", () => {
  const { sim } = makeSim(1, 50, 50, 50, 12, 16);
  const step = Math.round(0.3 / sim.dt);
  near(voltageAt(0.5, step, sim), 0, "V(0.5) = 0 before wave", 1e-9);
});

// ────────────────────────────────────────────────────────────────────────────
// SPICE comparison: shunt-R at midpoint between two T-line segments
//
// Circuit: Rg=10Ω → T1(Z0=50,τ=1ns) → Rshunt=100Ω → T2(Z0=50,τ=1ns) → open
// Netlist: sim/tline-shunt-r.sp   TSV columns: t v(a) t v(b) t v(c)
//
// Tolerance is 1e-2 (vs 1e-4 for the bounce-series SPICE tests) because the
// MoC simulator is discrete-time: linear interpolation between steps introduces
// a bounded error at ramp knees proportional to dt × slope ≈ 8 mV at oversample=256.
// This confirms the physics is correct, not floating-point exact.
// ────────────────────────────────────────────────────────────────────────────
// SPICE comparison: shunt-C at midpoint between two T-line segments
//
// Circuit: Rg=10Ω → T1(Z0=50,τ=1ns) → C1=2pF → T2(Z0=50,τ=1ns) → open
// Netlist: sim/tline-shunt-c.sp   TSV columns: t v(a) t v(b) t v(c)
//
// The 2pF cap creates an RC time constant τ_RC = C * Z_eq = 2pF * 25Ω = 50ps.
// The SPICE netlist uses tmax=2ps so the adaptive stepper resolves the 50ps RC
// transient accurately (SPICE error <1mV).  The dominant error source is the
// backward-Euler ODE in the MoC sim: dt_SI/(2·τ_RC) ≈ 1ps/100ps ≈ 1%.
// With oversample=1024 (dt≈1ps) this bounds the MoC error to ~10mV.
// ────────────────────────────────────────────────────────────────────────────
{
  const tsvPath = path.join(__dirname, "sim", "results-tline-shunt-c.tsv");
  const skip    = !existsSync(tsvPath);

  test("SPICE comparison: shunt-C at midpoint (tline-shunt-c.sp)",
    { skip: skip ? "TSV not found — run: cd sim && ngspice tline-shunt-c.sp" : false },
    () => {
      const TAU_D   = 1e-9;
      const TR_NORM = 0.1;
      const TOL     = 2e-2;
      const model = {
        Vg: 1, Rg: 10,
        tau_d: TAU_D,
        riseTimeTr: TR_NORM, riseShape: "linear",
        blocks: [
          { type: "tl", Z0: 50, tau: 1 },
          { type: "C",  value: 2e-12 },
          { type: "tl", Z0: 50, tau: 1 },
        ],
        terminal: { type: "open" },
      };
      const sim = simulateTimeDomain(model, { tEnd: 14, oversample: 1024 });
      const { dt, nSteps, nodeV } = sim;

      // Linear interpolation between adjacent sim steps to reduce quantisation error.
      function interp(arr, tn) {
        const f  = tn / dt;
        const k0 = Math.floor(f);
        const k1 = Math.min(k0 + 1, nSteps - 1);
        return arr[k0] * (1 - (f - k0)) + arr[k1] * (f - k0);
      }

      const rows = readFileSync(tsvPath, "utf8").trim().split("\n")
        .map((line) => line.trim().split(/\s+/).map(Number))
        .filter((cols) => cols.length >= 6);
      assert.ok(rows.length > 10, "TSV should have many rows");

      // SPICE produced ~7000 rows at 2ps spacing; sample every 25th to keep the
      // test fast while still covering the full 14ns window.
      for (let ri = 0; ri < rows.length; ri += 25) {
        const [t_s, va_spice, , vb_spice, , vc_spice] = rows[ri];
        const tn = t_s / TAU_D;
        near(interp(nodeV[0], tn), va_spice, `v(a) at t=${t_s.toExponential(3)}s`, TOL);
        near(interp(nodeV[1], tn), vb_spice, `v(b) at t=${t_s.toExponential(3)}s`, TOL);
        near(interp(nodeV[2], tn), vc_spice, `v(c) at t=${t_s.toExponential(3)}s`, TOL);
      }
    }
  );
}

// ────────────────────────────────────────────────────────────────────────────
// SPICE comparison: shunt-L at midpoint between two T-line segments
//
// Circuit: Rg=10Ω → T1(Z0=50,τ=1ns) → L1=10nH → T2(Z0=50,τ=1ns) → open
// Netlist: sim/tline-shunt-l.sp   TSV columns: t v(a) t v(b) t v(c)
//
// L=10nH, Z_eq=25Ω → τ_L = L/Z_eq = 400ps.  V_B rises then decays toward 0
// as the inductor builds up current and shunts the node.
// ────────────────────────────────────────────────────────────────────────────
{
  const tsvPath = path.join(__dirname, "sim", "results-tline-shunt-l.tsv");
  const skip    = !existsSync(tsvPath);

  test("SPICE comparison: shunt-L at midpoint (tline-shunt-l.sp)",
    { skip: skip ? "TSV not found — run: cd sim && ngspice tline-shunt-l.sp" : false },
    () => {
      const TAU_D   = 1e-9;
      const TR_NORM = 0.1;
      const TOL     = 2e-2;
      const model = {
        Vg: 1, Rg: 10,
        tau_d: TAU_D,
        riseTimeTr: TR_NORM, riseShape: "linear",
        blocks: [
          { type: "tl", Z0: 50, tau: 1 },
          { type: "L",  value: 10e-9 },
          { type: "tl", Z0: 50, tau: 1 },
        ],
        terminal: { type: "open" },
      };
      const sim = simulateTimeDomain(model, { tEnd: 14, oversample: 1024 });
      const { dt, nSteps, nodeV } = sim;

      function interp(arr, tn) {
        const f  = tn / dt;
        const k0 = Math.floor(f);
        const k1 = Math.min(k0 + 1, nSteps - 1);
        return arr[k0] * (1 - (f - k0)) + arr[k1] * (f - k0);
      }

      const rows = readFileSync(tsvPath, "utf8").trim().split("\n")
        .map((line) => line.trim().split(/\s+/).map(Number))
        .filter((cols) => cols.length >= 6);
      assert.ok(rows.length > 10, "TSV should have many rows");

      for (let ri = 0; ri < rows.length; ri += 25) {
        const [t_s, va_spice, , vb_spice, , vc_spice] = rows[ri];
        const tn = t_s / TAU_D;
        near(interp(nodeV[0], tn), va_spice, `v(a) at t=${t_s.toExponential(3)}s`, TOL);
        near(interp(nodeV[1], tn), vb_spice, `v(b) at t=${t_s.toExponential(3)}s`, TOL);
        near(interp(nodeV[2], tn), vc_spice, `v(c) at t=${t_s.toExponential(3)}s`, TOL);
      }
    }
  );
}

// ────────────────────────────────────────────────────────────────────────────
{
  const tsvPath = path.join(__dirname, "sim", "results-tline-shunt-r.tsv");
  const skip    = !existsSync(tsvPath);

  test("SPICE comparison: shunt-R at midpoint (tline-shunt-r.sp)",
    { skip: skip ? "TSV not found — run: cd sim && ngspice tline-shunt-r.sp" : false },
    () => {
      const TAU_D   = 1e-9;
      const TR_NORM = 0.1;
      const TOL     = 1e-2;
      const model = {
        Vg: 1, Rg: 10,
        riseTimeTr: TR_NORM, riseShape: "linear",
        blocks: [
          { type: "tl", Z0: 50, tau: 1 },
          { type: "R",  value: 100 },
          { type: "tl", Z0: 50, tau: 1 },
        ],
        terminal: { type: "open" },
      };
      const sim = simulateTimeDomain(model, { tEnd: 14, oversample: 256 });
      const { dt, nSteps, nodeV } = sim;

      // Linear interpolation between adjacent sim steps to reduce quantisation error.
      function interp(arr, tn) {
        const f  = tn / dt;
        const k0 = Math.floor(f);
        const k1 = Math.min(k0 + 1, nSteps - 1);
        return arr[k0] * (1 - (f - k0)) + arr[k1] * (f - k0);
      }

      const rows = readFileSync(tsvPath, "utf8").trim().split("\n")
        .map((line) => line.trim().split(/\s+/).map(Number))
        .filter((cols) => cols.length >= 6);
      assert.ok(rows.length > 10, "TSV should have many rows");

      for (const [t_s, va_spice, , vb_spice, , vc_spice] of rows) {
        const tn = t_s / TAU_D;
        near(interp(nodeV[0], tn), va_spice, `v(a) at t=${t_s.toExponential(3)}s`, TOL);
        near(interp(nodeV[1], tn), vb_spice, `v(b) at t=${t_s.toExponential(3)}s`, TOL);
        near(interp(nodeV[2], tn), vc_spice, `v(c) at t=${t_s.toExponential(3)}s`, TOL);
      }
    }
  );
}