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 return {
315 BOUNCE_EPS,
316 MAX_BOUNCES,
317 computeWaveParams,
318 buildBounceSeries,
319 sumEventsAtTime,
320 sumEventsWithRise,
321 sumEventsWithLinearRamp,
322 riseShape,
323 riseShapeLinear,
324 waveVoltageAt,
325 totalVoltageAt,
326 segmentForZ,
327 computeDynamicState,
328 };
329})();