this repo has no description
Replace exponential rise shape with centered Gaussian erf
The first-order RC exponential (1-exp(-dt/tau)) had infinite slope at
dt=0 and a lazy asymptotic tail — not physically realistic. Replace with
0.5*(1+erf(dt/tau)), which models skin-effect losses (alpha proportional
to sqrt(f)) and gives zero slope at dt=0 with a smooth symmetric sigmoid.
- Add erf() rational approximation (Abramowitz & Stegun 7.1.26)
- Extend waveVoltageAt spatial window by 3*tau for precursor foot
- Update sumEventsWithRise to include events up to 3*tau ahead
- Rename UI option to "Gaussian rise (skin-effect)"
Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
+44
-27
+1
-1
index.html
+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">Exponential rise (RC)</option>
48
+
<option value="exp">Gaussian rise (skin-effect)</option>
49
49
</select>
50
50
</label>
51
51
<label><span class="lbl">Rise time constant τ<sub>r</sub> / τ<sub>d</sub></span>
+29
-12
physics.js
+29
-12
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 (RC/finite bandwidth): model.riseTimeTau > 0 — see riseShape/waveVoltageAt.
33
+
// • Rise time (skin-effect / Gaussian erf): model.riseTimeTau > 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
···
181
181
// ---- rise-time wave shapes ----
182
182
// Fraction of final amplitude reached after time dt since the wave front passed.
183
183
184
-
// Exponential: tau = 0 → hard step. tau > 0 → 1 − exp(−dt/tau).
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
+
// 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.
185
197
function riseShape(dt, tau) {
186
-
if (dt <= 0) return 0;
187
-
if (tau <= 0) return 1;
188
-
return 1 - Math.exp(-dt / tau);
198
+
if (tau <= 0) return dt <= 0 ? 0 : 1;
199
+
return 0.5 * (1 + erf(dt / tau));
189
200
}
190
201
191
202
// Linear ramp (trapezoidal / SPICE PULSE): tr = 0 → hard step. tr > 0 → dt/tr clamped to [0,1].
···
198
209
// Voltage contribution of a single wave packet wf at position z.
199
210
// wf must carry the fields added by computeDynamicState: u, front.
200
211
// Returns 0 if z is outside this packet's z-range or ahead of the front.
201
-
// Pass either tau (exponential) or tr (linear); whichever is > 0 takes priority (tr first).
212
+
// Pass either tau (erf/Gaussian) or tr (linear); whichever is > 0 takes priority (tr first).
202
213
function waveVoltageAt(wf, z, tau, tr = 0) {
203
214
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;
204
218
if (dir > 0) {
205
-
// Rightward: nonzero in [zStart, front].
219
+
// Rightward: nonzero in [zStart, front + margin].
206
220
// Time since front passed z: u − (z − zStart).
207
-
if (z < zStart - 1e-9 || z > front + 1e-9) return 0;
221
+
if (z < zStart - 1e-9 || z > front + margin) return 0;
208
222
const dt = u - (z - zStart);
209
223
return A * (tr > 0 ? riseShapeLinear(dt, tr) : riseShape(dt, tau));
210
224
} else {
211
-
// Leftward: nonzero in [front, zStart].
225
+
// Leftward: nonzero in [front - margin, zStart].
212
226
// Time since front passed z: u − (zStart − z).
213
-
if (z > zStart + 1e-9 || z < front - 1e-9) return 0;
227
+
if (z > zStart + 1e-9 || z < front - margin) return 0;
214
228
const dt = u - (zStart - z);
215
229
return A * (tr > 0 ? riseShapeLinear(dt, tr) : riseShape(dt, tau));
216
230
}
···
238
252
return V;
239
253
}
240
254
241
-
// Node voltage as smooth sum of exponential-rise events.
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).
242
258
function sumEventsWithRise(events, tn, tau) {
243
259
let v = 0;
260
+
const lookAhead = 3 * tau;
244
261
for (const e of events) {
245
-
if (tn >= e.t) v += e.dV * riseShape(tn - e.t, tau);
262
+
if (tn >= e.t - lookAhead) v += e.dV * riseShape(tn - e.t, tau);
246
263
}
247
264
return v;
248
265
}
+14
-14
physics.test.js
+14
-14
physics.test.js
···
186
186
// riseShape
187
187
// ────────────────────────────────────────────────────────────────────────────
188
188
189
-
test("riseShape: dt≤0 → 0", () => {
190
-
near(riseShape(0, 0), 0, "dt=0 tau=0");
191
-
near(riseShape(-1, 1), 0, "dt<0");
192
-
});
193
-
194
-
test("riseShape: tau=0 (step) → 1 for dt>0", () => {
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");
195
192
near(riseShape(0.001, 0), 1, "tiny dt, tau=0");
196
193
near(riseShape(100, 0), 1, "large dt, tau=0");
197
194
});
198
195
199
-
test("riseShape: tau>0 — known exponential values", () => {
200
-
// riseShape(τ, τ) = 1 − e^{−1} ≈ 0.63212
201
-
near(riseShape(1, 1), 1 - Math.exp(-1), "1 time constant", 1e-12);
202
-
// riseShape(2τ, τ) = 1 − e^{−2} ≈ 0.86466
203
-
near(riseShape(2, 1), 1 - Math.exp(-2), "2 time constants", 1e-12);
204
-
// Small dt → roughly dt/tau (first-order approx)
205
-
const small = riseShape(0.001, 1);
206
-
near(small, 0.001, "linear regime approx", 1e-5);
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);
207
207
});
208
208
209
209
// ────────────────────────────────────────────────────────────────────────────