oscillatory.net
this repo has no description
1"use strict";
2
3// Transmission line physics — pure functions, no DOM or canvas dependencies.
4//
5// Spatial convention: z ∈ [0, 1] maps to [0, ℓ].
6// Time convention: tNorm is dimensionless time in units of τ_d (one-way delay).
7//
8// Single-segment model (N=1):
9// model.segments = [{Z0: 50}] — one segment, the whole line.
10// Degenerates exactly to the original two-terminal bounce series.
11//
12// Multi-segment model (N>1):
13// model.segments = [{Z0: Z1}, {Z0: Z2}, …] — N equal-length segments.
14// Each internal boundary produces a reflected and a transmitted wave (lattice diagram).
15// buildBounceSeries uses a BFS priority queue ordered by tDie.
16//
17// Wave packet struct (new):
18// { n, A, dir, zStart, zEnd, tBorn, tDie, segIdx }
19// n — sequential index (for display)
20// A — voltage amplitude
21// dir — +1 rightward, -1 leftward
22// zStart — z where this packet's front begins its journey
23// zEnd — z where this packet's front ends (next boundary, or 0/1)
24// tBorn — tNorm when the front was at zStart
25// tDie — tNorm when the front reaches zEnd (= tBorn + |zEnd − zStart|)
26// segIdx — index into model.segments for the segment this packet lives in
27//
28// After computeDynamicState, each wave also carries:
29// u — tNorm − tBorn (age of the packet)
30// front — current position of the leading edge ∈ [zStart, zEnd] (or reversed for leftward)
31//
32// Extension points:
33// • Rise time (causal shifted erf or linear ramp): model.riseTimeTr > 0 — see riseShape/waveVoltageAt.
34// • Non-equal segment lengths: store lengths[] summing to 1; the BFS solver already works
35// by position, so changing zBnds is all that's needed.
36
37const TLPhysics = (() => {
38 const BOUNCE_EPS = 1e-6;
39 const MAX_BOUNCES = 5000;
40
41 // ---- reflection coefficients and initial voltage step ----
42 // Uses the first segment's Z0 for the source divider and source-end Γ,
43 // and the last segment's Z0 for the load-end Γ.
44 function computeWaveParams(model) {
45 const { Vg, Rg, RL } = model;
46 const Z0_src = model.segments[0].Z0;
47 const Z0_load = model.segments[model.segments.length - 1].Z0;
48 const V1 = Vg * Z0_src / (Rg + Z0_src);
49 const gL = isFinite(RL) ? (RL - Z0_load) / (RL + Z0_load) : 1; // RL=∞ → open circuit
50 const gS = (Rg - Z0_src) / (Rg + Z0_src);
51 return { V1, gL, gS };
52 }
53
54 // ---- BFS lattice-diagram solver ----
55 // Builds the full set of wave packets produced by successive reflections and
56 // transmissions at every impedance boundary (source, internal, load).
57 //
58 // Returns {
59 // series: [{n, A, dir, zStart, zEnd, tBorn, tDie, segIdx}]
60 // srcEvents: [{t, dV}] — cumulative voltage steps at z=0 (VS)
61 // loadEvents: [{t, dV}] — cumulative voltage steps at z=ℓ (VL)
62 // gSEffective, tolPct, ampTol, tEnd
63 // }
64 function buildBounceSeries(model, waves) {
65 const segs = model.segments;
66 const N = segs.length;
67 const segLen = 1 / N;
68 // z positions of all N+1 boundaries: [0, 1/N, 2/N, ..., 1]
69 const zBnds = Array.from({ length: N + 1 }, (_, i) => i * segLen);
70
71 const gS = waves.gS;
72 const gL = waves.gL;
73 const tolPct = Number.isFinite(model.reflectTol) ? model.reflectTol : 1;
74 const ampTol = Math.max(BOUNCE_EPS, Math.abs(waves.V1) * (tolPct / 100));
75
76 const series = [];
77 const srcEvents = [{ t: 0, dV: waves.V1 }]; // initial step seen at source at t=0
78 const loadEvents = [];
79
80 // Queue of pending packets; each entry is a wave packet description.
81 // We process in order of tDie (earliest-completing first).
82 const queue = [];
83
84 function enqueue(pkt) {
85 if (Math.abs(pkt.A) <= ampTol) return;
86 if (series.length + queue.length >= MAX_BOUNCES) return;
87 queue.push(pkt);
88 }
89
90 // Seed: first incident wave in segment 0, rightward from z=0 to z=zBnds[1].
91 enqueue({ A: waves.V1, dir: +1, zStart: 0, zEnd: zBnds[1], tBorn: 0, segIdx: 0 });
92
93 while (queue.length > 0) {
94 // Pop earliest-completing packet (smallest tDie = tBorn + |zEnd − zStart|).
95 queue.sort((a, b) =>
96 (a.tBorn + Math.abs(a.zEnd - a.zStart)) -
97 (b.tBorn + Math.abs(b.zEnd - b.zStart))
98 );
99 const pkt = queue.shift();
100 const { A, dir, zStart, zEnd, tBorn, segIdx } = pkt;
101 const tDie = tBorn + Math.abs(zEnd - zStart);
102 const n = series.length + 1;
103 series.push({ n, A, dir, zStart, zEnd, tBorn, tDie, segIdx });
104
105 const atSource = zEnd < 1e-9;
106 const atLoad = zEnd > 1 - 1e-9;
107
108 if (atLoad) {
109 // Load terminal: add load event, spawn reflected wave back in last segment.
110 loadEvents.push({ t: tDie, dV: (1 + gL) * A });
111 enqueue({
112 A: gL * A, dir: -1,
113 zStart: 1, zEnd: zBnds[N - 1],
114 tBorn: tDie, segIdx: N - 1,
115 });
116 } else if (atSource) {
117 // Source terminal: add source event, spawn reflected wave into first segment.
118 srcEvents.push({ t: tDie, dV: (1 + gS) * A });
119 enqueue({
120 A: gS * A, dir: +1,
121 zStart: 0, zEnd: zBnds[1],
122 tBorn: tDie, segIdx: 0,
123 });
124 } else {
125 // Internal boundary between two segments.
126 // For a rightward wave in segIdx: left=segIdx, right=segIdx+1.
127 // For a leftward wave in segIdx: left=segIdx-1, right=segIdx.
128 const leftSeg = dir > 0 ? segIdx : segIdx - 1;
129 const rightSeg = dir > 0 ? segIdx + 1 : segIdx;
130 const Zl = segs[leftSeg].Z0;
131 const Zr = segs[rightSeg].Z0;
132 // Γ for a wave arriving from the left side of this boundary:
133 const gBound = (Zr - Zl) / (Zr + Zl);
134
135 if (dir > 0) {
136 // Rightward hitting boundary from left:
137 // reflected → Γ·A leftward, back into leftSeg
138 // transmitted → (1+Γ)·A rightward, into rightSeg
139 enqueue({
140 A: gBound * A, dir: -1,
141 zStart: zEnd, zEnd: zBnds[leftSeg],
142 tBorn: tDie, segIdx: leftSeg,
143 });
144 enqueue({
145 A: (1 + gBound) * A, dir: +1,
146 zStart: zEnd, zEnd: zBnds[rightSeg + 1],
147 tBorn: tDie, segIdx: rightSeg,
148 });
149 } else {
150 // Leftward hitting boundary from right:
151 // reflected → −Γ·A rightward, back into rightSeg
152 // transmitted → (1−Γ)·A leftward, into leftSeg
153 enqueue({
154 A: -gBound * A, dir: +1,
155 zStart: zEnd, zEnd: zBnds[rightSeg + 1],
156 tBorn: tDie, segIdx: rightSeg,
157 });
158 enqueue({
159 A: (1 - gBound) * A, dir: -1,
160 zStart: zEnd, zEnd: zBnds[leftSeg],
161 tBorn: tDie, segIdx: leftSeg,
162 });
163 }
164 }
165 }
166
167 const tEnd = series.reduce((m, w) => Math.max(m, w.tDie), 2);
168 return { series, srcEvents, loadEvents, gSEffective: gS, tolPct, ampTol, tEnd };
169 }
170
171 // ---- time query: node voltage (step model) ----
172 function sumEventsAtTime(events, tn) {
173 const eps = 1e-9;
174 let v = 0;
175 for (const e of events) {
176 if (tn >= e.t - eps) v += e.dV;
177 }
178 return v;
179 }
180
181 // ---- rise-time wave shapes ----
182 // Fraction of final amplitude reached after time dt since the wave front passed.
183
184 // Abramowitz & Stegun 7.1.26 rational approximation (max error < 1.5e-7).
185 function erf(x) {
186 const sign = x < 0 ? -1 : 1;
187 x = Math.abs(x);
188 const t = 1 / (1 + 0.3275911 * x);
189 const y = 1 - (((((1.061405429 * t - 1.453152027) * t) + 1.421413741) * t - 0.284496736) * t + 0.254829592) * t * Math.exp(-x * x);
190 return sign * y;
191 }
192
193 // Causal erf S-curve: tr = 0 → hard step. tr > 0 → smooth S from 0 to 1 over [0, tr].
194 // Uses shifted erf: 0.5·(1 + erf(k·(2·dt/tr − 1))) with k chosen so that
195 // riseShape(0, tr) ≈ 0 and riseShape(tr, tr) ≈ 1. k = 1.8 gives < 0.1% undershoot.
196 // Causal: strictly 0 for dt ≤ 0, strictly 1 for dt ≥ tr.
197 const ERF_K = 1.8;
198 function riseShape(dt, tr) {
199 if (dt <= 0) return 0;
200 if (tr <= 0 || dt >= tr) return 1;
201 return 0.5 * (1 + erf(ERF_K * (2 * dt / tr - 1)));
202 }
203
204 // Linear ramp (trapezoidal / SPICE PULSE): tr = 0 → hard step. tr > 0 → dt/tr clamped to [0,1].
205 function riseShapeLinear(dt, tr) {
206 if (dt <= 0) return 0;
207 if (tr <= 0 || dt >= tr) return 1;
208 return dt / tr;
209 }
210
211 // Voltage contribution of a single wave packet wf at position z.
212 // wf must carry the fields added by computeDynamicState: u, front.
213 // Returns 0 if z is outside this packet's z-range or ahead of the front.
214 // shape: "linear" → riseShapeLinear, "erf" → riseShape (causal shifted erf).
215 // tr: rise time in τ_d units (used by both shapes).
216 function waveVoltageAt(wf, z, tr = 0, shape = "erf") {
217 const { dir, zStart, A, u, front } = wf;
218 const eps = 1e-9;
219 const rise = shape === "linear" ? riseShapeLinear : riseShape;
220 if (dir > 0) {
221 // Rightward: nonzero in [zStart, front].
222 if (z < zStart - eps || z > front + eps) return 0;
223 const dt = u - (z - zStart);
224 return A * rise(dt, tr);
225 } else {
226 // Leftward: nonzero in [front, zStart].
227 if (z > zStart + eps || z < front - eps) return 0;
228 const dt = u - (zStart - z);
229 return A * rise(dt, tr);
230 }
231 }
232
233 // Map z ∈ [0,1] to the index of the segment that owns it.
234 // Boundaries are assigned to the right-hand (higher-index) segment; z=1 goes to the last.
235 // The +1e-9 nudge handles floating-point imprecision when z is exactly i/N.
236 function segmentForZ(z, N) {
237 if (z >= 1 - 1e-9) return N - 1;
238 return Math.min(N - 1, Math.floor(z * N + 1e-9));
239 }
240
241 // Sum of all wave contributions at position z.
242 // For multi-segment models (N > 1), only waves whose segIdx matches the segment
243 // that owns z are summed. This prevents double-counting at segment boundaries
244 // where adjacent waves share a z value (a plotting artifact, not a physics error).
245 function totalVoltageAt(z, launchedWaves, tr = 0, shape = "erf", N = 1) {
246 const segIdx = N > 1 ? segmentForZ(z, N) : 0;
247 let V = 0;
248 for (const wf of launchedWaves) {
249 if (N > 1 && wf.segIdx !== segIdx) continue;
250 V += waveVoltageAt(wf, z, tr, shape);
251 }
252 return V;
253 }
254
255 // Node voltage as smooth sum of causal erf-rise events.
256 function sumEventsWithRise(events, tn, tr) {
257 let v = 0;
258 for (const e of events) {
259 const dt = tn - e.t;
260 if (dt > 0) v += e.dV * riseShape(dt, tr);
261 }
262 return v;
263 }
264
265 // Node voltage as smooth sum of linear-ramp events (matches SPICE PULSE with finite TR).
266 function sumEventsWithLinearRamp(events, tn, tr) {
267 let v = 0;
268 for (const e of events) {
269 const dt = tn - e.t;
270 if (dt > 0) v += e.dV * riseShapeLinear(dt, tr);
271 }
272 return v;
273 }
274
275 // ---- dynamic state at a given tNorm ----
276 // Annotates each launched wave with:
277 // u: tNorm − tBorn
278 // front: current leading-edge position
279 // A wave is "active" while u < |zEnd − zStart| (front still in transit).
280 // shape: "step" | "linear" | "erf"
281 function computeDynamicState(tn, bounce, riseTimeTr = 0, riseShape_ = "erf") {
282 const { clamp } = TLUtils;
283 const launchedWaves = [];
284 const activeWaves = [];
285
286 for (const w of bounce.series) {
287 const u = tn - w.tBorn;
288 if (u < 0) continue;
289 let front;
290 if (w.dir > 0) {
291 front = clamp(w.zStart + u, w.zStart, w.zEnd);
292 } else {
293 front = clamp(w.zStart - u, w.zEnd, w.zStart);
294 }
295 const ww = { ...w, u, front };
296 launchedWaves.push(ww);
297 if (u < Math.abs(w.zEnd - w.zStart) - 1e-9) activeWaves.push(ww);
298 }
299
300 const sumNode = (events) => {
301 if (riseTimeTr <= 0) return sumEventsAtTime(events, tn);
302 if (riseShape_ === "linear") return sumEventsWithLinearRamp(events, tn, riseTimeTr);
303 return sumEventsWithRise(events, tn, riseTimeTr);
304 };
305
306 return {
307 launchedWaves,
308 activeWaves,
309 VS: sumNode(bounce.srcEvents),
310 VL: sumNode(bounce.loadEvents),
311 };
312 }
313
314 // ---- Method of Characteristics time-domain simulator ----
315 //
316 // model.blocks = [
317 // { type: 'tl', Z0, tau }, // T-line segment; tau = one-way delay in τ_d units
318 // { type: 'R', value }, // shunt resistor to ground (Ω)
319 // { type: 'C', value }, // shunt capacitor (F, SI) — Step 3
320 // { type: 'L', value }, // shunt inductor (H, SI) — Step 4
321 // { type: 'short' }, // shunt short to ground
322 // ]
323 // model.terminal = { type: 'R'|'open'|'short'|'C'|'L', value }
324 // model.Rg, model.Vg, model.riseTimeTr, model.riseShape: unchanged
325 //
326 // Consecutive shunt blocks (no 'tl' between them) share one electrical node;
327 // their admittances accumulate in the node equation.
328 // Consecutive 'tl' blocks share a pure impedance-mismatch boundary (no shunt).
329 //
330 // opts: { tEnd, oversample }
331 // tEnd — simulation end time in τ_d units (default 10)
332 // oversample — delay steps per shortest segment (default 8)
333 //
334 // Returns { dt, nSteps, tGrid, segs, nDelay, nodeV, vPlusHist, vMinusHist }
335 // dt — time step in τ_d units
336 // tGrid — Float64Array of length nSteps: tGrid[k] = k * dt
337 // segs — [{Z0, tau}] after parsing (T-line segments only)
338 // nDelay — Int32Array: nDelay[i] = delay of segment i in time steps
339 // nodeV — (N+1) × nSteps: nodeV[j][k] = voltage at node j at step k
340 // vPlusHist — N × nSteps: V⁺ written into left end of segment i at step k
341 // vMinusHist — N × nSteps: V⁻ written into right end of segment i at step k
342 //
343 // Spatial query — see voltageAt().
344 function simulateTimeDomain(model, opts) {
345 const { Vg: Vg_final, Rg } = model;
346 const tEnd_ = (opts && opts.tEnd) != null ? opts.tEnd : 10;
347 const oversample = (opts && opts.oversample) != null ? opts.oversample : 8;
348
349 // --- parse model.blocks into T-line segments + per-node shunt lists ---
350 // nodeShunts[j] are the shunt elements sitting at node j.
351 // Node j sits between segment j-1 (on the left) and segment j (on the right).
352 // Node 0 = source side; node N = load side.
353 const segs = []; // {Z0, tau}
354 const nodeShunts = []; // nodeShunts[j] = [{type, value?}]
355 let pending = []; // shunts accumulating until the next 'tl' block
356 for (const blk of model.blocks) {
357 if (blk.type === 'tl') {
358 nodeShunts.push(pending);
359 pending = [];
360 segs.push({ Z0: blk.Z0, tau: blk.tau });
361 } else {
362 pending.push(blk);
363 }
364 }
365 nodeShunts.push(pending); // shunts at load-side node (right of last T-line)
366 const N = segs.length;
367 if (N === 0) throw new Error('simulateTimeDomain: no tl blocks in model.blocks');
368
369 // --- time step and per-segment delay lengths ---
370 // dt chosen so the shortest segment spans exactly `oversample` steps.
371 const minTau = segs.reduce((m, s) => Math.min(m, s.tau), Infinity);
372 const dt = minTau / oversample;
373 const nSteps = Math.ceil(tEnd_ / dt) + 1;
374 // dt_SI: real-time duration of one time step (seconds).
375 // model.tau_d must be provided (SI seconds per tau unit) when any shunt C or L
376 // is present; for purely resistive circuits it is unused (default 1 is harmless).
377 const tau_d_SI = model.tau_d != null ? model.tau_d : 1;
378 const dt_SI = dt * tau_d_SI;
379 // Each segment rounded to the nearest integer number of steps (≥ 1).
380 const nDelay = new Int32Array(N);
381 for (let i = 0; i < N; i++) nDelay[i] = Math.max(1, Math.round(segs[i].tau / dt));
382
383 // --- ring buffers ---
384 // delayFwd[i]: V⁺ entering left end of segment i.
385 // Written at step k, read (= arrives at right end) at step k + nDelay[i].
386 // delayBwd[i]: V⁻ entering right end of segment i.
387 // Written at step k, arrives at left end at step k + nDelay[i].
388 // Both use the same ring-buffer idiom: fwdIdx[i] / bwdIdx[i] is the current
389 // write position; that same position holds the oldest (= arriving) value.
390 const delayFwd = Array.from({ length: N }, (_, i) => new Float64Array(nDelay[i]));
391 const delayBwd = Array.from({ length: N }, (_, i) => new Float64Array(nDelay[i]));
392 const fwdIdx = new Int32Array(N);
393 const bwdIdx = new Int32Array(N);
394
395 // --- reactive state per node (C: vCap; L: iL) ---
396 // Each element of nodeShuntState[j] is a mutable object parallel to nodeShunts[j].
397 const nodeShuntState = nodeShunts.map(shunts =>
398 shunts.map(sh => ({ vCap: 0, iL: 0 }))
399 );
400 const termState = { vCap: 0, iL: 0 };
401
402 // --- output storage ---
403 const nodeV = Array.from({ length: N + 1 }, () => new Float64Array(nSteps));
404 const vPlusHist = Array.from({ length: N }, () => new Float64Array(nSteps));
405 const vMinusHist = Array.from({ length: N }, () => new Float64Array(nSteps));
406
407 // --- source waveform ---
408 const riseTimeTr_ = model.riseTimeTr || 0;
409 const riseMode_ = model.riseShape || 'step';
410 function vgAt(tn) {
411 if (riseTimeTr_ <= 0) return tn > -1e-12 ? Vg_final : 0;
412 if (riseMode_ === 'linear') return Vg_final * riseShapeLinear(tn, riseTimeTr_);
413 return Vg_final * riseShape(tn, riseTimeTr_);
414 }
415
416 // --- helper: accumulate shunt contributions into node equation ---
417 // Each shunt contributes conductance G and Norton current I_src such that:
418 // V * (G_base + ΣG) = I_base + ΣI_src
419 // where G_base / I_base come from the Thevenin equivalent(s) of adjacent T-lines.
420 // Returns true if the node is hard-shorted to ground.
421 function accumulateShunts(shunts, states, Vnode_prev_fn) {
422 let G = 0, I = 0, shorted = false;
423 for (let s = 0; s < shunts.length; s++) {
424 const sh = shunts[s], st = states[s];
425 if (sh.type === 'R') {
426 G += 1 / sh.value;
427 } else if (sh.type === 'C') {
428 const Gc = sh.value / dt_SI;
429 G += Gc;
430 I += Gc * st.vCap;
431 } else if (sh.type === 'L') {
432 const Gl = dt_SI / sh.value;
433 G += Gl;
434 I -= st.iL; // see derivation: inductor history current subtracts
435 } else if (sh.type === 'short') {
436 shorted = true;
437 }
438 }
439 return { G, I, shorted };
440 }
441
442 function updateShuntStates(shunts, states, V) {
443 for (let s = 0; s < shunts.length; s++) {
444 const sh = shunts[s], st = states[s];
445 if (sh.type === 'C') st.vCap = V;
446 if (sh.type === 'L') st.iL += (dt_SI / sh.value) * V;
447 }
448 }
449
450 // --- main simulation loop ---
451 const vArr = {
452 right: new Float64Array(N), // V⁺ arriving at right end of each segment
453 left: new Float64Array(N), // V⁻ arriving at left end of each segment
454 };
455 const Vnodes = new Float64Array(N + 1);
456
457 for (let step = 0; step < nSteps; step++) {
458 const tn = step * dt;
459
460 // 1. Read arriving waves from delay lines.
461 for (let i = 0; i < N; i++) {
462 vArr.right[i] = delayFwd[i][fwdIdx[i]];
463 vArr.left[i] = delayBwd[i][bwdIdx[i]];
464 }
465
466 // 2. Solve node voltages.
467
468 // Source node (j=0): Thevenin from {Vg, Rg} on left; segment 0 on right.
469 // KCL: (vg - V)/Rg + (2·aR - V)/Z0 = 0 → V = (vg·Z0 + 2·aR·Rg)/(Rg + Z0)
470 // (Source-node shunts, if any, are rare but handled.)
471 {
472 const vg = vgAt(tn);
473 const aR = vArr.left[0]; // V⁻ arriving at left end of seg 0
474 const Z0 = segs[0].Z0;
475 let G = 1 / Rg + 1 / Z0;
476 let I = vg / Rg + 2 * aR / Z0;
477 const { G: Gs, I: Is, shorted } = accumulateShunts(
478 nodeShunts[0], nodeShuntState[0]);
479 if (shorted) {
480 Vnodes[0] = 0;
481 } else {
482 Vnodes[0] = (I + Is) / (G + Gs);
483 }
484 updateShuntStates(nodeShunts[0], nodeShuntState[0], Vnodes[0]);
485 }
486
487 // Internal nodes (j = 1 … N-1).
488 for (let j = 1; j < N; j++) {
489 const ZL = segs[j - 1].Z0;
490 const ZR = segs[j].Z0;
491 const aL = vArr.right[j - 1]; // V⁺ arriving at right end of seg j-1
492 const aR = vArr.left[j]; // V⁻ arriving at left end of seg j
493 let G = 1 / ZL + 1 / ZR;
494 let I = 2 * aL / ZL + 2 * aR / ZR;
495 const { G: Gs, I: Is, shorted } = accumulateShunts(
496 nodeShunts[j], nodeShuntState[j]);
497 if (shorted) {
498 Vnodes[j] = 0;
499 } else {
500 Vnodes[j] = (I + Is) / (G + Gs);
501 }
502 updateShuntStates(nodeShunts[j], nodeShuntState[j], Vnodes[j]);
503 }
504
505 // Load node (j=N): segment N-1 on left; terminal on right.
506 {
507 const aL = vArr.right[N - 1]; // V⁺ arriving at right end of last seg
508 const Z0 = segs[N - 1].Z0;
509 const term = model.terminal;
510 // First apply any mid-node shunts at the load-side node (nodeShunts[N]).
511 let G = 1 / Z0;
512 let I = 2 * aL / Z0;
513 const { G: Gs, I: Is, shorted: sh0 } = accumulateShunts(
514 nodeShunts[N], nodeShuntState[N]);
515 if (sh0) {
516 Vnodes[N] = 0;
517 } else {
518 // Then add the terminal element.
519 let Gt = 0, It = 0, shorted = false;
520 if (term.type === 'open') {
521 // no terminal conductance
522 } else if (term.type === 'short') {
523 shorted = true;
524 } else if (term.type === 'R') {
525 Gt = 1 / term.value;
526 } else if (term.type === 'C') {
527 Gt = term.value / dt_SI;
528 It = Gt * termState.vCap;
529 } else if (term.type === 'L') {
530 Gt = dt_SI / term.value;
531 It = -termState.iL;
532 }
533 Vnodes[N] = shorted ? 0 : (I + Is + It) / (G + Gs + Gt);
534 }
535 updateShuntStates(nodeShunts[N], nodeShuntState[N], Vnodes[N]);
536 if (term.type === 'C') termState.vCap = Vnodes[N];
537 if (term.type === 'L') termState.iL += (dt_SI / term.value) * Vnodes[N];
538 }
539
540 // 3. Compute outgoing waves and write into delay lines.
541 // At node j: V = V⁺_out + V⁻_in → V⁺_out = V - V⁻_in (outgoing into seg j rightward)
542 // At node j: V = V⁺_in + V⁻_out → V⁻_out = V - V⁺_in (outgoing into seg j-1 leftward)
543 for (let i = 0; i < N; i++) {
544 const vOutFwd = Vnodes[i] - vArr.left[i]; // V⁺ leaving node i into seg i
545 const vOutBwd = Vnodes[i + 1] - vArr.right[i]; // V⁻ leaving node i+1 into seg i
546 vPlusHist[i][step] = vOutFwd;
547 vMinusHist[i][step] = vOutBwd;
548 delayFwd[i][fwdIdx[i]] = vOutFwd;
549 delayBwd[i][bwdIdx[i]] = vOutBwd;
550 fwdIdx[i] = (fwdIdx[i] + 1) % nDelay[i];
551 bwdIdx[i] = (bwdIdx[i] + 1) % nDelay[i];
552 }
553
554 // 4. Record node voltages.
555 for (let j = 0; j <= N; j++) nodeV[j][step] = Vnodes[j];
556 }
557
558 const tGrid = Float64Array.from({ length: nSteps }, (_, k) => k * dt);
559 return { dt, nSteps, tGrid, segs, nDelay, nodeV, vPlusHist, vMinusHist };
560 }
561
562 // ---- Spatial voltage query for simulateTimeDomain output ----
563 //
564 // Returns V(z, tIdx) by reading from vPlusHist / vMinusHist at the appropriate
565 // delay offsets. z ∈ [0, 1] is normalized by total T-line delay (sum of taus).
566 //
567 // Within segment i, at fractional position f ∈ [0, 1]:
568 // V⁺ that is at position f at step tIdx was written into the segment
569 // kFwd = round(f * D) steps ago → look up vPlusHist[i][tIdx - kFwd]
570 // V⁻ that is at position f at step tIdx was written into the segment
571 // kBwd = round((1-f) * D) steps ago → look up vMinusHist[i][tIdx - kBwd]
572 function voltageAt(z, tIdx, sim) {
573 const { segs, nDelay, vPlusHist, vMinusHist, nodeV } = sim;
574 const N = segs.length;
575 const totalTau = segs.reduce((s, seg) => s + seg.tau, 0);
576 let cumZ = 0;
577 for (let i = 0; i < N; i++) {
578 const segZ = segs[i].tau / totalTau;
579 if (z <= cumZ + segZ + 1e-9) {
580 const f = Math.min(1, Math.max(0, (z - cumZ) / segZ));
581 const D = nDelay[i];
582 const kFwd = Math.round(f * D);
583 const kBwd = D - kFwd;
584 const tF = tIdx - kFwd;
585 const tB = tIdx - kBwd;
586 const vP = tF >= 0 ? vPlusHist[i][tF] : 0;
587 const vM = tB >= 0 ? vMinusHist[i][tB] : 0;
588 return vP + vM;
589 }
590 cumZ += segZ;
591 }
592 // z at or past the load end: return load node voltage
593 return nodeV[N][tIdx];
594 }
595
596 // ---- Convenience: convert old-style model (model.segments, model.RL) to blocks ----
597 // This lets existing callers construct a blocks-based model without rewriting app.js.
598 function segmentsToBlocks(model) {
599 const N = model.segments.length;
600 const tau = 1 / N; // each segment = 1/N of total normalized delay
601 const blocks = model.segments.map(seg => ({ type: 'tl', Z0: seg.Z0, tau }));
602 const RL = model.RL;
603 const terminal = isFinite(RL)
604 ? { type: 'R', value: RL }
605 : { type: 'open' };
606 return { ...model, blocks, terminal };
607 }
608
609 return {
610 BOUNCE_EPS,
611 MAX_BOUNCES,
612 computeWaveParams,
613 buildBounceSeries,
614 sumEventsAtTime,
615 sumEventsWithRise,
616 sumEventsWithLinearRamp,
617 riseShape,
618 riseShapeLinear,
619 waveVoltageAt,
620 totalVoltageAt,
621 segmentForZ,
622 computeDynamicState,
623 simulateTimeDomain,
624 voltageAt,
625 segmentsToBlocks,
626 };
627})();