+293

physics.js

reviewed

··· 311 311 }; 312 312 } 313 313 314 314 + // ---- Method of Characteristics time-domain simulator ---- 315 315 + // 316 316 + // model.blocks = [ 317 317 + // { type: 'tl', Z0, tau }, // T-line segment; tau = one-way delay in τ_d units 318 318 + // { type: 'R', value }, // shunt resistor to ground (Ω) 319 319 + // { type: 'C', value }, // shunt capacitor (F, SI) — Step 3 320 320 + // { type: 'L', value }, // shunt inductor (H, SI) — Step 4 321 321 + // { type: 'short' }, // shunt short to ground 322 322 + // ] 323 323 + // model.terminal = { type: 'R'|'open'|'short'|'C'|'L', value } 324 324 + // model.Rg, model.Vg, model.riseTimeTr, model.riseShape: unchanged 325 325 + // 326 326 + // Consecutive shunt blocks (no 'tl' between them) share one electrical node; 327 327 + // their admittances accumulate in the node equation. 328 328 + // Consecutive 'tl' blocks share a pure impedance-mismatch boundary (no shunt). 329 329 + // 330 330 + // opts: { tEnd, oversample } 331 331 + // tEnd — simulation end time in τ_d units (default 10) 332 332 + // oversample — delay steps per shortest segment (default 8) 333 333 + // 334 334 + // Returns { dt, nSteps, tGrid, segs, nDelay, nodeV, vPlusHist, vMinusHist } 335 335 + // dt — time step in τ_d units 336 336 + // tGrid — Float64Array of length nSteps: tGrid[k] = k * dt 337 337 + // segs — [{Z0, tau}] after parsing (T-line segments only) 338 338 + // nDelay — Int32Array: nDelay[i] = delay of segment i in time steps 339 339 + // nodeV — (N+1) × nSteps: nodeV[j][k] = voltage at node j at step k 340 340 + // vPlusHist — N × nSteps: V⁺ written into left end of segment i at step k 341 341 + // vMinusHist — N × nSteps: V⁻ written into right end of segment i at step k 342 342 + // 343 343 + // Spatial query — see voltageAt(). 344 344 + function simulateTimeDomain(model, opts) { 345 345 + const { Vg: Vg_final, Rg } = model; 346 346 + const tEnd_ = (opts && opts.tEnd) != null ? opts.tEnd : 10; 347 347 + const oversample = (opts && opts.oversample) != null ? opts.oversample : 8; 348 348 + 349 349 + // --- parse model.blocks into T-line segments + per-node shunt lists --- 350 350 + // nodeShunts[j] are the shunt elements sitting at node j. 351 351 + // Node j sits between segment j-1 (on the left) and segment j (on the right). 352 352 + // Node 0 = source side; node N = load side. 353 353 + const segs = []; // {Z0, tau} 354 354 + const nodeShunts = []; // nodeShunts[j] = [{type, value?}] 355 355 + let pending = []; // shunts accumulating until the next 'tl' block 356 356 + for (const blk of model.blocks) { 357 357 + if (blk.type === 'tl') { 358 358 + nodeShunts.push(pending); 359 359 + pending = []; 360 360 + segs.push({ Z0: blk.Z0, tau: blk.tau }); 361 361 + } else { 362 362 + pending.push(blk); 363 363 + } 364 364 + } 365 365 + nodeShunts.push(pending); // shunts at load-side node (right of last T-line) 366 366 + const N = segs.length; 367 367 + if (N === 0) throw new Error('simulateTimeDomain: no tl blocks in model.blocks'); 368 368 + 369 369 + // --- time step and per-segment delay lengths --- 370 370 + // dt chosen so the shortest segment spans exactly `oversample` steps. 371 371 + const minTau = segs.reduce((m, s) => Math.min(m, s.tau), Infinity); 372 372 + const dt = minTau / oversample; 373 373 + const nSteps = Math.ceil(tEnd_ / dt) + 1; 374 374 + // Each segment rounded to the nearest integer number of steps (≥ 1). 375 375 + const nDelay = new Int32Array(N); 376 376 + for (let i = 0; i < N; i++) nDelay[i] = Math.max(1, Math.round(segs[i].tau / dt)); 377 377 + 378 378 + // --- ring buffers --- 379 379 + // delayFwd[i]: V⁺ entering left end of segment i. 380 380 + // Written at step k, read (= arrives at right end) at step k + nDelay[i]. 381 381 + // delayBwd[i]: V⁻ entering right end of segment i. 382 382 + // Written at step k, arrives at left end at step k + nDelay[i]. 383 383 + // Both use the same ring-buffer idiom: fwdIdx[i] / bwdIdx[i] is the current 384 384 + // write position; that same position holds the oldest (= arriving) value. 385 385 + const delayFwd = Array.from({ length: N }, (_, i) => new Float64Array(nDelay[i])); 386 386 + const delayBwd = Array.from({ length: N }, (_, i) => new Float64Array(nDelay[i])); 387 387 + const fwdIdx = new Int32Array(N); 388 388 + const bwdIdx = new Int32Array(N); 389 389 + 390 390 + // --- reactive state per node (C: vCap; L: iL) --- 391 391 + // Each element of nodeShuntState[j] is a mutable object parallel to nodeShunts[j]. 392 392 + const nodeShuntState = nodeShunts.map(shunts => 393 393 + shunts.map(sh => ({ vCap: 0, iL: 0 })) 394 394 + ); 395 395 + const termState = { vCap: 0, iL: 0 }; 396 396 + 397 397 + // --- output storage --- 398 398 + const nodeV = Array.from({ length: N + 1 }, () => new Float64Array(nSteps)); 399 399 + const vPlusHist = Array.from({ length: N }, () => new Float64Array(nSteps)); 400 400 + const vMinusHist = Array.from({ length: N }, () => new Float64Array(nSteps)); 401 401 + 402 402 + // --- source waveform --- 403 403 + const riseTimeTr_ = model.riseTimeTr || 0; 404 404 + const riseMode_ = model.riseShape || 'step'; 405 405 + function vgAt(tn) { 406 406 + if (riseTimeTr_ <= 0) return tn > -1e-12 ? Vg_final : 0; 407 407 + if (riseMode_ === 'linear') return Vg_final * riseShapeLinear(tn, riseTimeTr_); 408 408 + return Vg_final * riseShape(tn, riseTimeTr_); 409 409 + } 410 410 + 411 411 + // --- helper: accumulate shunt contributions into node equation --- 412 412 + // Each shunt contributes conductance G and Norton current I_src such that: 413 413 + // V * (G_base + ΣG) = I_base + ΣI_src 414 414 + // where G_base / I_base come from the Thevenin equivalent(s) of adjacent T-lines. 415 415 + // Returns true if the node is hard-shorted to ground. 416 416 + function accumulateShunts(shunts, states, Vnode_prev_fn) { 417 417 + let G = 0, I = 0, shorted = false; 418 418 + for (let s = 0; s < shunts.length; s++) { 419 419 + const sh = shunts[s], st = states[s]; 420 420 + if (sh.type === 'R') { 421 421 + G += 1 / sh.value; 422 422 + } else if (sh.type === 'C') { 423 423 + const Gc = sh.value / dt; 424 424 + G += Gc; 425 425 + I += Gc * st.vCap; 426 426 + } else if (sh.type === 'L') { 427 427 + const Gl = dt / sh.value; 428 428 + G += Gl; 429 429 + I -= st.iL; // see derivation: inductor history current subtracts 430 430 + } else if (sh.type === 'short') { 431 431 + shorted = true; 432 432 + } 433 433 + } 434 434 + return { G, I, shorted }; 435 435 + } 436 436 + 437 437 + function updateShuntStates(shunts, states, V) { 438 438 + for (let s = 0; s < shunts.length; s++) { 439 439 + const sh = shunts[s], st = states[s]; 440 440 + if (sh.type === 'C') st.vCap = V; 441 441 + if (sh.type === 'L') st.iL += (dt / sh.value) * V; 442 442 + } 443 443 + } 444 444 + 445 445 + // --- main simulation loop --- 446 446 + const vArr = { 447 447 + right: new Float64Array(N), // V⁺ arriving at right end of each segment 448 448 + left: new Float64Array(N), // V⁻ arriving at left end of each segment 449 449 + }; 450 450 + const Vnodes = new Float64Array(N + 1); 451 451 + 452 452 + for (let step = 0; step < nSteps; step++) { 453 453 + const tn = step * dt; 454 454 + 455 455 + // 1. Read arriving waves from delay lines. 456 456 + for (let i = 0; i < N; i++) { 457 457 + vArr.right[i] = delayFwd[i][fwdIdx[i]]; 458 458 + vArr.left[i] = delayBwd[i][bwdIdx[i]]; 459 459 + } 460 460 + 461 461 + // 2. Solve node voltages. 462 462 + 463 463 + // Source node (j=0): Thevenin from {Vg, Rg} on left; segment 0 on right. 464 464 + // KCL: (vg - V)/Rg + (2·aR - V)/Z0 = 0 → V = (vg·Z0 + 2·aR·Rg)/(Rg + Z0) 465 465 + // (Source-node shunts, if any, are rare but handled.) 466 466 + { 467 467 + const vg = vgAt(tn); 468 468 + const aR = vArr.left[0]; // V⁻ arriving at left end of seg 0 469 469 + const Z0 = segs[0].Z0; 470 470 + let G = 1 / Rg + 1 / Z0; 471 471 + let I = vg / Rg + 2 * aR / Z0; 472 472 + const { G: Gs, I: Is, shorted } = accumulateShunts( 473 473 + nodeShunts[0], nodeShuntState[0]); 474 474 + if (shorted) { 475 475 + Vnodes[0] = 0; 476 476 + } else { 477 477 + Vnodes[0] = (I + Is) / (G + Gs); 478 478 + } 479 479 + updateShuntStates(nodeShunts[0], nodeShuntState[0], Vnodes[0]); 480 480 + } 481 481 + 482 482 + // Internal nodes (j = 1 … N-1). 483 483 + for (let j = 1; j < N; j++) { 484 484 + const ZL = segs[j - 1].Z0; 485 485 + const ZR = segs[j].Z0; 486 486 + const aL = vArr.right[j - 1]; // V⁺ arriving at right end of seg j-1 487 487 + const aR = vArr.left[j]; // V⁻ arriving at left end of seg j 488 488 + let G = 1 / ZL + 1 / ZR; 489 489 + let I = 2 * aL / ZL + 2 * aR / ZR; 490 490 + const { G: Gs, I: Is, shorted } = accumulateShunts( 491 491 + nodeShunts[j], nodeShuntState[j]); 492 492 + if (shorted) { 493 493 + Vnodes[j] = 0; 494 494 + } else { 495 495 + Vnodes[j] = (I + Is) / (G + Gs); 496 496 + } 497 497 + updateShuntStates(nodeShunts[j], nodeShuntState[j], Vnodes[j]); 498 498 + } 499 499 + 500 500 + // Load node (j=N): segment N-1 on left; terminal on right. 501 501 + { 502 502 + const aL = vArr.right[N - 1]; // V⁺ arriving at right end of last seg 503 503 + const Z0 = segs[N - 1].Z0; 504 504 + const term = model.terminal; 505 505 + // First apply any mid-node shunts at the load-side node (nodeShunts[N]). 506 506 + let G = 1 / Z0; 507 507 + let I = 2 * aL / Z0; 508 508 + const { G: Gs, I: Is, shorted: sh0 } = accumulateShunts( 509 509 + nodeShunts[N], nodeShuntState[N]); 510 510 + if (sh0) { 511 511 + Vnodes[N] = 0; 512 512 + } else { 513 513 + // Then add the terminal element. 514 514 + let Gt = 0, It = 0, shorted = false; 515 515 + if (term.type === 'open') { 516 516 + // no terminal conductance 517 517 + } else if (term.type === 'short') { 518 518 + shorted = true; 519 519 + } else if (term.type === 'R') { 520 520 + Gt = 1 / term.value; 521 521 + } else if (term.type === 'C') { 522 522 + Gt = term.value / dt; 523 523 + It = Gt * termState.vCap; 524 524 + } else if (term.type === 'L') { 525 525 + Gt = dt / term.value; 526 526 + It = -termState.iL; 527 527 + } 528 528 + Vnodes[N] = shorted ? 0 : (I + Is + It) / (G + Gs + Gt); 529 529 + } 530 530 + updateShuntStates(nodeShunts[N], nodeShuntState[N], Vnodes[N]); 531 531 + if (term.type === 'C') termState.vCap = Vnodes[N]; 532 532 + if (term.type === 'L') termState.iL += (dt / term.value) * Vnodes[N]; 533 533 + } 534 534 + 535 535 + // 3. Compute outgoing waves and write into delay lines. 536 536 + // At node j: V = V⁺_out + V⁻_in → V⁺_out = V - V⁻_in (outgoing into seg j rightward) 537 537 + // At node j: V = V⁺_in + V⁻_out → V⁻_out = V - V⁺_in (outgoing into seg j-1 leftward) 538 538 + for (let i = 0; i < N; i++) { 539 539 + const vOutFwd = Vnodes[i] - vArr.left[i]; // V⁺ leaving node i into seg i 540 540 + const vOutBwd = Vnodes[i + 1] - vArr.right[i]; // V⁻ leaving node i+1 into seg i 541 541 + vPlusHist[i][step] = vOutFwd; 542 542 + vMinusHist[i][step] = vOutBwd; 543 543 + delayFwd[i][fwdIdx[i]] = vOutFwd; 544 544 + delayBwd[i][bwdIdx[i]] = vOutBwd; 545 545 + fwdIdx[i] = (fwdIdx[i] + 1) % nDelay[i]; 546 546 + bwdIdx[i] = (bwdIdx[i] + 1) % nDelay[i]; 547 547 + } 548 548 + 549 549 + // 4. Record node voltages. 550 550 + for (let j = 0; j <= N; j++) nodeV[j][step] = Vnodes[j]; 551 551 + } 552 552 + 553 553 + const tGrid = Float64Array.from({ length: nSteps }, (_, k) => k * dt); 554 554 + return { dt, nSteps, tGrid, segs, nDelay, nodeV, vPlusHist, vMinusHist }; 555 555 + } 556 556 + 557 557 + // ---- Spatial voltage query for simulateTimeDomain output ---- 558 558 + // 559 559 + // Returns V(z, tIdx) by reading from vPlusHist / vMinusHist at the appropriate 560 560 + // delay offsets. z ∈ [0, 1] is normalized by total T-line delay (sum of taus). 561 561 + // 562 562 + // Within segment i, at fractional position f ∈ [0, 1]: 563 563 + // V⁺ that is at position f at step tIdx was written into the segment 564 564 + // kFwd = round(f * D) steps ago → look up vPlusHist[i][tIdx - kFwd] 565 565 + // V⁻ that is at position f at step tIdx was written into the segment 566 566 + // kBwd = round((1-f) * D) steps ago → look up vMinusHist[i][tIdx - kBwd] 567 567 + function voltageAt(z, tIdx, sim) { 568 568 + const { segs, nDelay, vPlusHist, vMinusHist, nodeV } = sim; 569 569 + const N = segs.length; 570 570 + const totalTau = segs.reduce((s, seg) => s + seg.tau, 0); 571 571 + let cumZ = 0; 572 572 + for (let i = 0; i < N; i++) { 573 573 + const segZ = segs[i].tau / totalTau; 574 574 + if (z <= cumZ + segZ + 1e-9) { 575 575 + const f = Math.min(1, Math.max(0, (z - cumZ) / segZ)); 576 576 + const D = nDelay[i]; 577 577 + const kFwd = Math.round(f * D); 578 578 + const kBwd = D - kFwd; 579 579 + const tF = tIdx - kFwd; 580 580 + const tB = tIdx - kBwd; 581 581 + const vP = tF >= 0 ? vPlusHist[i][tF] : 0; 582 582 + const vM = tB >= 0 ? vMinusHist[i][tB] : 0; 583 583 + return vP + vM; 584 584 + } 585 585 + cumZ += segZ; 586 586 + } 587 587 + // z at or past the load end: return load node voltage 588 588 + return nodeV[N][tIdx]; 589 589 + } 590 590 + 591 591 + // ---- Convenience: convert old-style model (model.segments, model.RL) to blocks ---- 592 592 + // This lets existing callers construct a blocks-based model without rewriting app.js. 593 593 + function segmentsToBlocks(model) { 594 594 + const N = model.segments.length; 595 595 + const tau = 1 / N; // each segment = 1/N of total normalized delay 596 596 + const blocks = model.segments.map(seg => ({ type: 'tl', Z0: seg.Z0, tau })); 597 597 + const RL = model.RL; 598 598 + const terminal = isFinite(RL) 599 599 + ? { type: 'R', value: RL } 600 600 + : { type: 'open' }; 601 601 + return { ...model, blocks, terminal }; 602 602 + } 603 603 + 314 604 return { 315 605 BOUNCE_EPS, 316 606 MAX_BOUNCES, ··· 325 615 totalVoltageAt, 326 616 segmentForZ, 327 617 computeDynamicState, 618 618 + simulateTimeDomain, 619 619 + voltageAt, 620 620 + segmentsToBlocks, 328 621 }; 329 622 })();

+89

physics.test.js

reviewed

··· 346 346 spiceCompare("short circuit (tline-sc.sp)", "results-tline-sc.tsv", 0); 347 347 // RL=100 (resistive load) ΓL=+1/3, ΓS=−2/3 348 348 spiceCompare("resistive load 100Ω (tline-rl.sp)", "results-tline-rl.tsv", 100); 349 349 + 350 350 + // ──────────────────────────────────────────────────────────────────────────── 351 351 + // simulateTimeDomain: validate against buildBounceSeries (resistive cases) 352 352 + // 353 353 + // Strategy: compare node voltages at the source (node 0) and load (node N) at 354 354 + // several tNorm values. We sample at half-integer multiples of τ_d so that 355 355 + // we land between bounce events, where both models agree on a constant value. 356 356 + // With integer delay steps the MoC simulation is exact for resistive circuits. 357 357 + // ──────────────────────────────────────────────────────────────────────────── 358 358 + 359 359 + const { simulateTimeDomain, voltageAt, segmentsToBlocks } = globalThis.TLPhysics; 360 360 + 361 361 + // Helper: build a blocks-format model from old-style params and run the sim. 362 362 + function makeSim(Vg, Rg, Z0, RL, tEnd = 12, oversample = 16) { 363 363 + const oldModel = { Vg, Rg, RL, segments: [{ Z0 }], reflectTol: 0.001 }; 364 364 + const newModel = segmentsToBlocks(oldModel); 365 365 + return { sim: simulateTimeDomain(newModel, { tEnd, oversample }), oldModel }; 366 366 + } 367 367 + 368 368 + // Compare VS (node 0) and VL (node N) between MoC sim and bounce series. 369 369 + function checkSimVsBounce(label, Vg, Rg, Z0, RL, tNorms, tol = 1e-9) { 370 370 + test(`simulateTimeDomain vs bounce: ${label}`, () => { 371 371 + const { sim, oldModel } = makeSim(Vg, Rg, Z0, RL); 372 372 + const waves = computeWaveParams(oldModel); 373 373 + const bounce = buildBounceSeries(oldModel, waves); 374 374 + const { dt, nSteps, nodeV } = sim; 375 375 + for (const tn of tNorms) { 376 376 + const step = Math.round(tn / dt); 377 377 + if (step >= nSteps) continue; 378 378 + near(nodeV[0][step], sumEventsAtTime(bounce.srcEvents, tn), `${label} VS @ tNorm=${tn}`, tol); 379 379 + near(nodeV[1][step], sumEventsAtTime(bounce.loadEvents, tn), `${label} VL @ tNorm=${tn}`, tol); 380 380 + } 381 381 + }); 382 382 + } 383 383 + 384 384 + // Sample at half-integer τ values: between bounce events both models are constant. 385 385 + const tSamples = [0.5, 1.5, 2.5, 3.5, 4.5, 5.5]; 386 386 + 387 387 + checkSimVsBounce("matched (50/50/50)", 1, 50, 50, 50, tSamples); 388 388 + checkSimVsBounce("open circuit (50/50/∞)", 1, 50, 50, Infinity, tSamples); 389 389 + checkSimVsBounce("short circuit (50/50/0)", 1, 50, 50, 0, tSamples); 390 390 + checkSimVsBounce("mismatched Rg=100 RL=200", 1, 100, 50, 200, tSamples); 391 391 + checkSimVsBounce("Vg=5 Rg=20 Z0=50 RL=100", 5, 20, 50, 100, tSamples); 392 392 + 393 393 + test("simulateTimeDomain: 2 equal segments match bounce series VL", () => { 394 394 + const oldModel = { Vg: 1, Rg: 50, RL: 100, segments: [{ Z0: 50 }, { Z0: 50 }], reflectTol: 0.001 }; 395 395 + const newModel = segmentsToBlocks(oldModel); 396 396 + const sim = simulateTimeDomain(newModel, { tEnd: 12, oversample: 16 }); 397 397 + const bounce = buildBounceSeries(oldModel, computeWaveParams(oldModel)); 398 398 + const { dt, nSteps, nodeV } = sim; 399 399 + for (const tn of tSamples) { 400 400 + const step = Math.round(tn / dt); 401 401 + if (step >= nSteps) continue; 402 402 + near(nodeV[2][step], sumEventsAtTime(bounce.loadEvents, tn), `VL @ tNorm=${tn}`, 1e-9); 403 403 + } 404 404 + }); 405 405 + 406 406 + test("simulateTimeDomain: 2 segments Z0=50/75, DC = Vg·RL/(Rg+RL)", () => { 407 407 + const oldModel = { Vg: 1, Rg: 50, RL: 100, segments: [{ Z0: 50 }, { Z0: 75 }], reflectTol: 0.001 }; 408 408 + const newModel = segmentsToBlocks(oldModel); 409 409 + const sim = simulateTimeDomain(newModel, { tEnd: 20, oversample: 8 }); 410 410 + const step = Math.min(Math.round(18 / sim.dt), sim.nSteps - 1); 411 411 + near(sim.nodeV[2][step], 1 * 100 / (50 + 100), "VL DC steady state", 1e-3); 412 412 + }); 413 413 + 414 414 + test("voltageAt: z=0 and z=1 match source/load node voltages", () => { 415 415 + const { sim } = makeSim(1, 50, 50, Infinity); 416 416 + const { dt, nSteps } = sim; 417 417 + for (const tn of [0.5, 1.5, 2.5]) { 418 418 + const step = Math.round(tn / dt); 419 419 + if (step >= nSteps) continue; 420 420 + near(voltageAt(0, step, sim), sim.nodeV[0][step], `z=0 @ tNorm=${tn}`, 1e-9); 421 421 + near(voltageAt(1, step, sim), sim.nodeV[1][step], `z=1 @ tNorm=${tn}`, 1e-9); 422 422 + } 423 423 + }); 424 424 + 425 425 + test("voltageAt: midpoint of matched line = V1 once front has passed (tNorm=0.7)", () => { 426 426 + // Matched line Rg=Z0=RL=50: single rightward wave V1=0.5, no reflections. 427 427 + // At tNorm=0.7 the front is past z=0.5, so V(0.5) = V1 = 0.5. 428 428 + const { sim } = makeSim(1, 50, 50, 50, 12, 16); 429 429 + const step = Math.round(0.7 / sim.dt); 430 430 + near(voltageAt(0.5, step, sim), 0.5, "V(0.5) = 0.5", 1e-9); 431 431 + }); 432 432 + 433 433 + test("voltageAt: midpoint = 0 before wave arrives (tNorm=0.3)", () => { 434 434 + const { sim } = makeSim(1, 50, 50, 50, 12, 16); 435 435 + const step = Math.round(0.3 / sim.dt); 436 436 + near(voltageAt(0.5, step, sim), 0, "V(0.5) = 0 before wave", 1e-9); 437 437 + });