| Waveform | Qpos (nC) | Qneg (nC) | W (µJ) | Ipk,pos (mA) | Ipk,neg (mA) |
|---|---|---|---|---|---|
| Check one or more waveforms, then click "Run simulation". | |||||
| Waveform | Ipk/Q (mA/nC) | Wfar/Wcap | EV,out (%) | EI,out (%) | Regime |
|---|---|---|---|---|---|
| Advanced metrics appear here after a run. | |||||
Across the 64 tissue × waveform combinations of the reference benchmark, this simulator agrees with the ngspice reference to within a few percent on 62 of them, and to better than 0.1 % on the low-impedance tissues (muscle, blood, CSF, heart). Two cases sit outside that band: Skin × rectangular_asymmetric on the residual charge Qpos (~85 %) and its associated energy W (~30 %). These are not solver errors: the 1:9 asymmetric rectangle is designed to nearly cancel the net charge, so its Qpos lands at ~0.7 nC while every other Skin metric sits above 10 nC. Any sub-nanosecond timing difference between the two solvers is amplified when the true value is that close to zero. The trajectory itself agrees; only the accounting on that particular vanishing quantity does not.
A notable side-observation: this simulator converges to a perfectly balanced
steady state (Qpos + Qneg = 0 within 10⁻⁴ nC), as
dictated by the symmetry of the biphasic rectangle. The ngspice reference
retains a residual ~5 % imbalance even after 100 periods
(is_steady_state = False in its CSV), attributable to the
accumulation of its adaptive step error.
When a single waveform is selected, the ngspice reference trajectory is drawn beneath the simulator output as a thin grey line, so the agreement can be verified visually as well as numerically. The reference data are pre-computed offline (ngspice 46, adaptive step) and embedded into the page at ~200 KB total via a per-curve adaptive mesh (RDP with vertical-distance metric, N ≈ 150-300 points per curve, worst-case reinterpolation error < 0.2 %).
Numerical scheme details
Topology. Vin → Rsrc (50 Ω) → node a → tissue → GND.
The tissue is a Foster II network: Rdc in parallel with Cinf in parallel
with N series Rk-Ck branches (N between 18 for CSF and 45 for heart muscle,
time constants τk = Rk·Ck spanning roughly 10 orders of magnitude,
from ~1 ps to ~10 s). The state variables are the N capacitor voltages Vck and
the terminal voltage Va.
Per-branch exponential integrator. Each branch obeys the linear scalar ODE dVck/dt = (Va(t) − Vck) / τk. Assuming Va(t) varies linearly across the step (first-order hold), the exact update on the interval [t, t+dt] is
Vckn+1 = φk · Vckn + αk · Van + βk · Van+1
with φk = exp(−dt/τk),
γk = (1 − φk)·τk/dt,
αk = γk − φk,
βk = 1 − γk.
Van+1 is then resolved from KCL at the terminal by a single
algebraic division. The three coefficients per branch depend only on the tissue and
the step size, so they are precomputed once. For very small dt/τ ratios the coefficients
are computed by a truncated Taylor series to sidestep catastrophic cancellation
in αk = γk − φk. The scheme is exact for each
individual branch given Va linear across the step, and therefore has no
stiffness restriction: the fastest τk ~ 1 ps and the slowest τk ~ 10 s
are handled with identical accuracy.
Modelled source rise time. The canonical rectangular waveforms have transitions that a fixed-step integrator cannot resolve, and a real hardware pulse generator does not produce instantaneous edges either. This simulator models the source as an ngspice-style PULSE with TR = 200 ns (matching the default when TR is unspecified in ngspice, namely TR = TSTEP), so each transition in Vin is replaced by a linear ramp of that duration centred on the transition instant. This one modelling choice reconciles the simulator with the ngspice reference; without it the current peaks are either diluted (fixed 200 ns step) or artificially amplified (arbitrary sub-nanosecond stepping on an idealised discontinuous step).
Adaptive time step. Away from any transition the solver uses the requested step (Fast: 200 ns · Fine: 50 ns). Inside a 5 µs window after each detected transition (identified by a pre-scan of Vin for jumps > 30 mV) the step is refined to 5 ns, which resolves both the source ramp and the fast Debye branches while the transient decays. Output samples are emitted at the fine resolution inside that window so the plots and the peak-current metrics both capture the pic. The two-scale scheme adds ~5 % to runtime.
Sign convention. Positive I = physical current entering the tissue. ngspice reports I(Vin) with the internal source convention (negative when the source drives current); this simulator flips it for comparison.
No external dependency. No JavaScript library, no network request, no local storage. The page works offline once loaded.
Physical tissue model (Foster / Cole-Cole)
Every tissue is represented by its complex admittance:
Y(ω) = 1/Rdc + jωCinf + Σn,k jωCk,n / (1 + jωRk,nCk,n)
The sum runs over the Cole-Cole dispersions (n = 1 to 4 depending on tissue) and over the 12 Debye sub-branches per dispersion that approximate them, obtained by least-squares fitting (trust-region-reflective) on a log-frequency grid from 10⁻¹ to 10¹¹ Hz. Target: relative error < 2 % on log|Y(f)| in the band 10 Hz to 10⁸ Hz.
For skin, 19 active branches are sufficient (7 for dispersion 1 at τ ≈ 7 ps, 12 for dispersion 2 at τ ≈ 32 ns).