Resources Numerical accuracy Calculation policy
Operating point, DC sweep, AC response and transient of a circuit of resistors, capacitors, inductors, sources, diodes, transistors and switches.
The equations are assembled and solved in Float64Array with dense LU and partial pivoting. Nonlinear circuits are iterated until every unknown changes by less than a relative 1e-3 plus a small absolute allowance, so the stopping test, not the arithmetic, sets the digits. The transient uses the step the person gives, with the trapezoidal rule and no error estimate. The thermal voltage is a fixed constant for 300 K.
Open item. Compute the thermal voltage from k, q and the temp option instead of the fixed 0.025852 V.
18 calculations in 2 source files. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
Tests bin/test-engines.mjs, tools/em-test/test.mjs
kT/q, the voltage scale of every diode and transistor junction.
VT = 0.025852 V
A constant with 5 significant digits. It equals k·T/q at 300.00 K to 8 digits; the comment beside it and the temp option both say 300.15 K, where the value is 0.025865 V, 0.05 % higher. The temp option is stored but not used.
const VT = 0.025852; // kT/q at 300.15 K4.7k, 100n, 1meg: a component value read from text.
value = number × factor of the suffix (t, g, meg, k, m, u, n, p, f); mil = 25.4e-6
parseFloat of the digits, then one multiplication by the suffix factor. Factors such as 1e-9 are not exact in binary, so the product can differ from the decimal value in the last bit: 100n gives 1.0000000000000001e-7. Letters after the suffix are ignored.
/** Engineering-suffix number: "4.7k", "100n", "1meg", "2u". Returns NaN when it is not a number. */
export function parseValue(s) {
const m = /^\s*([+-]?(?:\d+\.?\d*|\.\d+)(?:e[+-]?\d+)?)\s*(t|g|meg|k|mil|m|u|µ|n|p|f)?[a-z]*\s*$/i.exec(String(s));
if (!m) return NaN;
const sc = { t: 1e12, g: 1e9, meg: 1e6, k: 1e3, mil: 25.4e-6, m: 1e-3, u: 1e-6, 'µ': 1e-6, n: 1e-9, p: 1e-12, f: 1e-15 };
return parseFloat(m[1]) * (m[2] ? sc[m[2].toLowerCase()] : 1);
}The value of a DC, sine or pulse source at a time.
sin: off + amp × e^(−damp (t − delay)) × sin(2π f (t − delay) + phase); pulse: linear edges over tr and tf
Direct 64-bit evaluation. The pulse position in its period is (t − td) mod period, which loses digits as t grows against the period. A rise or fall time of zero is replaced by 1e-12 s.
/** Value of a source's waveform at time t (t = null: the DC value used for the operating point). */
export function waveAt(w, t) {
if (!w) return 0;
if (t === null || t === undefined) return w.sin ? w.sin.off : w.pulse ? w.pulse.v1 : (w.dc || 0);
if (w.sin) { const s = w.sin; if (t < (s.delay || 0)) return s.off; return s.off + s.amp * Math.exp(-(s.damp || 0) * (t - (s.delay || 0))) * Math.sin(2 * Math.PI * s.freq * (t - (s.delay || 0)) + (s.phase || 0) * Math.PI / 180); }
if (w.pulse) {
const p = w.pulse, td = p.td || 0; if (t < td) return p.v1;
const per = p.per || Infinity, tt = (t - td) % per, tr = p.tr || 1e-12, tf = p.tf || 1e-12, pw = p.pw ?? per;
if (tt < tr) return p.v1 + (p.v2 - p.v1) * tt / tr;
if (tt < tr + pw) return p.v2;
if (tt < tr + pw + tf) return p.v2 + (p.v1 - p.v2) * (tt - tr - pw) / tf;
return p.v1;
}
return w.dc || 0;
}The current through a diode and its slope at a junction voltage.
i = Is (e^(v ÷ (n VT)) − 1); g = Is e^(v ÷ (n VT)) ÷ (n VT)
The exponent is capped at 80. The current is formed as exp(x) − 1, not with expm1, so for |v| far below n·VT it loses digits in proportion to 1 ÷ |x|. The junction voltage is limited between iterations, and the step is flagged while the limit moves it by more than 1e-9 V.
if (e.kind === 'D') {
const Is = md.Is ?? 1e-14, nn = md.n ?? 1, vt = nn * VT, vcrit = vt * Math.log(vt / (Math.SQRT2 * Is));
const vraw = vf(t[0]) - vf(t[1]); let vd = vraw; if (st.vj !== undefined) vd = pnjlim(vraw, st.vj, vt, vcrit); st.vj = vd; st.lim = Math.abs(vd - vraw) > 1e-9;
const ex = Math.exp(Math.min(vd / vt, 80)), id = Is * (ex - 1), gd = Is * ex / vt;
return { controls: [{ p: t[0], m: t[1], v0: vd }], rows: [{ node: t[0], i0: id, d: [gd] }, { node: t[1], i0: -id, d: [-gd] }] };
}Collector and base current of a bipolar transistor.
ic = If − Ir − Ir ÷ βr; ib = If ÷ βf + Ir ÷ βr; If = Is (e^(vbe ÷ VT) − 1)
The transport model with the same capped exponent and exp(x) − 1 form as the diode. The model has no Early effect, capacitances or series resistance.
const s = md.type === 'pnp' ? -1 : 1, Is = md.Is ?? 1e-16, bf = md.bf ?? 100, br = md.br ?? 1, vt = VT, vcrit = vt * Math.log(vt / (Math.SQRT2 * Is));
const [c, b, em] = t, be0 = s * (vf(b) - vf(em)), bc0 = s * (vf(b) - vf(c));
let be = be0, bc = bc0; if (st.vbe !== undefined) { be = pnjlim(be0, st.vbe, vt, vcrit); bc = pnjlim(bc0, st.vbc, vt, vcrit); }
st.vbe = be; st.vbc = bc; st.lim = Math.abs(be - be0) > 1e-9 || Math.abs(bc - bc0) > 1e-9;
const ef = Math.exp(Math.min(be / vt, 80)), er = Math.exp(Math.min(bc / vt, 80)), If = Is * (ef - 1), Ir = Is * (er - 1), gf = Is * ef / vt, gr = Is * er / vt;
const ic = If - Ir - Ir / br, ib = If / bf + Ir / br;
const dic = [gf, -gr - gr / br], dib = [gf / bf, gr / br]; // wrt (vbe, vbc) in transformed sense
// actual currents = s * transformed; actual controls v = s * transformed ⇒ derivative wrt actual (vb-ve),(vb-vc) is unchanged
return { controls: [{ p: b, m: em, v0: s * be }, { p: b, m: c, v0: s * bc }], rows: [
{ node: c, i0: s * ic, d: dic }, { node: b, i0: s * ib, d: dib }, { node: em, i0: -s * (ic + ib), d: [-(dic[0] + dib[0]), -(dic[1] + dib[1])] }] };
}Drain current of a MOSFET in the triode and saturation regions.
triode: id = kp (vov vds − vds² ÷ 2)(1 + λ vds); saturation: id = kp vov² (1 + λ vds) ÷ 2
Closed forms in 64-bit with the conductances as their derivatives. Drain and source are exchanged when vds is negative. There is no subthreshold current: id is 0 at or below threshold.
if (e.kind === 'M') {
const s = md.type === 'pmos' ? -1 : 1, kp = (md.kp ?? 2e-5) * ((md.w ?? 1) / (md.l ?? 1)), vto = md.vto ?? (s > 0 ? 1 : -1), lam = md.lambda ?? 0;
let [d, g, sr] = t; const vds0 = s * (vf(d) - vf(sr)); if (vds0 < 0) [d, sr] = [sr, d];
const vgs = s * (vf(g) - vf(sr)), vds = s * (vf(d) - vf(sr)), vov = vgs - s * vto;
let id = 0, gm = 0, gds = 0;
if (vov > 0) {
if (vds < vov) { id = kp * (vov * vds - vds * vds / 2) * (1 + lam * vds); gm = kp * vds * (1 + lam * vds); gds = kp * (vov - vds) * (1 + lam * vds) + kp * (vov * vds - vds * vds / 2) * lam; }
else { id = 0.5 * kp * vov * vov * (1 + lam * vds); gm = kp * vov * (1 + lam * vds); gds = 0.5 * kp * vov * vov * lam; }
}
return { controls: [{ p: g, m: sr, v0: s * vgs }, { p: d, m: sr, v0: s * vds }], rows: [{ node: d, i0: s * id, d: [gm, gds] }, { node: sr, i0: -s * id, d: [-gm, -gds] }] };
}The conductance of a voltage-controlled switch.
ln g = ln(1 ÷ roff) + a (ln(1 ÷ ron) − ln(1 ÷ roff)); a = 1 ÷ (1 + e^(−(v − vt) ÷ vh))
The conductance moves smoothly between off and on on a log scale, so it approaches 1 ÷ ron and 1 ÷ roff without being set to them. The transition width has a floor of 1e-6 V.
_swG(e, vf) { const vt = e.p.vt ?? 0.5, vh = e.p.vh ?? 0.05, ron = e.p.ron ?? 1, roff = e.p.roff ?? 1e9, v = vf(e.nodes[2]) - vf(e.nodes[3]); const a = 1 / (1 + Math.exp(-(v - vt) / Math.max(vh, 1e-6))); return Math.exp(Math.log(1 / roff) + a * (Math.log(1 / ron) - Math.log(1 / roff))); }The matrix and right side for one Newton step, at DC or at one time step.
capacitor: geq = C ÷ h (backward Euler) or 2C ÷ h (trapezoid); inductor: r = L ÷ h or 2L ÷ h
Stamps are added into Float64Array rows in element order. gmin (1e-12 S by default) is added from every node to ground, which shifts the answer of a circuit whose node resistances approach 1e12 Ω. Each nonlinear element is linearized about the current guess.
/** Assemble A x = b for 'dc' | 'tran' (real) around the guess x. ctx: {t, h, trap, src (scale), gmin, state}. */
_assemble(ix, x, ctx) {
const n = ix.size, A = Array.from({ length: n }, () => new Float64Array(n)), b = new Float64Array(n), { map } = ix;
const ni = (nd) => (isGnd(nd) ? -1 : map.get(nd) - 1);
const vf = (nd) => { const i = ni(nd); return i < 0 ? 0 : x[i]; };
const g2 = (a, c, g) => { if (a >= 0) A[a][a] += g; if (c >= 0) A[c][c] += g; if (a >= 0 && c >= 0) { A[a][c] -= g; A[c][a] -= g; } };
const rhs = (a, v) => { if (a >= 0) b[a] += v; };
for (let i = 0; i < ix.nodeCount; i++) A[i][i] += ctx.gmin;
for (const e of this.elements) {
const a = ni(e.nodes[0]), c = ni(e.nodes[1]), st = ctx.state.get(e.name);
switch (e.kind) {
case 'R': g2(a, c, 1 / e.p.value); break;
case 'S': g2(a, c, this._swG(e, vf)); break;
case 'C':
if (ctx.h) { const geq = (ctx.trap ? 2 : 1) * e.p.value / ctx.h, ieq = geq * st.v + (ctx.trap ? st.i : 0); g2(a, c, geq); rhs(a, ieq); rhs(c, -ieq); }
break;
case 'V': case 'E': case 'L': {
const k = ix.br.get(e.name); if (a >= 0) { A[a][k] += 1; A[k][a] += 1; } if (c >= 0) { A[c][k] -= 1; A[k][c] -= 1; }
if (e.kind === 'V') b[k] += waveAt(e.p.wave, ctx.t) * ctx.src;
else if (e.kind === 'E') { const cp = ni(e.nodes[2]), cn = ni(e.nodes[3]); if (cp >= 0) A[k][cp] -= e.p.gain; if (cn >= 0) A[k][cn] += e.p.gain; }
else if (ctx.h) { const r = (ctx.trap ? 2 : 1) * e.p.value / ctx.h; A[k][k] -= r; b[k] -= r * st.i + (ctx.trap ? st.v : 0); }
break;
}
case 'I': { const v = waveAt(e.p.wave, ctx.t) * ctx.src; rhs(a, -v); rhs(c, v); break; }
case 'G': { const cp = ni(e.nodes[2]), cn = ni(e.nodes[3]), g = e.p.gm; for (const [r, sg] of [[a, 1], [c, -1]]) if (r >= 0) { if (cp >= 0) A[r][cp] += sg * g; if (cn >= 0) A[r][cn] -= sg * g; } break; }
default: {
const nl = this._nl(e, vf, st); if (!nl) break;
for (const r of nl.rows) { const row = ni(r.node); if (row < 0) continue; let sum = 0;
r.d.forEach((d, k) => { const p = ni(nl.controls[k].p), m = ni(nl.controls[k].m); if (p >= 0) A[row][p] += d; if (m >= 0) A[row][m] -= d; sum += d * nl.controls[k].v0; });
b[row] -= r.i0 - sum; }
}
}
}
return { A, b };
}Repeats the linear solve until the voltages and currents stop changing.
stop when |Δxᵢ| ≤ 1e-3 × max(|xᵢ new|, |xᵢ old|) + 1e-6 V (or 1e-12 A) for every unknown
The test is on the change between two iterations, with a relative part of 1e-3. It needs at least two iterations and no junction limiting in the last one. The digits of an operating point are set by this test and by how fast the iteration was converging, not by 64-bit rounding. At most 100 iterations.
_newton(ix, x0, ctx) {
let x = Float64Array.from(x0); const o = this.opts;
for (let it = 1; it <= o.maxIter; it++) {
for (const st of ctx.state.values()) st.lim = false;
const { A, b } = this._assemble(ix, x, ctx), xn = solveReal(A, b);
if (!xn) return { ok: false, singular: true, x, iterations: it };
let conv = true; for (let i = 0; i < xn.length; i++) { if (!Number.isFinite(xn[i])) return { ok: false, x, iterations: it, nan: true }; const tol = o.reltol * Math.max(Math.abs(xn[i]), Math.abs(x[i])) + (i < ix.nodeCount ? o.vntol : o.abstol); if (Math.abs(xn[i] - x[i]) > tol) { conv = false; break; } }
const limited = [...ctx.state.values()].some((s) => s.lim);
x = xn; if (conv && !limited && it > 1) return { ok: true, x, iterations: it };
}
return { ok: false, x, iterations: o.maxIter };
}The DC voltages and currents of the circuit.
direct Newton, then gmin stepping, then source stepping
Tries Newton from zero, then gmin stepping from 1e-2 S down to the set gmin in six steps, then source stepping in 20 equal steps. The method that worked is returned with the result.
/** DC operating point. Tries direct Newton, then gmin stepping, then source stepping. */
op(opts = {}) {
const pre = this.check().filter((p) => p.severity === 'error' && p.id !== 'em:circuit.floating-node');
if (pre.length) return { ok: false, problems: pre };
const ix = this._index(), state = this._freshState(), n = ix.size, o = this.opts;
const base = { t: null, h: 0, trap: false, src: 1, gmin: o.gmin, state };
let r = this._newton(ix, new Float64Array(n), base), how = 'direct';
if (!r.ok) { // gmin stepping
let x = new Float64Array(n), okAll = true; for (const g of [1e-2, 1e-4, 1e-6, 1e-8, 1e-10, o.gmin]) { const rr = this._newton(ix, x, { ...base, gmin: g, state: this._freshState() }); if (!rr.ok) { okAll = false; break; } x = rr.x; }
if (okAll) { r = this._newton(ix, x, { ...base, state }); how = 'gmin stepping'; }
}
if (!r.ok) { // source stepping
let x = new Float64Array(n), okAll = true; for (let k = 1; k <= 20; k++) { const rr = this._newton(ix, x, { ...base, src: k / 20, state: this._freshState() }); if (!rr.ok) { okAll = false; break; } x = rr.x; }
if (okAll) { r = this._newton(ix, x, { ...base, state }); how = 'source stepping'; }
}
if (!r.ok) {
const problems = this.check();
problems.push(problem(r.singular ? 'circuit.singular' : 'circuit.no-convergence', r.singular ? 'The circuit equations are singular (a node or loop is not determined).' : 'The operating point did not converge after direct, gmin and source stepping.', { fix: r.singular ? 'See the floating-node and loop problems.' : 'Check diode and transistor polarity, add series resistance, or give a start state.' }));
return { ok: false, problems };
}
this._x0 = r.x; this._ix = ix;
return { ok: true, ...this._pack(ix, r.x), iterations: r.iterations, method: how, problems: this.check().filter((p) => p.severity !== 'error') };
}The circuit solved for each value of a swept source.
Each point starts from the previous solution. If that fails, the point is solved again from zero with gmin raised to 1e-9 S, and the result of that point then includes the larger leak without a message.
/** DC sweep of a V or I source: values array. Returns {sweep, v: {node:[…]}, i:{name:[…]}}. Continuation from the previous point. */
dc(srcName, values) {
const e = this.elements.find((x) => x.name === srcName && (x.kind === 'V' || x.kind === 'I'));
if (!e) return { ok: false, problems: [problem('circuit.no-source', 'No source named ' + srcName + '.')] };
const keep = e.p.wave, ix = this._index(), state = this._freshState(), res = { ok: true, sweep: [...values], v: {}, i: {}, problems: [] };
let x = new Float64Array(ix.size);
for (const val of values) {
e.p.wave = { dc: val }; let r = this._newton(ix, x, { t: null, h: 0, trap: false, src: 1, gmin: this.opts.gmin, state });
if (!r.ok) r = this._newton(ix, new Float64Array(ix.size), { t: null, h: 0, trap: false, src: 1, gmin: 1e-9, state: this._freshState() });
if (!r.ok) { res.ok = false; res.problems.push(problem('circuit.no-convergence', 'The sweep stopped at ' + srcName + ' = ' + val + ': no convergence.')); break; }
x = r.x; const p = this._pack(ix, x); for (const k in p.v) (res.v[k] ||= []).push(p.v[k]); for (const k in p.i) (res.i[k] ||= []).push(p.i[k]);
}
e.p.wave = keep; return res;
}Complex node voltages and branch currents at each frequency.
capacitor: jωC; inductor: −jωL on its branch row; ω = 2π f
Real and imaginary parts are kept in separate Float64Array matrices and solved by complex LU at each frequency. Nonlinear elements enter with their conductances at the operating point. Source phase is converted from degrees with π ÷ 180.
/** AC small-signal at the operating point. freqs in Hz. Sources use wave.ac = {mag, phase°}. Returns {f, v: {node: {re, im}[]}, i}. */
ac(freqs) {
const o = this.op(); if (!o.ok) return o;
const ix = this._index(), x = this._x0, n = ix.size, { map } = ix, state = this._freshState();
const ni = (nd) => (isGnd(nd) ? -1 : map.get(nd) - 1), vf = (nd) => { const i = ni(nd); return i < 0 ? 0 : x[i]; };
for (const e of this.elements) if ('DQM'.includes(e.kind)) { const s = state.get(e.name); this._nl(e, vf, s); } // sets limiting reference to the op point
const res = { ok: true, f: [...freqs], v: {}, i: {}, problems: [] };
for (const f of freqs) {
const w = 2 * Math.PI * f, Ar = Array.from({ length: n }, () => new Float64Array(n)), Ai = Array.from({ length: n }, () => new Float64Array(n)), br = new Float64Array(n), bi = new Float64Array(n);
const add = (M, r, c, v) => { if (r >= 0 && c >= 0) M[r][c] += v; };
const g2 = (M, a, c, g) => { add(M, a, a, g); add(M, c, c, g); add(M, a, c, -g); add(M, c, a, -g); };
for (let i = 0; i < ix.nodeCount; i++) Ar[i][i] += this.opts.gmin;
for (const e of this.elements) {
const a = ni(e.nodes[0]), c = ni(e.nodes[1]);
switch (e.kind) {
case 'R': g2(Ar, a, c, 1 / e.p.value); break;
case 'S': g2(Ar, a, c, this._swG(e, vf)); break;
case 'C': g2(Ai, a, c, w * e.p.value); break;
case 'V': case 'E': case 'L': {
const k = ix.br.get(e.name); add(Ar, a, k, 1); add(Ar, k, a, 1); add(Ar, c, k, -1); add(Ar, k, c, -1);
if (e.kind === 'V') { const ac = e.p.wave?.ac; if (ac) { const ph = (ac.phase || 0) * Math.PI / 180; br[k] += ac.mag * Math.cos(ph); bi[k] += ac.mag * Math.sin(ph); } }
else if (e.kind === 'E') { add(Ar, k, ni(e.nodes[2]), -e.p.gain); add(Ar, k, ni(e.nodes[3]), e.p.gain); }
else Ai[k][k] -= w * e.p.value;
break;
}
case 'I': { const ac = e.p.wave?.ac; if (ac) { const ph = (ac.phase || 0) * Math.PI / 180, vr = ac.mag * Math.cos(ph), vi = ac.mag * Math.sin(ph); if (a >= 0) { br[a] -= vr; bi[a] -= vi; } if (c >= 0) { br[c] += vr; bi[c] += vi; } } break; }
case 'G': { const cp = ni(e.nodes[2]), cn = ni(e.nodes[3]), g = e.p.gm; for (const [r, sg] of [[a, 1], [c, -1]]) { add(Ar, r, cp, sg * g); add(Ar, r, cn, -sg * g); } break; }
default: { const nl = this._nl(e, vf, state.get(e.name)); if (!nl) break; for (const r of nl.rows) { const row = ni(r.node); r.d.forEach((d, k) => { add(Ar, row, ni(nl.controls[k].p), d); add(Ar, row, ni(nl.controls[k].m), -d); }); } }
}
}
const sol = solveComplex(Ar, Ai, br, bi);
if (!sol) { res.ok = false; res.problems.push(problem('circuit.singular', 'The AC matrix is singular at ' + f + ' Hz.')); break; }
for (const [nd, i] of map) (res.v[nd] ||= []).push({ re: sol.re[i - 1], im: sol.im[i - 1] });
for (const [nm, k] of ix.br) (res.i[nm] ||= []).push({ re: sol.re[k], im: sol.im[k] });
}
return res;
}Voltages and currents over time.
trapezoidal rule after a first backward-Euler step; output at t = k × tstep
Output times are k × tstep, not a running sum. The step is the print step; it is only halved (up to 12 times) when Newton fails, never because of an error estimate, so how close a waveform is depends on the step the person chose.
/** Transient. {tstop, tstep, uic:false, method:'trap'|'be'}. Returns {t:[…], v:{node:[…]}, i:{…}}. */
tran({ tstop, tstep, uic = false, method = 'trap' }) {
const ix = this._index(), state = this._freshState(), n = ix.size, o = this.opts, res = { ok: true, t: [], v: {}, i: {}, problems: [] };
let x = new Float64Array(n);
if (!uic) { const r = this.op(); if (!r.ok) return r; x = Float64Array.from(this._x0); }
const pull = () => { for (const e of this.elements) { const st = state.get(e.name), a = e.nodes[0], c = e.nodes[1], vf = (nd) => (isGnd(nd) ? 0 : x[ix.map.get(nd) - 1]);
if (e.kind === 'C') { st.v = uic && e.p.ic !== undefined && res.t.length === 0 ? e.p.ic : vf(a) - vf(c); st.i = 0; }
if (e.kind === 'L') { const k = ix.br.get(e.name); st.i = uic && e.p.ic !== undefined && res.t.length === 0 ? e.p.ic : x[k]; st.v = 0; } } };
pull();
const record = (t) => { res.t.push(t); const p = this._pack(ix, x); for (const k in p.v) (res.v[k] ||= []).push(p.v[k]); for (const k in p.i) (res.i[k] ||= []).push(p.i[k]); };
record(0);
// dynamic-element history: after a solved step, update v,i per the companion formulas
const commit = (h, trap) => { for (const e of this.elements) { const st = state.get(e.name), vf = (nd) => (isGnd(nd) ? 0 : x[ix.map.get(nd) - 1]);
if (e.kind === 'C') { const vn = vf(e.nodes[0]) - vf(e.nodes[1]); const geq = (trap ? 2 : 1) * e.p.value / h; st.i = geq * (vn - st.v) - (trap ? st.i : 0); st.v = vn; }
if (e.kind === 'L') { const k = ix.br.get(e.name), vn = vf(e.nodes[0]) - vf(e.nodes[1]); st.v = vn; st.i = x[k]; } } };
let t = 0, first = true, k = 0;
while (k * tstep < tstop - tstep * 1e-9) {
k++; const target = Math.min(k * tstep, tstop);
let guard = 0;
while (target - t > tstep * 1e-9) {
let h = target - t, ok = false;
for (let tries = 0; tries < 12 && !ok; tries++) {
const trap = method === 'trap' && !first, snapshot = new Map([...state].map(([kk, v]) => [kk, { ...v }]));
const r = this._newton(ix, x, { t: t + h, h, trap, src: 1, gmin: o.gmin, state });
if (r.ok) { x = r.x; commit(h, trap); t += h; ok = true; first = false; }
else { for (const [kk, v] of snapshot) state.set(kk, v); h /= 2; }
}
if (!ok || ++guard > 4096) { res.ok = false; res.problems.push(problem('circuit.no-convergence', 'The transient stopped at t = ' + t + ' s: no convergence even with the step cut a thousandfold.', { fix: 'Smaller step, or soften the fastest edge.' })); return res; }
}
t = target; record(t);
}
return res;
}Keeps a diode or transistor voltage from jumping too far in one iteration.
v = v_old + vt ln(1 + (v_new − v_old) ÷ vt) when the step exceeds 2 vt above vcrit
It changes the path of the iteration, not the converged answer, because convergence is refused while a limit is moving a voltage.
function pnjlim(vnew, vold, vt, vcrit) {
if (vnew > vcrit && Math.abs(vnew - vold) > 2 * vt) {
if (vold > 0) { const arg = 1 + (vnew - vold) / vt; return arg > 0 ? vold + vt * Math.log(arg) : vcrit; }
return vt * Math.log(vnew / vt);
}
return vnew;
}The source values of a DC sweep and the frequencies of an AC analysis.
dc: start, start + step, …; ac dec or oct: f = fstart × (fstop ÷ fstart)^(k ÷ N); lin: fstart + (fstop − fstart) k ÷ (points − 1)
DC sweep values are built by repeated addition of the step, so later values carry accumulated rounding; the end test allows 1e-9 of a step. AC frequencies are computed from the index, with no accumulation.
/** Run a deck's analyses in order; returns [{type, result}]. */
export function runDeck(deck) {
const out = [];
for (const a of deck.analyses) {
if (a.type === 'op') out.push({ type: 'op', result: deck.circuit.op() });
else if (a.type === 'dc') { const v = []; const dir = a.step > 0 ? 1 : -1; for (let x = a.start; dir * (x - a.stop) <= Math.abs(a.step) * 1e-9; x += a.step) v.push(x); out.push({ type: 'dc', result: deck.circuit.dc(a.source, v) }); }
else if (a.type === 'ac') { const f = []; if (a.mode === 'lin') for (let k = 0; k < a.points; k++) f.push(a.fstart + (a.fstop - a.fstart) * k / Math.max(a.points - 1, 1)); else { const per = a.mode === 'dec' ? Math.log10(a.fstop / a.fstart) : Math.log2(a.fstop / a.fstart), N = Math.max(1, Math.round(per * a.points)); for (let k = 0; k <= N; k++) f.push(a.fstart * Math.pow(a.fstop / a.fstart, k / N)); } out.push({ type: 'ac', result: deck.circuit.ac(f) }); }
else if (a.type === 'tran') out.push({ type: 'tran', result: deck.circuit.tran({ tstop: a.tstop, tstep: a.tstep, uic: a.uic }) });
}
return out;
}The size and angle of a complex result.
|z| = hypot(re, im); phase = atan2(im, re) × 180 ÷ π
Math.hypot avoids overflow of the squares. One 64-bit expression each.
export const mag = (z) => Math.hypot(z.re, z.im);
export const phaseDeg = (z) => Math.atan2(z.im, z.re) * 180 / Math.PI;Tests bin/test-engines.mjs, tools/em-test/test.mjs
Solves the circuit equations of one Newton step.
LU with partial pivoting, then back substitution
Works on Float64Array copies. The pivot is the largest entry in the column; rows are not scaled first and there is no iterative refinement. The matrix is called singular only when the pivot is not above 1e-300, an absolute value, so a nearly singular circuit is solved and returned without a warning.
/** Solve A x = b (A: n×n array of Float64Array rows or flat), in place on copies. Returns x or null when singular. */
export function solveReal(A, b, tiny = 1e-300) {
const n = b.length, a = A.map((r) => Float64Array.from(r)), x = Float64Array.from(b);
for (let k = 0; k < n; k++) {
let p = k, mx = Math.abs(a[k][k]);
for (let i = k + 1; i < n; i++) if (Math.abs(a[i][k]) > mx) { mx = Math.abs(a[i][k]); p = i; }
if (!(mx > tiny)) return null;
if (p !== k) { [a[p], a[k]] = [a[k], a[p]]; [x[p], x[k]] = [x[k], x[p]]; }
for (let i = k + 1; i < n; i++) {
const f = a[i][k] / a[k][k]; if (f === 0) continue;
for (let j = k; j < n; j++) a[i][j] -= f * a[k][j];
x[i] -= f * x[k];
}
}
for (let i = n - 1; i >= 0; i--) { let s = x[i]; for (let j = i + 1; j < n; j++) s -= a[i][j] * x[j]; x[i] = s / a[i][i]; }
return x;
}Solves the AC equations and the antenna's impedance matrix.
LU with partial pivoting on separate real and imaginary arrays
The same elimination as the real solve. Each complex division multiplies by the conjugate and divides by re² + im², which can overflow or underflow for entries near the ends of the 64-bit range.
/** Complex solve: Ar, Ai are arrays of Float64Array rows; br, bi vectors. Returns {re, im} or null when singular. */
export function solveComplex(Ar, Ai, br, bi, tiny = 1e-300) {
const n = br.length, ar = Ar.map((r) => Float64Array.from(r)), ai = Ai.map((r) => Float64Array.from(r)), xr = Float64Array.from(br), xi = Float64Array.from(bi);
for (let k = 0; k < n; k++) {
let p = k, mx = Math.hypot(ar[k][k], ai[k][k]);
for (let i = k + 1; i < n; i++) { const m = Math.hypot(ar[i][k], ai[i][k]); if (m > mx) { mx = m; p = i; } }
if (!(mx > tiny)) return null;
if (p !== k) { [ar[p], ar[k]] = [ar[k], ar[p]]; [ai[p], ai[k]] = [ai[k], ai[p]]; [xr[p], xr[k]] = [xr[k], xr[p]]; [xi[p], xi[k]] = [xi[k], xi[p]]; }
const d = ar[k][k] * ar[k][k] + ai[k][k] * ai[k][k];
for (let i = k + 1; i < n; i++) {
// f = a[i][k] / a[k][k]
const fr = (ar[i][k] * ar[k][k] + ai[i][k] * ai[k][k]) / d, fi = (ai[i][k] * ar[k][k] - ar[i][k] * ai[k][k]) / d;
if (fr === 0 && fi === 0) continue;
for (let j = k; j < n; j++) { ar[i][j] -= fr * ar[k][j] - fi * ai[k][j]; ai[i][j] -= fr * ai[k][j] + fi * ar[k][j]; }
xr[i] -= fr * xr[k] - fi * xi[k]; xi[i] -= fr * xi[k] + fi * xr[k];
}
}
for (let i = n - 1; i >= 0; i--) {
let sr = xr[i], si = xi[i];
for (let j = i + 1; j < n; j++) { sr -= ar[i][j] * xr[j] - ai[i][j] * xi[j]; si -= ar[i][j] * xi[j] + ai[i][j] * xr[j]; }
const d = ar[i][i] * ar[i][i] + ai[i][i] * ai[i][i];
xr[i] = (sr * ar[i][i] + si * ai[i][i]) / d; xi[i] = (si * ar[i][i] - sr * ai[i][i]) / d;
}
return { re: xr, im: xi };
}The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.