Make erf rise shape causal by shifting to [0, tr] interval

+7 -8

app.js

··· 42 42 RL: 30, 43 43 secPerTau: 2.5, 44 44 reflectTol: 1, 45 - riseTimeTau: 0, 46 45 riseTimeTr: 0, 46 + riseShape: "step", // "step" | "linear" | "erf" 47 47 }; 48 48 49 49 // ---- segment input management ---- ··· 92 92 model.RL = rlIsOpen ? Infinity : parseFloat(el.RL.value); 93 93 model.secPerTau = parseFloat(el.secPerTau.value); 94 94 model.reflectTol = Math.max(0, parseFloat(el.reflectTol.value)); 95 - model.riseTimeTau = el.riseMode.value === "exp" 96 - ? Math.max(0.001, parseFloat(el.riseTau.value) || 0.1) 97 - : 0; 98 - model.riseTimeTr = el.riseMode.value === "linear" 95 + model.riseShape = el.riseMode.value === "step" ? "step" 96 + : el.riseMode.value === "linear" ? "linear" : "erf"; 97 + model.riseTimeTr = model.riseShape !== "step" 99 98 ? Math.max(0.001, parseFloat(el.riseTau.value) || 0.1) 100 99 : 0; 101 100 model.segments = segZ0Inputs.map((inp) => ({ ··· 157 156 syncModelFromInputs(); 158 157 const waves = computeWaveParams(model); 159 158 const bounce = buildBounceSeries(model, waves); 160 - const dyn = computeDynamicState(tNorm, bounce, model.riseTimeTau, model.riseTimeTr); 159 + const dyn = computeDynamicState(tNorm, bounce, model.riseTimeTr, model.riseShape); 161 160 timeHorizon = Math.max(2.2, Math.min(16, bounce.tEnd + 0.4)); 162 161 163 162 el.secPerTauRead.textContent = fmt(model.secPerTau, 1); ··· 173 172 } 174 173 175 174 const d = resizeCanvasToCSS(el.tdr); 176 - drawTDR(d.ctx, d.w, d.h, tNorm, bounce, model.riseTimeTau, model.riseTimeTr, timeHorizon, theme); 175 + drawTDR(d.ctx, d.w, d.h, tNorm, bounce, model.riseTimeTr, model.riseShape, timeHorizon, theme); 177 176 178 177 const c = resizeCanvasToCSS(el.circuit); 179 178 drawCircuit(c.ctx, c.w, c.h, tNorm, dyn, theme, model.segments, model.RL); 180 179 181 180 ensurePlotCanvasHeight(el.plot, 2); 182 181 const p = resizeCanvasToCSS(el.plot); 183 - drawPlot(p.ctx, p.w, p.h, tNorm, dyn, theme, model.riseTimeTau, model.riseTimeTr, model.segments); 182 + drawPlot(p.ctx, p.w, p.h, tNorm, dyn, theme, model.riseTimeTr, model.riseShape, model.segments); 184 183 185 184 if (!hasStarted) { 186 185 el.pauseBtn.textContent = "Pause";

+1 -1

index.html

··· 45 45 <select id="riseMode"> 46 46 <option value="step">Step (instantaneous)</option> 47 47 <option value="linear">Linear ramp (trapezoidal)</option> 48 - <option value="exp">Gaussian rise (skin-effect)</option> 48 + <option value="exp">Smooth S-curve (erf)</option> 49 49 </select> 50 50 </label> 51 51 <label><span class="lbl">Rise time constant τ<sub>r</sub> / τ<sub>d</sub></span>

+37 -41

physics.js

··· 30 30 // front — current position of the leading edge ∈ [zStart, zEnd] (or reversed for leftward) 31 31 // 32 32 // Extension points: 33 - // • Rise time (skin-effect / Gaussian erf): model.riseTimeTau > 0 — see riseShape/waveVoltageAt. 33 + // • Rise time (causal shifted erf or linear ramp): model.riseTimeTr > 0 — see riseShape/waveVoltageAt. 34 34 // • Non-equal segment lengths: store lengths[] summing to 1; the BFS solver already works 35 35 // by position, so changing zBnds is all that's needed. 36 36 ··· 190 190 return sign * y; 191 191 } 192 192 193 - // Gaussian (erf): tau = 0 → hard step. tau > 0 → 0.5·(1 + erf(dt/tau)). 194 - // Models skin-effect losses (α ∝ √f) whose step response is an error function. 195 - // Centered at dt=0 (50% at front arrival), with symmetric smooth sigmoid — 196 - // zero slope at dt=0 from both sides, small precursor for dt < 0. 197 - function riseShape(dt, tau) { 198 - if (tau <= 0) return dt <= 0 ? 0 : 1; 199 - return 0.5 * (1 + erf(dt / tau)); 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))); 200 202 } 201 203 202 204 // Linear ramp (trapezoidal / SPICE PULSE): tr = 0 → hard step. tr > 0 → dt/tr clamped to [0,1]. ··· 209 211 // Voltage contribution of a single wave packet wf at position z. 210 212 // wf must carry the fields added by computeDynamicState: u, front. 211 213 // Returns 0 if z is outside this packet's z-range or ahead of the front. 212 - // Pass either tau (erf/Gaussian) or tr (linear); whichever is > 0 takes priority (tr first). 213 - function waveVoltageAt(wf, z, tau, tr = 0) { 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") { 214 217 const { dir, zStart, A, u, front } = wf; 215 - // For the centered erf, the wave has a small precursor ahead of the front. 216 - // Extend the spatial window by 3·tau (erf(3)≈0.9999) to capture it. 217 - const margin = (tr > 0 || tau <= 0) ? 1e-9 : 3 * tau; 218 + const eps = 1e-9; 219 + const rise = shape === "linear" ? riseShapeLinear : riseShape; 218 220 if (dir > 0) { 219 - // Rightward: nonzero in [zStart, front + margin]. 220 - // Time since front passed z: u − (z − zStart). 221 - if (z < zStart - 1e-9 || z > front + margin) return 0; 221 + // Rightward: nonzero in [zStart, front]. 222 + if (z < zStart - eps || z > front + eps) return 0; 222 223 const dt = u - (z - zStart); 223 - return A * (tr > 0 ? riseShapeLinear(dt, tr) : riseShape(dt, tau)); 224 + return A * rise(dt, tr); 224 225 } else { 225 - // Leftward: nonzero in [front - margin, zStart]. 226 - // Time since front passed z: u − (zStart − z). 227 - if (z > zStart + 1e-9 || z < front - margin) return 0; 226 + // Leftward: nonzero in [front, zStart]. 227 + if (z > zStart + eps || z < front - eps) return 0; 228 228 const dt = u - (zStart - z); 229 - return A * (tr > 0 ? riseShapeLinear(dt, tr) : riseShape(dt, tau)); 229 + return A * rise(dt, tr); 230 230 } 231 231 } 232 232 ··· 242 242 // For multi-segment models (N > 1), only waves whose segIdx matches the segment 243 243 // that owns z are summed. This prevents double-counting at segment boundaries 244 244 // where adjacent waves share a z value (a plotting artifact, not a physics error). 245 - function totalVoltageAt(z, launchedWaves, tau, tr = 0, N = 1) { 245 + function totalVoltageAt(z, launchedWaves, tr = 0, shape = "erf", N = 1) { 246 246 const segIdx = N > 1 ? segmentForZ(z, N) : 0; 247 247 let V = 0; 248 248 for (const wf of launchedWaves) { 249 249 if (N > 1 && wf.segIdx !== segIdx) continue; 250 - V += waveVoltageAt(wf, z, tau, tr); 250 + V += waveVoltageAt(wf, z, tr, shape); 251 251 } 252 252 return V; 253 253 } 254 254 255 - // Node voltage as smooth sum of erf-rise events. 256 - // The centered erf has a precursor, so events slightly in the future contribute. 257 - // We include events up to 3·tau ahead (erf(−3)≈−0.9999, giving ~0.02% amplitude). 258 - function sumEventsWithRise(events, tn, tau) { 255 + // Node voltage as smooth sum of causal erf-rise events. 256 + function sumEventsWithRise(events, tn, tr) { 259 257 let v = 0; 260 - const lookAhead = 3 * tau; 261 258 for (const e of events) { 262 - if (tn >= e.t - lookAhead) v += e.dV * riseShape(tn - e.t, tau); 259 + const dt = tn - e.t; 260 + if (dt > 0) v += e.dV * riseShape(dt, tr); 263 261 } 264 262 return v; 265 263 } ··· 279 277 // u: tNorm − tBorn 280 278 // front: current leading-edge position 281 279 // A wave is "active" while u < |zEnd − zStart| (front still in transit). 282 - function computeDynamicState(tn, bounce, riseTimeTau = 0, riseTimeTr = 0) { 280 + // shape: "step" | "linear" | "erf" 281 + function computeDynamicState(tn, bounce, riseTimeTr = 0, riseShape_ = "erf") { 283 282 const { clamp } = TLUtils; 284 283 const launchedWaves = []; 285 284 const activeWaves = []; ··· 291 290 if (w.dir > 0) { 292 291 front = clamp(w.zStart + u, w.zStart, w.zEnd); 293 292 } else { 294 - // zEnd < zStart for leftward 295 293 front = clamp(w.zStart - u, w.zEnd, w.zStart); 296 294 } 297 295 const ww = { ...w, u, front }; ··· 299 297 if (u < Math.abs(w.zEnd - w.zStart) - 1e-9) activeWaves.push(ww); 300 298 } 301 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 + 302 306 return { 303 307 launchedWaves, 304 308 activeWaves, 305 - VS: riseTimeTr > 0 306 - ? sumEventsWithLinearRamp(bounce.srcEvents, tn, riseTimeTr) 307 - : riseTimeTau > 0 308 - ? sumEventsWithRise(bounce.srcEvents, tn, riseTimeTau) 309 - : sumEventsAtTime(bounce.srcEvents, tn), 310 - VL: riseTimeTr > 0 311 - ? sumEventsWithLinearRamp(bounce.loadEvents, tn, riseTimeTr) 312 - : riseTimeTau > 0 313 - ? sumEventsWithRise(bounce.loadEvents, tn, riseTimeTau) 314 - : sumEventsAtTime(bounce.loadEvents, tn), 309 + VS: sumNode(bounce.srcEvents), 310 + VL: sumNode(bounce.loadEvents), 315 311 }; 316 312 } 317 313

+29 -24

physics.test.js

··· 186 186 // riseShape 187 187 // ──────────────────────────────────────────────────────────────────────────── 188 188 189 - test("riseShape: tau=0 (step) — hard step at dt=0", () => { 190 - near(riseShape(0, 0), 0, "dt=0 tau=0"); 191 - near(riseShape(-1, 0), 0, "dt<0 tau=0"); 192 - near(riseShape(0.001, 0), 1, "tiny dt, tau=0"); 193 - near(riseShape(100, 0), 1, "large dt, tau=0"); 189 + test("riseShape: tr=0 (step) — hard step at dt=0", () => { 190 + near(riseShape(0, 0), 0, "dt=0 tr=0"); 191 + near(riseShape(-1, 0), 0, "dt<0 tr=0"); 192 + near(riseShape(0.001, 0), 1, "tiny dt, tr=0"); 193 + near(riseShape(100, 0), 1, "large dt, tr=0"); 194 194 }); 195 195 196 - test("riseShape: tau>0 — centered erf: 0.5·(1+erf(dt/tau))", () => { 197 - // At dt=0: 0.5 (midpoint of sigmoid); erf(0) approximation has ~1e-9 error 198 - near(riseShape(0, 1), 0.5, "dt=0 → 0.5", 1e-7); 199 - // 0.5·(1+erf(1)) ≈ 0.9214 200 - near(riseShape(1, 1), 0.5 * (1 + 0.8427007929497149), "dt=tau", 1.5e-7); 201 - // 0.5·(1+erf(2)) ≈ 0.9977 202 - near(riseShape(2, 1), 0.5 * (1 + 0.9953222650189527), "dt=2·tau", 1.5e-7); 203 - // Negative dt: 0.5·(1+erf(−1)) ≈ 0.0786 (precursor) 204 - near(riseShape(-1, 1), 0.5 * (1 - 0.8427007929497149), "dt=-tau", 1.5e-7); 205 - // Symmetry: riseShape(dt) + riseShape(-dt) = 1 206 - near(riseShape(0.7, 1) + riseShape(-0.7, 1), 1, "symmetry", 1e-6); 196 + test("riseShape: tr>0 — causal shifted erf over [0, tr]", () => { 197 + // Causal: dt <= 0 → exactly 0 198 + near(riseShape(0, 1), 0, "dt=0 → 0", 1e-3); 199 + near(riseShape(-1, 1), 0, "dt<0 → 0"); 200 + // dt = tr → essentially 1 201 + near(riseShape(1, 1), 1, "dt=tr → 1", 1e-3); 202 + // Midpoint dt = tr/2 → 0.5 (erf(0) = 0) 203 + near(riseShape(0.5, 1), 0.5, "dt=tr/2 → 0.5", 1e-7); 204 + // Monotonically increasing 205 + const v1 = riseShape(0.2, 1); 206 + const v2 = riseShape(0.5, 1); 207 + const v3 = riseShape(0.8, 1); 208 + assert(v1 < v2, `monotonic: r(0.2)=${v1} < r(0.5)=${v2}`); 209 + assert(v2 < v3, `monotonic: r(0.5)=${v2} < r(0.8)=${v3}`); 210 + // Clamped to 1 for dt > tr 211 + near(riseShape(2, 1), 1, "dt>tr → 1"); 207 212 }); 208 213 209 214 // ──────────────────────────────────────────────────────────────────────────── ··· 254 259 segments: [{ Z0: 40 }, { Z0: 60 }], reflectTol: 0.001 }; 255 260 const waves = computeWaveParams(model); 256 261 const bounce = buildBounceSeries(model, waves); 257 - const dyn = computeDynamicState(0.6, bounce, 0, 0); 262 + const dyn = computeDynamicState(0.6, bounce); 258 263 const N = model.segments.length; 259 264 const eps = 1e-4; 260 265 261 - const vL = totalVoltageAt(0.5 - eps, dyn.launchedWaves, 0, 0, N); 262 - const vB = totalVoltageAt(0.5, dyn.launchedWaves, 0, 0, N); 263 - const vR = totalVoltageAt(0.5 + eps, dyn.launchedWaves, 0, 0, N); 266 + const vL = totalVoltageAt(0.5 - eps, dyn.launchedWaves, 0, "step", N); 267 + const vB = totalVoltageAt(0.5, dyn.launchedWaves, 0, "step", N); 268 + const vR = totalVoltageAt(0.5 + eps, dyn.launchedWaves, 0, "step", N); 264 269 265 270 // Correct value: (1 + 0.2) * (4/9) = 8/15 266 271 near(vB, 8 / 15, "V at boundary = (1+Γ)·V1", 1e-6); ··· 280 285 const eps = 1e-4; 281 286 282 287 for (const tNorm of [0.3, 0.6, 1.2, 2.0]) { 283 - const dyn = computeDynamicState(tNorm, bounce, 0, 0); 288 + const dyn = computeDynamicState(tNorm, bounce); 284 289 for (const zB of [0.25, 0.5, 0.75]) { 285 - const vL = totalVoltageAt(zB - eps, dyn.launchedWaves, 0, 0, N); 286 - const vB = totalVoltageAt(zB, dyn.launchedWaves, 0, 0, N); 287 - const vR = totalVoltageAt(zB + eps, dyn.launchedWaves, 0, 0, N); 290 + const vL = totalVoltageAt(zB - eps, dyn.launchedWaves, 0, "step", N); 291 + const vB = totalVoltageAt(zB, dyn.launchedWaves, 0, "step", N); 292 + const vR = totalVoltageAt(zB + eps, dyn.launchedWaves, 0, "step", N); 288 293 // No spike: boundary value must lie between its two neighbours (within tolerance). 289 294 const lo = Math.min(vL, vR) - 1e-6; 290 295 const hi = Math.max(vL, vR) + 1e-6;

+9 -9

render.js

··· 333 333 334 334 // ---- plot canvas ---- 335 335 // segments = model.segments (used to draw boundary tick marks) 336 - function drawPlot(ctx, w, h, tn, dyn, theme, riseTimeTau = 0, riseTimeTr = 0, segments = [{ Z0: 50 }]) { 336 + function drawPlot(ctx, w, h, tn, dyn, theme, riseTimeTr = 0, riseShape = "step", segments = [{ Z0: 50 }]) { 337 337 ctx.clearRect(0, 0, w, h); 338 338 ctx.fillStyle = theme.panel; 339 339 ctx.fillRect(0, 0, w, h); ··· 439 439 440 440 const panelH = PLOT_PANEL_H; 441 441 const { waveVoltageAt, totalVoltageAt } = TLPhysics; 442 - const smooth = riseTimeTau > 0 || riseTimeTr > 0; 442 + const smooth = riseTimeTr > 0 && riseShape !== "step"; 443 443 444 444 // Panel 1 — sum of all waves 445 445 const top0 = PLOT_PAD_T; ··· 447 447 const y0 = (V) => top0 + 4 + (vHi - V) / (vHi - vLo) * (panelH - 8); 448 448 drawPanelFrame(top0, bot0, y0, "Sum (all waves)", vLo, vHi); 449 449 if (smooth) { 450 - drawSampledWave(ctx, xOfZ, y0, (z) => totalVoltageAt(z, launched, riseTimeTau, riseTimeTr, segments.length), theme.ok, 2.4); 450 + drawSampledWave(ctx, xOfZ, y0, (z) => totalVoltageAt(z, launched, riseTimeTr, riseShape, segments.length), theme.ok, 2.4); 451 451 } else { 452 452 drawPiecewise(ctx, xOfZ, y0, sumSegs, theme.ok, 2.4); 453 453 } ··· 483 483 ctx.globalAlpha = alpha; 484 484 ctx.setLineDash(style.dash); 485 485 if (smooth) { 486 - drawSampledWave(ctx, xOfZ, y1, (z) => waveVoltageAt(wf, z, riseTimeTau, riseTimeTr), style.color, 2.0); 486 + drawSampledWave(ctx, xOfZ, y1, (z) => waveVoltageAt(wf, z, riseTimeTr, riseShape), style.color, 2.0); 487 487 } else { 488 488 drawPiecewise(ctx, xOfZ, y1, segmentsForWave(wf), style.color, 2.0); 489 489 } ··· 502 502 // ---- TDR canvas ---- 503 503 // Draws VS(t) — voltage at the detection point (z=0, between Rg and T-line) versus time. 504 504 // The trace grows up to the current tNorm; a cursor line shows current time. 505 - function drawTDR(ctx, w, h, tn, bounce, riseTimeTau, riseTimeTr, timeHorizon, theme) { 505 + function drawTDR(ctx, w, h, tn, bounce, riseTimeTr, riseShape, timeHorizon, theme) { 506 506 ctx.clearRect(0, 0, w, h); 507 507 ctx.fillStyle = theme.panel; 508 508 ctx.fillRect(0, 0, w, h); ··· 516 516 517 517 const { sumEventsAtTime, sumEventsWithRise, sumEventsWithLinearRamp } = TLPhysics; 518 518 const vsAtTime = (t) => { 519 - if (riseTimeTr > 0) return sumEventsWithLinearRamp(bounce.srcEvents, t, riseTimeTr); 520 - if (riseTimeTau > 0) return sumEventsWithRise(bounce.srcEvents, t, riseTimeTau); 521 - return sumEventsAtTime(bounce.srcEvents, t); 519 + if (riseTimeTr <= 0 || riseShape === "step") return sumEventsAtTime(bounce.srcEvents, t); 520 + if (riseShape === "linear") return sumEventsWithLinearRamp(bounce.srcEvents, t, riseTimeTr); 521 + return sumEventsWithRise(bounce.srcEvents, t, riseTimeTr); 522 522 }; 523 523 524 524 // y range from all srcEvents (cumulative sum) ··· 609 609 610 610 // VS(t) trace up to tn 611 611 const tnClamped = Math.min(tn, tMax); 612 - const smooth = riseTimeTau > 0 || riseTimeTr > 0; 612 + const smooth = riseTimeTr > 0 && riseShape !== "step"; 613 613 ctx.strokeStyle = theme.tdr; 614 614 ctx.lineWidth = 2.2; 615 615 ctx.beginPath();

Configure Feed

Configure Feed