Resources Numerical accuracy Calculation policy
Input impedance, segment currents, radiation pattern and directivity of thin-wire antennas.
The impedance matrix is filled and solved in Float64Array with complex LU. How close the result is comes from the model: the number of segments, the thin-wire kernel, Gauss–Legendre rules of 3 or 4 points, and a midpoint sum over the sphere for the radiated power. The file's own header says its numbers are not identical to those of the reference code whose equation it solves.
8 calculations in 1 source file. 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
The speed of light, μ0, ε0 and the impedance of free space.
c = 299792458 m/s; μ0 = 4π × 10⁻⁷; ε0 = 1 ÷ (μ0 c²); η0 = μ0 c
c is the exact defined value. μ0 is the value fixed before 2019, and ε0 and η0 are derived from it, so ε0 here (8.854187817620e-12) differs from the field solver's 8.8541878128e-12 in the 10th digit.
export const C0 = 299792458, MU0 = 4e-7 * Math.PI, EPS0 = 1 / (MU0 * C0 * C0), ETA0 = MU0 * C0;The sample points and weights of the integrals along segments.
2, 3 and 4 points on [0, 1]
Nodes and weights are decimal constants with 15 significant digits; the 3-point weights are the fractions 5/18 and 4/9.
// Gauss–Legendre nodes/weights on [0,1]
const GL = (n) => { const t = { 2: [[0.5 - 0.288675134594813, 0.5], [0.5 + 0.288675134594813, 0.5]], 3: [[0.5 - 0.387298334620742, 5 / 18], [0.5, 4 / 9], [0.5 + 0.387298334620742, 5 / 18]],
4: [[0.0694318442029737, 0.173927422568727], [0.330009478207572, 0.326072577431273], [0.669990521792428, 0.326072577431273], [0.930568155797026, 0.173927422568727]] }; return t[n]; };How wires are cut into segments and joined at shared points.
segments per leg = ⌈L ÷ maxSeg⌉, maxSeg = λ ÷ 20 by default; merge distance = 1e-6 × max(largest coordinate, λ ÷ 1000)
Interior nodes are placed at a + (b − a) s ÷ n. Points closer than the merge distance become one node, which is how wires are joined. When a segment count is given it is divided evenly over the legs of the polyline and rounded.
_build(wires) {
let big = 0; for (const w of wires) for (const p of w.pts) big = Math.max(big, Math.abs(p[0]), Math.abs(p[1]), Math.abs(p[2])); this.tol = Math.max(big, this.lambda * 1e-3) * 1e-6;
for (const w of wires) {
const maxSeg = w.maxSeg || this.lambda / 20;
if (!(w.radius > 0)) this.problems.push(P('bad-radius', 'A wire has no radius; the thin-wire kernel needs one.'));
if (w.radius > this.lambda / 50) this.problems.push(P('thick-wire', 'A wire of radius ' + w.radius + ' m is thicker than λ/50; the thin-wire model is poor.', { severity: 'warning' }));
for (let i = 0; i + 1 < w.pts.length; i++) {
const a = w.pts[i], b = w.pts[i + 1], L = len(sub(b, a)); if (L === 0) { this.problems.push(P('zero-leg', 'A wire has a zero-length leg.')); continue; }
const n = w.segments ? Math.max(1, Math.round(w.segments / Math.max(1, w.pts.length - 1))) : Math.max(1, Math.ceil(L / maxSeg)), ids = [this._node(a)];
for (let s = 1; s < n; s++) ids.push(this._node(add(a, mul(sub(b, a), s / n))));
ids.push(this._node(b));
for (let s = 0; s < n; s++) { const pa = this.nodes[ids[s]], pb = this.nodes[ids[s + 1]], d = sub(pb, pa), l = len(d); this.segs.push({ a: ids[s], b: ids[s + 1], len: l, u: mul(d, 1 / l), radius: w.radius, mid: mul(add(pa, pb), 0.5) }); }The coupling between two segments.
∫∫ cos(kR) ÷ (4πR), R = √(|r − r′|² + a²)
The reduced kernel adds the wire radius squared under the root, so R never reaches zero. Near pairs use a 3-point rule on 4 sub-intervals (8 for a segment with itself); far pairs and images one 4-point rule. The self term is sharply peaked for thin wires, so its value depends on the sub-division, not on 64-bit rounding.
// ∫∫ G over segment pair, with end-ramp weights: returns {g, r: [[r00, r01],[r10, r11]]} where r[e][f] = ∫∫ ramp_e(s) ramp_f(s') G
_pair(i, j, image) {
const s = this.segs[i], t = this.segs[j], close = len(sub(s.mid, t.mid)) < 2.5 * Math.max(s.len, t.len) && !image;
const nsub = close ? (i === j ? 8 : 4) : 1, Q = close ? Q3 : Q4, a2 = t.radius * t.radius; let g = 0; const r = [[0, 0], [0, 0]];
const pa = this.nodes[s.a], tb = this.nodes[t.a];
for (let m = 0; m < nsub; m++) for (const [xs, ws] of Q) {
const fs = (m + xs) / nsub, ps = add(pa, mul(s.u, fs * s.len)), wsm = ws * s.len / nsub;
for (let n2 = 0; n2 < nsub; n2++) for (const [xt, wt] of Q) {
const ft = (n2 + xt) / nsub; let pt = add(tb, mul(t.u, ft * t.len)); if (image) pt = [pt[0], pt[1], -pt[2]];
const R = Math.sqrt(dot(sub(ps, pt), sub(ps, pt)) + a2), kr = this.k * R, w = wsm * wt * t.len / nsub, gr = Math.cos(kr) / (4 * Math.PI * R);
g += w * gr; r[0][0] += w * gr * (1 - fs) * (1 - ft); r[0][1] += w * gr * (1 - fs) * ft; r[1][0] += w * gr * fs * (1 - ft); r[1][1] += w * gr * fs * ft;
}
}
return { g, r, ...this._pairIm(i, j, image, nsub, Q) };
}The matrix that links the currents to the source voltage.
Z_mn = jωμ (u · u′) ∫∫ f_m f_n G + (1 ÷ jωε) ∫∫ f_m′ f_n′ G
Sums over the two halves of each basis function and, with a ground, their images. The code's own warning gives about ±25 % for the radiation resistance of an antenna smaller than λ ÷ 20, where it is a small remainder of large reactive terms.
/** Fill the impedance matrix. Called by solve(). */
_matrix() {
const nb = this.basis.length, Zr = Array.from({ length: nb }, () => new Float64Array(nb)), Zi = Array.from({ length: nb }, () => new Float64Array(nb)), cache = new Map();
const wm = this.omega * MU0, wi = 1 / (this.omega * EPS0); // jωμ ; 1/(jωε) = -j/(ωε)
const get = (i, j, image) => { const key = (image ? 'i' : 'r') + i + ',' + j; let v = cache.get(key); if (!v) { v = this._pair(i, j, image); cache.set(key, v); } return v; };
for (let m = 0; m < nb; m++) for (let n = 0; n < nb; n++) {
let zr = 0, zi = 0;
for (const pm of this.basis[m].parts) for (const pn of this.basis[n].parts) {
const s = this.segs[pm.seg], t = this.segs[pn.seg], dm = pm.sigma * (pm.end === 1 ? 1 : -1) / s.len, dn = pn.sigma * (pn.end === 1 ? 1 : -1) / t.len;
for (const image of this.ground === 'perfect' ? [false, true] : [false]) {
const v = get(pm.seg, pn.seg, image), ud = image ? [-t.u[0], -t.u[1], t.u[2]] : t.u, cs = dot(s.u, ud), dsign = image ? -1 : 1;
const rr = v.r[pm.end][pn.end], ri = v.ri[pm.end][pn.end];
// jωμ σσ (u·u') ⟨r r G⟩ → imaginary·(rr + j ri) = −ωμ ri·… : j·(rr + j ri) = −ri + j rr
zr += -wm * pm.sigma * pn.sigma * cs * ri; zi += wm * pm.sigma * pn.sigma * cs * rr;
// (1/(jωε)) d d' ⟨G⟩ : −j/(ωε)·(g + j gi) = gi/(ωε)·… = (gi + (−j)(g)) → real = gi·wi, imag = −g·wi
zr += wi * dm * dn * dsign * v.gi; zi += -wi * dm * dn * dsign * v.g;
}
}
Zr[m][n] = zr; Zi[m][n] = zi;
}
for (const L of this.loads) { const nd = this.nodes.findIndex((p) => len(sub(p, L.at)) < this.tol * 10); const bi = this.basis.findIndex((b) => b.node === nd); if (bi >= 0) { Zr[bi][bi] += L.z.re; Zi[bi][bi] += L.z.im; } else this.problems.push(P('load-nowhere', 'A load is not at a junction or ground point of the wires.', { severity: 'warning' })); }
return { Zr, Zi };
}The impedance seen at the feed point.
Z_in = V ÷ I at the feed
One complex LU solve, then a complex division written out with re² + im² in the denominator. The feed is the node nearest the given point, with no check on how far away that node is.
/** Solve with a 1 V delta-gap source at the node nearest `feed` (a point). Returns {Zin, I (per segment mid, complex), feedCurrent, problems}. */
solve({ feed, voltage = { re: 1, im: 0 } }) {
if (this.problems.some((p) => p.severity === 'error')) return { ok: false, problems: this.problems };
const nd = this.nodes.reduce((best, p, i) => (len(sub(p, feed)) < len(sub(this.nodes[best], feed)) ? i : best), 0);
const fb = this.basis.findIndex((b) => b.node === nd);
if (fb < 0) return { ok: false, problems: [...this.problems, P('feed-nowhere', 'The feed point is not at a junction, a wire-to-wire joint or the ground; a delta gap needs a node with current.', { fix: 'Place the feed at an interior point of a wire.' })] };
const { Zr, Zi } = this._matrix(), nb = this.basis.length, br = new Float64Array(nb), bi = new Float64Array(nb); br[fb] = voltage.re; bi[fb] = voltage.im;
const sol = solveComplex(Zr, Zi, br, bi); if (!sol) return { ok: false, problems: [...this.problems, P('singular', 'The impedance matrix is singular.')] };
this.I = { re: sol.re, im: sol.im }; this.feedBasis = fb;
const d = sol.re[fb] ** 2 + sol.im[fb] ** 2, Zin = { re: (voltage.re * sol.re[fb] + voltage.im * sol.im[fb]) / d, im: (voltage.im * sol.re[fb] - voltage.re * sol.im[fb]) / d };
this.V = voltage; this.Zin = Zin;
return { ok: true, Zin, feedCurrent: { re: sol.re[fb], im: sol.im[fb] }, problems: this.problems, currents: this.segmentCurrents() };
}Power per solid angle in a direction.
U = η0 k² ÷ (32π²) × (|N_θ|² + |N_φ|²)
The radiation vector is integrated along each segment with a 4-point rule. 64-bit sums of products of sines and cosines.
/** Radiation intensity U(θ,φ) in W/sr for the solved current (1 V source), plus the θ and φ components. */
intensity(th, ph) {
const { Nr, Ni } = this._N(th, ph), et = [Math.cos(th) * Math.cos(ph), Math.cos(th) * Math.sin(ph), -Math.sin(th)], ep = [-Math.sin(ph), Math.cos(ph), 0];
const nt = { re: dot(Nr, et), im: dot(Ni, et) }, np = { re: dot(Nr, ep), im: dot(Ni, ep) }, c = ETA0 * this.k * this.k / (32 * Math.PI * Math.PI);
return { U: c * (nt.re ** 2 + nt.im ** 2 + np.re ** 2 + np.im ** 2), Ut: c * (nt.re ** 2 + nt.im ** 2), Up: c * (np.re ** 2 + np.im ** 2) };
}Total radiated power, the strongest direction and the directivity.
P = Σ U sin θ Δθ Δφ; D = 4π U_max ÷ P
A midpoint sum on a 90 × 72 grid of directions (the upper half only over a ground). The peak is the largest of the sampled cells, so a narrow lobe between samples is missed. The ratio of radiated to input power is returned with the result.
/** Total radiated power by quadrature on the sphere (upper hemisphere only with a ground). Returns {P, Umax, Dmax, dirMax, Pin, balance}. */
pattern(nth = 90, nph = 72) {
const th1 = this.ground === 'perfect' ? Math.PI / 2 : Math.PI, dth = th1 / nth, dph = 2 * Math.PI / nph; let Prad = 0, Umax = -1, dir = [0, 0];
for (let i = 0; i < nth; i++) { const th = (i + 0.5) * dth; for (let j = 0; j < nph; j++) { const ph = (j + 0.5) * dph, u = this.intensity(th, ph).U; Prad += u * Math.sin(th) * dth * dph; if (u > Umax) { Umax = u; dir = [th, ph]; } } }
const Pin = 0.5 * (this.V.re * this.I.re[this.feedBasis] + this.V.im * this.I.im[this.feedBasis]); // ½ Re(V I*)
return { P: Prad, Umax, dirMax: dir, Dmax: 4 * Math.PI * Umax / Prad, Pin, balance: Prad / Pin };
}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.