diff --git a/CHANGELOG.md b/CHANGELOG.md index 8e50d64..3e9c43e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,4 +1,30 @@ -## Unreleased +## 0.9.0 - 2026-10-16 + +### Added + +* Conductivity-load (LD 5) regression coverage. Native and package tests pin + the following: + * `structureLossW` equals the segment-centre sum of NEC's `ZINT` resistance + in all three branches. + * Every selected segment is loaded with its own radius, feed segments + included. + * `prepare()` at a new frequency recomputes the load, on the stateful model + and on `NecArraySolver`. + * Perfect-ground loss counts physical segments only. + * The ring evaluator accepts loaded wires. + * Results agree with the original FORTRAN NEC-2 on copper, stainless, and + low-σ decks. + `docs/wasm-api.md` documents the formula and the accounting rule. + +### Bug Fixes + +* **Out-of-range load segments get an accurate error:** a tag-relative load + range that runs past the tag's last segment now fails at `addLoad` with + "SEGMENT TARGET IS OUT OF RANGE FOR THE TAG". The old message, "NO SEGMENT + HAS AN ITAG OF n", named a tag that does exist. The error class + (`NecGeometryError`) and the absent-tag message are unchanged. + +## 0.8.0 - 2026-09-27 ### Bug Fixes @@ -10,8 +36,6 @@ collision with upstream's C `nec_ground_connection`. Native C++ callers must update the type name and rebuild. C/WASM entry points and values are unchanged. -## 0.8.0 - 2026-09-27 - ### Added * Opt-in ring band-limited far-field evaluation for ordinary fields, both diff --git a/docs/wasm-api.md b/docs/wasm-api.md index 34babec..fce835d 100644 --- a/docs/wasm-api.md +++ b/docs/wasm-api.md @@ -120,6 +120,60 @@ tag-relative segment addressing follow the [EX card](https://www.nec2.org/part_3/cards/ex.html). Load units and targeting follow the [LD card](https://www.nec2.org/part_3/cards/ld.html). +### Conductivity loads (LD 5) + +A `{ kind: "conductivity", conductivitySPerM: σ }` load is NEC-2 LD type 5. +Perfect conductor is expressed by having no load. This section is normative +for consumers that compute conductor loss from exported currents. + +- **Scope.** Every segment the target selects is loaded, including source + (port) segments. A tag-only target selects every segment of the tag. +- **Radius.** Each segment uses its own radius, `radiiM` in current + distributions and `radiusM` in NECQ. +- **Impedance per metre.** NEC's `ZINT` is used, with \(\lambda\) the NEC + free-space wavelength and \(a\) the segment radius. Internally NEC uses + \(\omega\mu_0=2\pi\cdot120\pi/\lambda\), not \(\mu_0c_0\). Define + \(x=a\sqrt{2\pi\cdot120\pi\,\sigma/\lambda}=\sqrt2\,a/\delta\). Then + \[Z'=\frac{j}{a}\sqrt{\frac{60}{\sigma\lambda}}\; + \frac{\mathrm{ber}\,x+j\,\mathrm{bei}\,x}{\mathrm{ber}'x+j\,\mathrm{bei}'x} + \quad[\Omega/\mathrm{m}],\] + evaluated in three branches: + - \(x\le8\): NEC's Kelvin-function polynomials; + - \(8110\): the constant ratio \(0.70710678-0.70710678j\), which gives + \(R'=\sqrt{60/(\sigma\lambda)}\cdot0.70710678/a\). + + The first two branches match the exact \(J_0/J_1\) form to about + \(2\times10^{-7}\). The last branch drops the \(1/(\sqrt2\,x)\) correction, + so it reads \(0.35\%\) below the exact Bessel value at \(x\approx200\), and + \(R'\) jumps by about \(-0.64\%\) across \(x=110\). Consumers that must + reproduce `structureLossW` copy the three branches literally. The reference + implementation is `ld5ResistancePerMetre` in + `packages/necpp-wasm/test/conductivity-load.test.mjs`. +- **Frequency.** Each `prepare()` at a new frequency recomputes \(Z'\), on + both the stateful model and `NecArraySolver`. A model swept + \(f_1\to f_2\) is bitwise identical to a fresh model prepared at \(f_2\). +- **Loss accounting.** `structureLossW` is the segment-centre rule + \(\sum_i\tfrac12\,\mathrm{Re}(Z'_i)\,L_i\,|I_i(0)|^2\), where + \(I_i(0)=A_i+C_i\). It is not a quadrature of \(|I(s)|^2\) along the + segment. The sum runs over physical segments only, so perfect-ground images + are excluded. With four Gauss–Legendre NECQ nodes, \(I(0)\) is recovered + exactly from the even part \(A+C\cos ks\) at the two node pairs. The + Gauss–Legendre integral of \(|I|^2\) differs from `structureLossW` by about + \(10^{-5}\) to \(4\times10^{-4}\) on dipoles and helices. +- **Impedance shift.** \(\mathrm{Re}\,Z(\sigma)-\mathrm{Re}\,Z(\mathrm{PEC})\) + is not the conductor loss. The internal reactance \(\mathrm{Im}\,Z'\) + redistributes the current, which changes the radiated power by a comparable + amount: +9 % of the loss on a 0.47 λ copper dipole, −80 % on a helix. Use + `structureLossW` or the centre rule above. +- **Validation.** `addLoad` rejects \(\sigma\le0\) and non-finite values with + `NecInputError`. A tag-relative range that runs past the tag's last segment + fails with `NecGeometryError` "SEGMENT TARGET IS OUT OF RANGE FOR THE TAG"; + it is never silently clipped. Out-of-range absolute targets and absent tags + also fail at `addLoad`, not at `prepare`. +- **Field evaluators.** Loads do not change the field kernels. The ring + evaluator accepts loaded wires as ordinary wires. + ## Array layout and ownership `ComplexVector.real` and `.imag` have the same length. A vector supplied to a @@ -245,7 +299,7 @@ enumerates every operation/state pair plus both `prepare()` branches. | `addWire(wire)` | Positive integer tag/count; distinct finite endpoints and positive finite radius, all in m | `void`; copies the definition | `NecInputError` for shape/range errors; `NecGeometryError` for engine geometry limits | | `completeGeometry(options?)` | Ground connection plus optional finalized reflection/rotation descriptor | `GeometryCompletionResult`; `symmetry` is absent for ordinary completion | `NecInputError` for an invalid descriptor; `NecGeometryError` for intersections, invalid junctions, symmetry/ground conflicts, or a ground-incompatible structure | | `definePorts(ports)` | Nonempty ordered tag and one-based segment pairs | `void`; copies and freezes order | `NecPortError` for missing/duplicate ports or non-source-capable segments; `NecInputError` for malformed integers | -| `addLoad(load)` | Segment target and impedance, RLC, or conductivity values in the units above | `void`; invalidates prepared data | `NecInputError` for invalid values/ranges; `NecGeometryError` when no segment matches | +| `addLoad(load)` | Segment target and impedance, RLC, or conductivity values in the units above | `void`; invalidates prepared data | `NecInputError` for invalid values/ranges; `NecGeometryError` when no segment matches or a tag-relative range runs past the tag's last segment | | `clearLoads()` | None | `void`; removes every load and invalidates prepared data | No non-state failure | | `setGround(ground)` | Free space, perfect ground, or finite ground with method, relative permittivity, and S/m conductivity | `void`; invalidates prepared data | `NecInputError` for nonphysical values; `NecGeometryError` if inconsistent with geometry completion | | `prepare({ frequencyMHz })` | One positive finite MHz value | `void`; retains/factors the interaction matrix | `NecInputError` for frequency; `NecPortError` if ports are absent; `NecSolverError` for fill/factorization failure | diff --git a/packages/necpp-wasm/README.md b/packages/necpp-wasm/README.md index c47fa8a..5c9a6e6 100644 --- a/packages/necpp-wasm/README.md +++ b/packages/necpp-wasm/README.md @@ -864,7 +864,7 @@ appropriate CORS header. import { createNecModel } from "@necpp-engine/wasm"; const model = await createNecModel({ - wasmUrl: new URL("https://cdn.example.test/necpp/0.8.0/nec2pp.wasm"), + wasmUrl: new URL("https://cdn.example.test/necpp/0.9.0/nec2pp.wasm"), }); model.dispose(); ``` diff --git a/packages/necpp-wasm/fixtures/current-quadrature-v1/manifest.json b/packages/necpp-wasm/fixtures/current-quadrature-v1/manifest.json index 5c075c2..58aa715 100644 --- a/packages/necpp-wasm/fixtures/current-quadrature-v1/manifest.json +++ b/packages/necpp-wasm/fixtures/current-quadrature-v1/manifest.json @@ -3,7 +3,7 @@ "fixtureSchema": "current-quadrature-v1", "abiVersion": 1, "engineVersion": "2.5.0", - "packageVersion": "0.8.0", + "packageVersion": "0.9.0", "cases": [ { "id": "dipole", diff --git a/packages/necpp-wasm/package-lock.json b/packages/necpp-wasm/package-lock.json index b4379b4..9e32298 100644 --- a/packages/necpp-wasm/package-lock.json +++ b/packages/necpp-wasm/package-lock.json @@ -1,12 +1,12 @@ { "name": "@necpp-engine/wasm", - "version": "0.8.0", + "version": "0.9.0", "lockfileVersion": 3, "requires": true, "packages": { "": { "name": "@necpp-engine/wasm", - "version": "0.8.0", + "version": "0.9.0", "license": "GPL-2.0-or-later", "devDependencies": { "@types/node": "^24.13.3", diff --git a/packages/necpp-wasm/package.json b/packages/necpp-wasm/package.json index 65ed1fa..c4f1d8c 100644 --- a/packages/necpp-wasm/package.json +++ b/packages/necpp-wasm/package.json @@ -1,6 +1,6 @@ { "name": "@necpp-engine/wasm", - "version": "0.8.0", + "version": "0.9.0", "private": false, "type": "module", "description": "Stateful NEC2++ electromagnetic solver and symmetric-array API for Node and browsers", diff --git a/packages/necpp-wasm/src/versions.ts b/packages/necpp-wasm/src/versions.ts index 9267fcb..859867f 100644 --- a/packages/necpp-wasm/src/versions.ts +++ b/packages/necpp-wasm/src/versions.ts @@ -5,6 +5,6 @@ * CMake `project(necpp VERSION ...)` value compiled into the shipped WASM. * `abiVersion` is the stable C ABI prefix `necpp_wasm_v1`. */ -export const packageVersion = "0.8.0"; +export const packageVersion = "0.9.0"; export const abiVersion = 1; export const engineVersion = "2.5.0"; diff --git a/packages/necpp-wasm/test/conductivity-load.test.mjs b/packages/necpp-wasm/test/conductivity-load.test.mjs new file mode 100644 index 0000000..577c831 --- /dev/null +++ b/packages/necpp-wasm/test/conductivity-load.test.mjs @@ -0,0 +1,315 @@ +import assert from "node:assert/strict"; +import { existsSync } from "node:fs"; +import test from "node:test"; + +import { + NecGeometryError, + createNecArraySolver, + createNecModel, +} from "../.test-build/src/index.js"; + +const hasWasm = existsSync(new URL("../.test-build/src/nec2pp.wasm", import.meta.url)); + +const COPPER = 5.8e7; +const STAINLESS = 1.4e6; +const ONE_AMP = { real: Float64Array.of(1), imag: Float64Array.of(0) }; + +// Reference port of NEC-2 ZINT (src/nec_context.cpp, nec_context::zint) as +// applied by an LD 5 load: the real part of the internal impedance per metre +// of a round wire, in ohm/m. Constants, branch thresholds, and polynomial +// coefficients are NEC's own; see docs/wasm-api.md, "Conductivity loads". +function ld5ResistancePerMetre(sigmaSPerM, radiusM, wavelengthM) { + const sigl = sigmaSPerM * wavelengthM; + const rolam = radiusM / wavelengthM; + const x = Math.sqrt(2.368705e3 * sigl) * rolam; + // br1 = (ber + j bei) / (ber' + j bei'); Z'λ = j sqrt(60/sigl) br1 / rolam. + let br1; + if (x > 110) { + br1 = [0.70710678, -0.70710678]; + } else if (x > 8) { + const c = (re, im = 0) => [re, im]; + const add = (a, b) => [a[0] + b[0], a[1] + b[1]]; + const sub = (a, b) => [a[0] - b[0], a[1] - b[1]]; + const mul = (a, b) => [a[0] * b[0] - a[1] * b[1], a[0] * b[1] + a[1] * b[0]]; + const scale = (a, s) => [a[0] * s, a[1] * s]; + const exp = (a) => scale([Math.cos(a[1]), Math.sin(a[1])], Math.exp(a[0])); + const horner = (coefficients, d) => coefficients.slice(1) + .reduce((sum, value) => add(scale(sum, d), value), coefficients[0]); + const th = (d) => horner([c(6.0e-7, 1.9e-6), c(-3.4e-6, 5.1e-6), c(-2.52e-5), + c(-9.06e-5, -9.01e-5), c(0, -9.765e-4), c(0.0110486, -0.0110485), + c(0, -0.3926991)], d); + const ph = (d) => horner([c(1.6e-6, -3.2e-6), c(1.17e-5, -2.4e-6), + c(3.46e-5, 3.38e-5), c(5.0e-7, 2.452e-4), c(-1.3813e-3, 1.3811e-3), + c(-6.25001e-2, -1.0e-7), c(0.7071068, 0.7071068)], d); + const cn = c(0.7071068, 0.7071068); + const f = scale(exp(add(scale(cn, -x), th(-8 / x))), Math.sqrt((Math.PI / 2) / x)); + const g = scale(exp(add(scale(cn, x), th(8 / x))), 1 / Math.sqrt(2 * Math.PI * x)); + const jf = scale(mul([0, 1], f), 1 / Math.PI); + const numerator = add(g, jf); + const denominator = sub(mul(g, ph(8 / x)), mul(jf, ph(-8 / x))); + const d2 = denominator[0] ** 2 + denominator[1] ** 2; + br1 = [ + (numerator[0] * denominator[0] + numerator[1] * denominator[1]) / d2, + (numerator[1] * denominator[0] - numerator[0] * denominator[1]) / d2, + ]; + } else { + const y = (x / 8) ** 2; + const s = y * y; + const ber = ((((((-9.01e-6 * s + 1.22552e-3) * s - 0.08349609) * s + 2.6419140) * s + - 32.363456) * s + 113.77778) * s - 64) * s + 1; + const bei = ((((((1.1346e-4 * s - 0.01103667) * s + 0.52185615) * s - 10.567658) * s + + 72.817777) * s - 113.77778) * s + 16) * y; + const berp = (((((((-3.94e-6 * s + 4.5957e-4) * s - 0.02609253) * s + 0.66047849) * s + - 6.0681481) * s + 14.222222) * s - 4) * y) * x; + const beip = ((((((4.609e-5 * s - 3.79386e-3) * s + 0.14677204) * s - 2.3116751) * s + + 11.377778) * s - 10.666667) * s + 0.5) * x; + const d2 = berp ** 2 + beip ** 2; + br1 = [(ber * berp + bei * beip) / d2, (bei * berp - ber * beip) / d2]; + } + // Re(j * k * br1) = -k * Im(br1); divide λ·(Ω/m) by λ. + return -Math.sqrt(60 / sigl) * br1[1] / rolam / wavelengthM; +} + +// Segment-centre rule NEC uses for PowerBudget.structureLossW. +function centreRuleLossW(distribution, sigmaForTag) { + let loss = 0; + distribution.segments.forEach((segment, index) => { + const sigma = sigmaForTag(segment.tag); + if (!sigma) return; + const re = distribution.aReal[index] + distribution.cReal[index]; + const im = distribution.aImag[index] + distribution.cImag[index]; + loss += 0.5 * ld5ResistancePerMetre(sigma, distribution.radiiM[index], distribution.wavelengthM) + * distribution.lengthsM[index] * (re * re + im * im); + }); + return loss; +} + +const wavelength = (frequencyMHz) => 299.79563769321626 / frequencyMHz; + +async function loadedDipole({ frequencyMHz = 1000, sigma = COPPER, segments = 21 } = {}) { + const lambda = wavelength(frequencyMHz); + const model = await createNecModel(); + model.addWire({ + tag: 1, + segments, + start: [0, 0, -0.235 * lambda], + end: [0, 0, 0.235 * lambda], + radiusM: 0.001 * lambda, + }); + model.completeGeometry(); + model.definePorts([{ tag: 1, segment: (segments + 1) / 2 }]); + if (sigma) model.addLoad({ kind: "conductivity", target: { tag: 1 }, conductivitySPerM: sigma }); + model.setGround({ kind: "free-space" }); + return model; +} + +function helixWires(lambda) { + const radius = 0.16 * lambda; + const pitch = 0.23 * lambda; + const perTurn = 12; + const point = (index) => [ + radius * Math.cos(2 * Math.PI * index / perTurn), + radius * Math.sin(2 * Math.PI * index / perTurn), + 0.05 * lambda + pitch * index / perTurn, + ]; + const wires = [{ tag: 1, segments: 3, start: [radius, 0, 0], end: point(0), radiusM: 0.004 * lambda }]; + for (let index = 0; index < 3 * perTurn; index += 1) { + wires.push({ tag: index + 2, segments: 2, start: point(index), end: point(index + 1), radiusM: 0.004 * lambda }); + } + return wires; +} + +test("structureLossW is the LD5 segment-centre rule over every ZINT branch", + { skip: !hasWasm }, async () => { + for (const [sigma, frequencyMHz, label] of [ + [COPPER, 1000, "x > 110"], + [1.5e7, 1000, "8 < x <= 110 (brass)"], + [STAINLESS, 1000, "8 < x <= 110 (stainless)"], + [1e3, 1000, "x <= 8"], + ]) { + const model = await loadedDipole({ sigma, frequencyMHz }); + model.prepare({ frequencyMHz }); + const loss = model.solveCurrents(ONE_AMP).powerBudget.structureLossW; + const distribution = model.getCurrentDistribution({ kind: "latest-solution" }); + const expected = centreRuleLossW(distribution, () => sigma); + assert.ok(Math.abs(expected / loss - 1) < 1e-12, `${label}: ${expected} vs ${loss}`); + model.dispose(); + } + }); + +test("GL-4 NECQ samples reproduce structureLossW through the centre current", + { skip: !hasWasm }, async () => { + const nodes = Float64Array.of(-0.8611363115940526, -0.3399810435848563, + 0.3399810435848563, 0.8611363115940526); + const weights = Float64Array.of(0.34785484513745385, 0.6521451548625461, + 0.6521451548625461, 0.34785484513745385); + const lambda = wavelength(1000); + const model = await createNecModel(); + for (const wire of helixWires(lambda)) model.addWire(wire); + model.completeGeometry(); + model.definePorts([{ tag: 1, segment: 2 }]); + for (let tag = 1; tag <= 37; tag += 1) { + model.addLoad({ kind: "conductivity", target: { tag }, conductivitySPerM: COPPER }); + } + model.setGround({ kind: "free-space" }); + model.prepare({ frequencyMHz: 1000 }); + const characterization = model.characterizeIsolatedElement({ + quadrature: { nodes, weights, images: "physical-only", modes: "unit-current" }, + field: { theta: { startDeg: 0, count: 3, stepDeg: 90 }, phi: { startDeg: 0, count: 1, stepDeg: 0 } }, + }); + const view = new DataView(characterization.quadrature.buffer); + const nSegments = view.getUint32(12, true); + const nNodes = view.getUint32(16, true); + const wavelengthM = view.getFloat64(40, true); + let offset = 64 + 12 * nSegments; + offset += (8 - (offset % 8)) % 8; + const sampleCount = nSegments * nNodes; + const plane = (index) => new Float64Array(characterization.quadrature.buffer, + offset + 8 * sampleCount * index, sampleCount); + const radiusM = plane(6); + const lengthM = plane(7); + const dsWeight = plane(8); + const iReal = plane(9); + const iImag = plane(10); + const k = 2 * Math.PI / wavelengthM; + let centreRule = 0; + let gaussRule = 0; + for (let segment = 0; segment < nSegments; segment += 1) { + const base = segment * nNodes; + const r = ld5ResistancePerMetre(COPPER, radiusM[base], wavelengthM); + for (let node = 0; node < nNodes; node += 1) { + gaussRule += 0.5 * r * dsWeight[base + node] + * (iReal[base + node] ** 2 + iImag[base + node] ** 2); + } + // Even part of A + B sin(ks) + C cos(ks) at the two node pairs gives A + C. + const c1 = Math.cos(k * nodes[3] * lengthM[base] / 2); + const c2 = Math.cos(k * nodes[2] * lengthM[base] / 2); + const centre = [iReal, iImag].map((values) => { + const outer = (values[base] + values[base + 3]) / 2; + const inner = (values[base + 1] + values[base + 2]) / 2; + return outer + (outer - inner) * (1 - c1) / (c1 - c2); + }); + centreRule += 0.5 * r * lengthM[base] * (centre[0] ** 2 + centre[1] ** 2); + } + const native = model.solveCurrents(ONE_AMP).powerBudget.structureLossW; + assert.ok(Math.abs(centreRule / native - 1) < 1e-12, `${centreRule} vs ${native}`); + // The plain quadrature integral is a different, slightly larger rule here. + assert.ok(Math.abs(gaussRule / native - 1) > 1e-5); + assert.ok(Math.abs(gaussRule / native - 1) < 2e-3); + model.dispose(); + }); + +test("re-prepare at a new frequency recomputes LD5, stateful and array solver", + { skip: !hasWasm }, async () => { + const reused = await loadedDipole(); + reused.prepare({ frequencyMHz: 900 }); + reused.prepare({ frequencyMHz: 1100 }); + const zReused = reused.computeImpedanceMatrix().impedance; + const fresh = await loadedDipole(); + fresh.prepare({ frequencyMHz: 1100 }); + const zFresh = fresh.computeImpedanceMatrix().impedance; + assert.deepEqual([zReused.real[0], zReused.imag[0]], [zFresh.real[0], zFresh.imag[0]]); + reused.dispose(); + fresh.dispose(); + + const lambda = wavelength(1000); + const description = { + elements: [ + { id: "a", positionM: [-0.25 * lambda, 0], patternId: "dipole" }, + { id: "b", positionM: [0.25 * lambda, 0], patternId: "dipole" }, + ], + patterns: [{ + id: "dipole", + kind: "straight-wire-pattern", + wires: [{ + id: "radiator", + segments: 21, + startM: [0, 0, -0.235 * lambda], + endM: [0, 0, 0.235 * lambda], + radiusM: 0.001 * lambda, + }], + ports: [{ wireId: "radiator", segment: 11 }], + loads: [{ kind: "conductivity", target: { wireId: "radiator" }, conductivitySPerM: COPPER }], + }], + ground: { kind: "free-space" }, + }; + for (const symmetry of ["off", "auto"]) { + const options = { symmetry, symmetrizer: { positionEpsilonM: 1e-9 } }; + const sweep = await createNecArraySolver(description, options); + const single = await createNecArraySolver(description, options); + try { + await sweep.prepare({ frequencyMHz: 900 }); + await sweep.computeImpedanceMatrix(); + await sweep.prepare({ frequencyMHz: 1100 }); + const swept = (await sweep.computeImpedanceMatrix()).impedance; + await single.prepare({ frequencyMHz: 1100 }); + const direct = (await single.computeImpedanceMatrix()).impedance; + assert.deepEqual(Array.from(swept.real), Array.from(direct.real), symmetry); + assert.deepEqual(Array.from(swept.imag), Array.from(direct.imag), symmetry); + } finally { + await sweep.dispose(); + await single.dispose(); + } + } + }); + +test("perfect-ground structureLossW excludes image segments", + { skip: !hasWasm }, async () => { + const lambda = wavelength(1000); + const model = await createNecModel(); + model.addWire({ tag: 1, segments: 11, start: [0, 0, 0], end: [0, 0, 0.235 * lambda], radiusM: 0.001 * lambda }); + model.completeGeometry({ groundConnection: "interpolate" }); + model.definePorts([{ tag: 1, segment: 1 }]); + model.addLoad({ kind: "conductivity", target: { tag: 1 }, conductivitySPerM: COPPER }); + model.setGround({ kind: "perfect" }); + model.prepare({ frequencyMHz: 1000 }); + const loss = model.solveCurrents(ONE_AMP).powerBudget.structureLossW; + const distribution = model.getCurrentDistribution({ kind: "latest-solution" }); + assert.equal(distribution.segments.length, 11); + assert.ok(Math.abs(centreRuleLossW(distribution, () => COPPER) / loss - 1) < 1e-12); + model.dispose(); + }); + +test("addLoad rejects a segment range past the end of the tag", + { skip: !hasWasm }, async () => { + const model = await loadedDipole({ sigma: 0 }); + for (const target of [ + { tag: 1, firstSegment: 1, lastSegment: 22 }, + { tag: 1, firstSegment: 22 }, + ]) { + assert.throws( + () => model.addLoad({ kind: "conductivity", target, conductivitySPerM: COPPER }), + (error) => error instanceof NecGeometryError + && /OUT OF RANGE FOR THE TAG/.test(error.message), + ); + } + assert.throws( + () => model.addLoad({ kind: "conductivity", target: { tag: 7 }, conductivitySPerM: COPPER }), + (error) => error instanceof NecGeometryError && /NO SEGMENT HAS AN ITAG OF 7/.test(error.message), + ); + model.addLoad({ kind: "conductivity", target: { tag: 1, firstSegment: 1, lastSegment: 21 }, conductivitySPerM: COPPER }); + model.prepare({ frequencyMHz: 1000 }); + assert.ok(model.solveCurrents(ONE_AMP).powerBudget.structureLossW > 0); + model.dispose(); + }); + +test("the ring evaluator accepts loaded wires", + { skip: !hasWasm }, async () => { + const model = await loadedDipole({ sigma: STAINLESS }); + model.prepare({ frequencyMHz: 1000 }); + model.solveVoltages(ONE_AMP); + const request = { theta: { startDeg: 0, count: 91, stepDeg: 2 }, phi: { startDeg: 0, count: 180, stepDeg: 2 } }; + const exact = model.computeFarField(request); + const ring = model.computeFarField({ ...request, evaluator: "ring" }); + assert.equal(ring.fieldEvaluation?.execution, "ring"); + let delta = 0; + let scale = 0; + for (const key of ["eThetaReal", "eThetaImag", "ePhiReal", "ePhiImag"]) { + for (let index = 0; index < exact[key].length; index += 1) { + delta += (exact[key][index] - ring[key][index]) ** 2; + scale += exact[key][index] ** 2; + } + } + assert.ok(Math.sqrt(delta / scale) < 1e-9); + model.dispose(); + }); diff --git a/packages/necpp-wasm/test/facade-runtime.test.mjs b/packages/necpp-wasm/test/facade-runtime.test.mjs index cda017d..f552bd5 100644 --- a/packages/necpp-wasm/test/facade-runtime.test.mjs +++ b/packages/necpp-wasm/test/facade-runtime.test.mjs @@ -54,7 +54,7 @@ function addDipole(model) { } test("package, engine, and ABI versions are exported", () => { - assert.equal(packageVersion, "0.8.0"); + assert.equal(packageVersion, "0.9.0"); assert.equal(engineVersion, "2.5.0"); assert.equal(abiVersion, 1); }); diff --git a/packages/necpp-wasm/test/pack/manifest.test.mjs b/packages/necpp-wasm/test/pack/manifest.test.mjs index eb97293..f2bdb1a 100644 --- a/packages/necpp-wasm/test/pack/manifest.test.mjs +++ b/packages/necpp-wasm/test/pack/manifest.test.mjs @@ -23,7 +23,7 @@ test("npm pack contains only the documented publish files", { skip }, () => { assert.equal(packed.version, packageJson.version); const filenamePrefix = packageJson.name.slice(1).replace("/", "-"); assert.equal(packed.filename, `${filenamePrefix}-${packageJson.version}.tgz`); - assert.equal(packageJson.version, "0.8.0"); + assert.equal(packageJson.version, "0.9.0"); assert.equal(fixtureManifest.packageVersion, packageJson.version); assert.equal(packageJson.engines.node, ">=24"); assert.deepEqual(packageJson.publishConfig, { diff --git a/src/nec_stateful_model.cpp b/src/nec_stateful_model.cpp index b7a40ec..8fefbe7 100644 --- a/src/nec_stateful_model.cpp +++ b/src/nec_stateful_model.cpp @@ -387,9 +387,17 @@ void nec_stateful_model::validate_load_target( fail("ADD LOAD", "ABSOLUTE SEGMENT TARGET IS OUT OF RANGE"); } else { const int first = load.first_segment == 0 ? 1 : load.first_segment; - geometry->get_segment_number(load.tag, first); - if (load.last_segment != 0) - geometry->get_segment_number(load.tag, load.last_segment); + const int last = load.last_segment == 0 ? first : load.last_segment; + int64_t tag_segment_count = 0; + for (int64_t i = 0; i < geometry->n_segments; ++i) + if (geometry->segment_tags[i] == load.tag) + ++tag_segment_count; + // An absent tag keeps the legacy LD-card diagnostic. A present tag with a + // segment position past its end is reported as a range error instead. + if (tag_segment_count == 0) + geometry->get_segment_number(load.tag, first); + if (last > tag_segment_count) + fail("ADD LOAD", "SEGMENT TARGET IS OUT OF RANGE FOR THE TAG"); } } diff --git a/src/nec_stateful_model_conductivity_tb.cpp b/src/nec_stateful_model_conductivity_tb.cpp new file mode 100644 index 0000000..d50b977 --- /dev/null +++ b/src/nec_stateful_model_conductivity_tb.cpp @@ -0,0 +1,324 @@ +/* + Conductivity (LD 5) load regression tests for the stateful model. + + These pin the behaviour that downstream consumers rely on when they compute + conductor loss themselves from exported segment currents: + + - LD 5 applies NEC-2 ZINT per segment, using that segment's own radius, to + every selected segment (feed segments included). + - structure_loss_w is the segment-centre rule + sum 0.5 * Re(Z'_i) * L_i * |I(0)_i|^2 over physical segments only. + - prepare() at a new frequency recomputes the internal impedance. + - Out-of-range segment targets are rejected by add_load. + - Results agree with the original FORTRAN NEC-2 (nec2dx) on the same decks. +*/ +#include +#include + +#include "nec_exception.h" +#include "nec_stateful_model.h" + +#include +#include +#include +#include + +namespace { + +constexpr nec_float k_copper = 5.8e7; +constexpr nec_float k_dipole_half_length_m = 0.07045; +constexpr nec_float k_radius_m = 2.998e-4; +// NEC free-space wavelength at 1000 MHz (c0 = 299,795,637.69 m/s). +constexpr nec_float k_wavelength_1000_m = 0.299795637693; + +// ZINT's x > 110 branch: Re of j sqrt(60/(sigma lambda)) (0.70710678 - j 0.70710678) +// / (a/lambda), per metre. The truncated constant makes this 1.68e-9 below +// sqrt(30/(sigma lambda))/a. +nec_float skin_limit_r_per_metre(nec_float sigma, nec_float radius_m, + nec_float wavelength_m) +{ + return std::sqrt(60.0 / (sigma * wavelength_m)) * 0.70710678 / radius_m; +} + +nec_load_definition conductivity_load(int tag, nec_float sigma, + int first = 0, int last = 0) +{ + return {nec_load_kind::conductivity, tag, first, last, sigma, 0.0, 0.0}; +} + +void build_dipole(nec_stateful_model& model, int segments = 21) +{ + model.add_wire({1, segments, 0.0, 0.0, -k_dipole_half_length_m, + 0.0, 0.0, k_dipole_half_length_m, k_radius_m}); + model.complete_geometry(); + model.define_ports({{1, (segments + 1) / 2}}); +} + +void build_monopole(nec_stateful_model& model) +{ + model.add_wire({1, 11, 0.0, 0.0, 0.0, + 0.0, 0.0, k_dipole_half_length_m, k_radius_m}); + model.complete_geometry(nec_model_ground_connection::interpolate); + model.define_ports({{1, 1}}); +} + +nec_complex port_impedance(nec_stateful_model& model) +{ + const nec_impedance_result& result = model.compute_impedance_matrix(); + return result.impedance.at(0, 0); +} + +// Sum of 0.5 * L_i * |I(0)_i|^2 over the exported physical segments, with +// I(0) = A + C, the centre current NEC uses for its power budget. +nec_float half_l_centre_current_squared(const nec_current_distribution& d, + int only_tag = 0) +{ + nec_float sum = 0.0; + for (size_t i = 0; i < d.segment_count(); ++i) { + if (only_tag != 0 && d.segments[i].tag != only_tag) + continue; + const nec_complex centre = d.a_at(0, i) + d.c_at(0, i); + sum += 0.5 * d.lengths_m[i] * std::norm(centre); + } + return sum; +} + +// Measured per-metre resistance of a uniformly loaded single-radius wire. +nec_float measured_r_per_metre(nec_stateful_model& model) +{ + const nec_port_solution& solution = + model.solve_port_currents({nec_complex(1.0, 0.0)}); + const nec_float loss = solution.power_budget.structure_loss_w; + const nec_current_distribution d = + model.get_current_distribution(nec_current_mode_kind::latest_solution); + return loss / half_l_centre_current_squared(d); +} + +void require_close(nec_complex actual, nec_complex expected, nec_float rel) +{ + CAPTURE(actual, expected); + REQUIRE(std::abs(actual - expected) <= rel * std::abs(expected)); +} + +} // namespace + +TEST_CASE("LD5 per-metre resistance follows the NEC-2 ZINT branches", + "[conductivity]") +{ + struct row { + nec_float sigma; + nec_float expected_ohm_per_m; + nec_float rel_tol; + }; + // x = a * sqrt(omega mu0 sigma) with NEC's omega mu0 = 2 pi 120 pi / lambda. + // x > 110: ZINT uses the pure skin limit sqrt(30/(sigma lambda))/a. That + // omits the +1/(sqrt(2) x) term of the Bessel form (0.35 % low + // at x = 203), so the exact Bessel value is not the reference. + // x <= 110: polynomial/asymptotic Kelvin forms, within 2e-7 of the exact + // J0/J1 value (mpmath, 40 digits) listed here. + const std::vector rows = { + {5.8e7, 4.38129350226955, 1e-8}, // x = 202.95, copper + {1.5e7, 8.67463434860497, 1e-6}, // x = 103.21, brass + {1.4e6, 28.8432406473203, 1e-6}, // x = 31.53, stainless steel + {1.0e3, 3550.77958750132, 1e-6}, // x = 0.843 + }; + for (const row& r : rows) { + CAPTURE(r.sigma); + nec_stateful_model model; + build_dipole(model); + model.add_load(conductivity_load(1, r.sigma)); + model.prepare(1000.0); + const nec_current_distribution d = + model.get_current_distribution(nec_current_mode_kind::unit_current); + REQUIRE(d.wavelength_m == Catch::Approx(k_wavelength_1000_m).epsilon(1e-11)); + const nec_float measured = measured_r_per_metre(model); + REQUIRE(measured == Catch::Approx(r.expected_ohm_per_m).epsilon(r.rel_tol)); + if (r.sigma == k_copper) { + REQUIRE(measured == Catch::Approx( + skin_limit_r_per_metre(r.sigma, k_radius_m, d.wavelength_m)).epsilon(1e-13)); + } + } +} + +TEST_CASE("LD5 uses each segment's own radius and loads feed segments", + "[conductivity]") +{ + const nec_float lambda = k_wavelength_1000_m; + nec_stateful_model model; + model.add_wire({1, 7, 0.0, 0.0, -0.235 * lambda, 0.0, 0.0, -0.02 * lambda, 0.002 * lambda}); + model.add_wire({2, 1, 0.0, 0.0, -0.02 * lambda, 0.0, 0.0, 0.02 * lambda, 0.0006 * lambda}); + model.add_wire({3, 7, 0.0, 0.0, 0.02 * lambda, 0.0, 0.0, 0.235 * lambda, 0.001 * lambda}); + model.complete_geometry(); + model.define_ports({{2, 1}}); + + // All three radii are in the x > 110 branch for copper at 1 GHz + // (x = 406, 122, 203). + auto r_per_metre = [&](nec_float radius_m) { + return skin_limit_r_per_metre(k_copper, radius_m, lambda); + }; + + SECTION("only the feed segment") { + model.add_load(conductivity_load(2, k_copper)); + model.prepare(1000.0); + const nec_float loss = model.solve_port_currents({nec_complex(1.0, 0.0)}) + .power_budget.structure_loss_w; + const nec_current_distribution d = + model.get_current_distribution(nec_current_mode_kind::latest_solution); + REQUIRE(loss > 0.0); + REQUIRE(loss == Catch::Approx( + r_per_metre(0.0006 * lambda) * half_l_centre_current_squared(d, 2)) + .epsilon(1e-12)); + } + + SECTION("every wire, three radii") { + for (int tag = 1; tag <= 3; ++tag) + model.add_load(conductivity_load(tag, k_copper)); + model.prepare(1000.0); + const nec_float loss = model.solve_port_currents({nec_complex(1.0, 0.0)}) + .power_budget.structure_loss_w; + const nec_current_distribution d = + model.get_current_distribution(nec_current_mode_kind::latest_solution); + const nec_float expected = + r_per_metre(0.002 * lambda) * half_l_centre_current_squared(d, 1) + + r_per_metre(0.0006 * lambda) * half_l_centre_current_squared(d, 2) + + r_per_metre(0.001 * lambda) * half_l_centre_current_squared(d, 3); + REQUIRE(loss == Catch::Approx(expected).epsilon(1e-12)); + } +} + +TEST_CASE("prepare at a new frequency recomputes LD5 internal impedance", + "[conductivity][cache]") +{ + nec_stateful_model reused; + build_dipole(reused); + reused.add_load(conductivity_load(1, k_copper)); + reused.prepare(900.0); + const nec_complex z900 = port_impedance(reused); + reused.prepare(1100.0); + const nec_complex z1100_reused = port_impedance(reused); + const nec_float loss1100_reused = + reused.solve_port_currents({nec_complex(1.0, 0.0)}).power_budget.structure_loss_w; + reused.prepare(900.0); + const nec_complex z900_again = port_impedance(reused); + + nec_stateful_model fresh; + build_dipole(fresh); + fresh.add_load(conductivity_load(1, k_copper)); + fresh.prepare(1100.0); + const nec_complex z1100_fresh = port_impedance(fresh); + const nec_float loss1100_fresh = + fresh.solve_port_currents({nec_complex(1.0, 0.0)}).power_budget.structure_loss_w; + + nec_stateful_model pec; + build_dipole(pec); + pec.prepare(1100.0); + const nec_complex z1100_pec = port_impedance(pec); + + REQUIRE(z1100_reused == z1100_fresh); + REQUIRE(loss1100_reused == loss1100_fresh); + REQUIRE(z900_again == z900); + // Guard against a vacuous pass: copper visibly moves Re Z at 1100 MHz. + REQUIRE(z1100_fresh.real() - z1100_pec.real() > 0.3); +} + +TEST_CASE("perfect-ground structure loss counts physical segments only", + "[conductivity]") +{ + nec_stateful_model model; + build_monopole(model); + model.add_load(conductivity_load(1, k_copper)); + model.set_ground({nec_ground_kind::perfect, 0.0, 0.0}); + model.prepare(1000.0); + const nec_float loss = model.solve_port_currents({nec_complex(1.0, 0.0)}) + .power_budget.structure_loss_w; + const nec_current_distribution d = + model.get_current_distribution(nec_current_mode_kind::latest_solution); + REQUIRE(d.segment_count() == 11); + const nec_float r_per_metre = + skin_limit_r_per_metre(k_copper, k_radius_m, d.wavelength_m); + REQUIRE(loss == Catch::Approx(r_per_metre * half_l_centre_current_squared(d)) + .epsilon(1e-12)); +} + +TEST_CASE("conductivity loads agree with FORTRAN NEC-2 decks", + "[conductivity][reference]") +{ + // Reference: testharness/FORTRAN/nec2dx.f built with gfortran + // (-std=legacy -O0 -fno-automatic -finit-local-zero) and CVEL set to + // 299.79563769321626 so both engines use the same wavelength. Decks: + // GW 1 21 0 0 -0.07045 0 0 0.07045 0.0002998 / GE 0 / LD 5 1 0 0 + // FR 0 2 0 0 1000 100 / EX 0 1 11 0 1 0 / XQ + // and the monopole GW 1 11 0 0 0 0 0 0.07045 0.0002998 / GE 1 / GN 1 / + // ... / EX 0 1 1 0 1 0. FORTRAN prints 6 significant digits for Z and 5 + // for STRUCTURE LOSS (watts with a 1 V source). + struct row { + const char* name; + nec_float sigma; + bool monopole; + nec_float f_mhz; + nec_complex z; + nec_float loss_w; + }; + const std::vector rows = { + {"copper dipole", 5.8e7, false, 1000.0, {7.00946e+01, -7.99689e+00}, 3.2336e-05}, + {"copper dipole", 5.8e7, false, 1100.0, {9.56520e+01, 8.05099e+01}, 1.2257e-05}, + {"stainless dipole", 1.4e6, false, 1000.0, {7.20581e+01, -6.43996e+00}, 2.0287e-04}, + {"stainless dipole", 1.4e6, false, 1100.0, {9.80647e+01, 8.22555e+01}, 7.7112e-05}, + {"sigma 1e3 dipole", 1.0e3, false, 1000.0, {3.21389e+02, -5.35592e+01}, 1.1979e-03}, + {"sigma 1e3 dipole", 1.0e3, false, 1100.0, {3.72580e+02, 9.50813e+00}, 1.0198e-03}, + {"copper monopole", 5.8e7, true, 1000.0, {3.52834e+01, -3.66555e+00}, 6.4372e-05}, + {"copper monopole", 5.8e7, true, 1100.0, {4.69929e+01, 4.06593e+01}, 2.4374e-05}, + }; + for (const row& r : rows) { + CAPTURE(r.name, r.f_mhz); + nec_stateful_model model; + if (r.monopole) { + build_monopole(model); + model.set_ground({nec_ground_kind::perfect, 0.0, 0.0}); + } else { + build_dipole(model); + } + model.add_load(conductivity_load(1, r.sigma)); + model.prepare(r.f_mhz); + require_close(port_impedance(model), r.z, 2e-5); + const nec_float loss = model.solve_port_voltages_detailed({nec_complex(1.0, 0.0)}) + .power_budget.structure_loss_w; + REQUIRE(loss == Catch::Approx(r.loss_w).epsilon(1e-4)); + } +} + +TEST_CASE("add_load rejects segment targets outside the tag", + "[conductivity][validation]") +{ + auto message_of = [](const nec_load_definition& load) { + nec_stateful_model model; + build_dipole(model); + try { + model.add_load(load); + } catch (const nec_exception& error) { + return error.get_message(); + } + return std::string(); + }; + + REQUIRE(message_of(conductivity_load(1, k_copper)).empty()); + REQUIRE(message_of(conductivity_load(1, k_copper, 1, 21)).empty()); + REQUIRE(message_of(conductivity_load(1, k_copper, 21, 0)).empty()); + REQUIRE(message_of(conductivity_load(1, k_copper, 1, 22)) + .find("OUT OF RANGE FOR THE TAG") != std::string::npos); + REQUIRE(message_of(conductivity_load(1, k_copper, 22, 0)) + .find("OUT OF RANGE FOR THE TAG") != std::string::npos); + REQUIRE(message_of(conductivity_load(7, k_copper)) + .find("NO SEGMENT HAS AN ITAG OF") != std::string::npos); + REQUIRE(message_of(conductivity_load(0, k_copper, 1, 22)) + .find("ABSOLUTE SEGMENT TARGET IS OUT OF RANGE") != std::string::npos); + + // A rejected load leaves the model configurable and unloaded. + nec_stateful_model model; + build_dipole(model); + REQUIRE_THROWS_AS(model.add_load(conductivity_load(1, k_copper, 1, 22)), + nec_exception); + model.prepare(1000.0); + REQUIRE(model.solve_port_currents({nec_complex(1.0, 0.0)}) + .power_budget.structure_loss_w == 0.0); +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 4f48964..a48c802 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -46,6 +46,7 @@ set(NECPP_TEST_SRCS ${CMAKE_SOURCE_DIR}/src/nec_stateful_model_wp2_tb.cpp ${CMAKE_SOURCE_DIR}/src/nec_stateful_model_wp3_tb.cpp ${CMAKE_SOURCE_DIR}/src/nec_stateful_model_symmetry_tb.cpp + ${CMAKE_SOURCE_DIR}/src/nec_stateful_model_conductivity_tb.cpp ${CMAKE_SOURCE_DIR}/src/current_quadrature_wp0_tb.cpp ${CMAKE_SOURCE_DIR}/src/current_quadrature_wp1_tb.cpp ${CMAKE_SOURCE_DIR}/src/current_quadrature_wp2_tb.cpp