"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),
    };
  }

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