// 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);