Fix gshank convergence-floor regression, nec2diff comparison bugs (#138), and add issue #139 regression tests

CHANGELOG.md

@@ -3,6 +3,8 @@
 ### Bug Fixes
 * **`zero-current` now maps to signed NEC `GE -1`:** the stable public ABI value remains `2` (`NEC_GROUND_CONNECTION_ZERO_CURRENT`), while the native geometry receives `-1`. Previously `GE -1` was documented but effectively treated as no ground: `build_connections()` guarded every ground check with `ignd > 0`, so zero-current geometries skipped the plane checks the manual mandates. Both signed ground modes now reject segments that extend below the ground plane or lie in it (a segment may still end on the plane), and a non-none connection without a ground model (`GN` card) fails at `simulate()` instead of silently running free space — NEC-2 Part 3: "A positive or negative value of I1 does not cause a ground to be included in the calculation... The ground parameters must be specified on a program control card following the geometry cards." The `GE -1` connection semantics are unchanged from the Fortran: a ground-touching end stays a free end, so the current goes to zero there, and the low-horizontal-wire escape hatch (height < 1e-3 × segment length) still applies. `testharness/data/discone.nec` gained the missing `GN 1`; `test_example1` now completes the geometry with 0, matching its free-space deck.
 * **`nec2diff` no longer reports spurious radiation-pattern differences:** at angles where the polarization is undefined (HORIZ gain `-999.99`, e.g. at the horizon) NEC leaves the SENSE column blank. The `RadiationInput` parser read the next numeric token as the sense string, shifting the remaining columns left by one, and `read_fixed`/`read_sci` returned **uninitialized memory** when the stream was exhausted — so comparing a file against itself reported a nonzero difference with garbage values (observed on `bruce_sommerfeld`). The parser now detects the blank SENSE column, and the readers return 0.0 on extraction failure. All 52 testharness decks now self-compare exactly clean; genuine differences are still flagged.
+* **`gshank` Shanks-tail convergence floor restored after it caused non-terminating Sommerfeld integrations:** issue #131 item 3 also raised gshank()'s convergence floor from `denm = 1e-3*den*CRIT` to `2e-1*den*CRIT_H` — 200× looser. On decks with wires close to a `GN 2` ground (e.g. a radial screen elevated 0.2 m over epsr 13 / sig 0.005 at 2.7 MHz), the truncated Shanks tail degrades the Sommerfeld tables enough that `rom2()` can never meet its `RX = 1e-4` Romberg test: the integration doubles to the 65536-subdivision cap at every step, emits an endless `ROM2 -- STEP SIZE LIMITED` stream (GB/minute into both stdout and the output file) and effectively never terminates. The same deck completes in ~14 ms with the previous floor, agreeing with `nec2c` (13.906 − j445.84 Ω). The floor is restored to `1e-3*den*CRIT_H` while the milder `CRIT_H = 1.1e-4` tail-test scaling is kept; output is byte-identical (modulo timing lines) to the pre-loosening build on `sommerfeld2/3` and `bruce_sommerfeld`. Regression decks: `testharness/data/wire_order_normal.nec`, `wire_order_reversed.nec`; unit test `[issue139]` (which also fails on this regression via timeout).
+* **`nec2diff` compares negative values, reports finite sign-flip differences, and exits non-zero when results differ (#138):** the real-valued `diff()` in `BaseInput.h` tested *signed* values against the near-zero floor (`a < 1e-8 && b < 1e-8`), so any two negative numbers compared equal — a 178 Ω negative-reactance change (`−179.39 → −1.00`) or a `−1.71 → −9.71 dBi` gain change reported `Difference = 0`. A pure sign flip (`−179.39 → +179.39`) divided by `|a+b| = 0` and printed `inf`. Both `diff()` overloads now normalize by the larger magnitude (`|a−b|/max(|a|,|b|)`), keeping the near-zero early-out but keyed on absolute values. Two behavioral notes: (1) this metric is close to twice the old `|a−b|/|a+b|` for near-equal values, so the `equalto()` thresholds (sum ≤ 1e-4) now bind about twice as tightly and previously hidden last-digit engine discrepancies are visible in `nec2diff` output; (2) `main()` now tracks a flag across the three sections and returns 1 if any section differs, so scripts and harnesses can rely on the exit status. `nec2diff` with fewer than two arguments prints usage and returns 1 instead of throwing `std::logic_error` from a null `std::string` construction. Regression fixtures and an exit-code test live in `tests/data/issue138/` (ctest `necpp_nec2diff_issue138`); the four cases report 0.994, 0.824, 0.167 (control) and a finite 2.0, matching the reporter's values.
 
 ### Performance
 * **`saoa` shares one `xl*xl` across all five uses:** the Sommerfeld integrand squared `xl` independently for each of `cgam1`, `cgam2` (the two complex `sqrt` arguments), the three `dgam` large-|xl| branches, and the `rho==0` transverse answer. One complex multiply per integrand evaluation is now shared (issue #131). Algebraic identity — byte-identical output on all five `GN 2` decks; benchmarked at ≤0.1% on a noisy laptop, consistent with the multiplication being a small fraction of the transcendental-dominant kernel.

nec2diff.1

@@ -6,10 +6,15 @@ antennas using the Method of Moments.
 
 .SH SYNOPSIS
 .B nec2diff
+.I file_a file_b
 
 .SH DESCRIPTION
 .B nec2diff
 evaluates differences between 2 result files given as parameters.
 
+.SH EXIT STATUS
+nec2diff returns 0 if no differences were found, and 1 if the files differ.
+With fewer than two arguments it prints a usage line and returns 1.
+
 .SH AUTHOR
 nec2diff was written by Tim Molteno <tim@physics.otago.ac.nz>

src/BaseInput.h

@@ -19,23 +19,24 @@
 
 #include "math_util.h"
 
+#include <algorithm>
+#include <cmath>
 #include <fstream>
 #include <sstream>
 
+/* Relative difference of two values, normalised by the larger magnitude.
+   The scale test must use magnitudes: comparing signed values against the
+   1e-8 floor treats any two negative numbers as "effectively zero". */
 nec_float diff(nec_float a, nec_float b)
 {
 	if (a == b)
 		return 0;
-	if ((a < 1e-8) && (b < 1e-8))
+
+	nec_float scale = std::max(std::fabs(a), std::fabs(b));
+	if (scale < 1e-8)
 		return 0;
 
-	nec_float sub = a - b;
-	nec_float sum = a + b;
-
-	nec_float ret = sub*sub;
-	if ((a != 0) && (b != 0))
-		ret /= sum*sum;
-	ret = sqrt(ret);
+	nec_float ret = std::fabs(a - b) / scale;
 	if (ret > 1e-2)
 	{
 		cout << "    diff(" << a << "," << b << ") = " << ret << endl;
@@ -47,16 +48,12 @@ nec_float diff(nec_complex a, nec_complex b)
 {
 	if (a == b)
 		return 0;
-	if ((abs(a) < 1e-8) && (abs(b) < 1e-8))
+
+	nec_float scale = std::max(abs(a), abs(b));
+	if (scale < 1e-8)
 		return 0;
 
-	nec_float ret = norm(a - b);
-	nec_float sum = norm(a + b);
-
-	if (sum != 0)
-		ret /= sum;
-
-	ret = sqrt(ret);
+	nec_float ret = abs(a - b) / scale;
 	if (ret > 1e-2)
 	{
 		stringstream ss;

src/c_evlcom.cpp

@@ -36,7 +36,9 @@ using namespace std;
    before the step-level tolerance binds. Measured to cut ~1.5% of saoa()
    integrand evaluations per Sommerfeld table with byte-identical output on
    every testharness deck and unchanged agreement with the nec2c/FORTRAN
-   references. */
+   references. The companion change that raised gshank()'s denm floor was
+   reverted (see the denm assignment below) after it broke convergence on
+   decks with wires close to the ground. */
 #define CRIT_H	1.1E-4
 #define	PTP	.6283185308
 #define NTS	4
@@ -181,12 +183,15 @@ void c_evlcom::gshank( nec_complex start, nec_complex dela, nec_complex* sum,
 				den=amg;		
 		} /* for ( i = 0; i < nans; i++ ) */
 
-		/* Convergence floor for the Shanks test (issue #131 item 3): the
-		   increment is judged against max(|q2|*CRIT_H, denm); at 2.e-1 of
-		   the largest accelerated term the loop stops once the increment
-		   has stalled relative to the tail, instead of chasing entries
-		   orders of magnitude smaller than the answer. */
-		denm=2.e-1*den*CRIT_H;
+		/* Convergence floor for the Shanks test: the increment is judged
+		   against max(|q2|*CRIT_H, denm) with denm at 1.e-3 of the largest
+		   accelerated term. Issue #131 item 3 initially raised this floor
+		   to 2.e-1, which lets gshank() return while the integrand is still
+		   converging; the resulting errors in the Sommerfeld tables make
+		   rom2()'s Romberg test (rx = 1e-4) unreachable for decks with
+		   wires close to the ground (e.g. elevated radial screens), so the
+		   integration grinds at the subdivision cap indefinitely. */
+		denm=1.e-3*den*CRIT_H;
 		jm=intx-3;
 		if (jm < 1)
 			jm = 1;

src/necDiff.cpp

@@ -40,28 +40,33 @@ using namespace std;
 
 int main(int argc, char** argv) {
 
-  if (argc < 2) {
-    cout << "Usage: nec2diff fila, file b" << endl;
+  if (argc < 3) {
+    cout << "Usage: nec2diff file_a file_b" << endl;
+    return 1;
   }
-  cout << "NecDiff 0.12" << endl;
+
+  cout << "NecDiff 0.13" << endl;
 
   string file1(argv[1]);
   string file2(argv[2]);
 
   cout << "File1: " << file1 << endl;
   cout << "File2: " << file2 << endl;
 
+  bool same = true;
+
   {
     AntennaInput ai1(file1);
     AntennaInput ai2(file2);
 
     cout << "AntennaInput Difference = " << ai1.difference(ai2) << endl;
     if (!ai1.equalto(ai2)) {
+      same = false;
       cerr << file1 << "::" << file2;
       cerr << ". Input Parameters Different " << ai1.difference(ai2) <<  endl;
     }
   }
-  
+
   // Compare Power Budget
 
   {
@@ -70,9 +75,9 @@ int main(int argc, char** argv) {
 
     cout << "PowerBudget Difference = " << r1.difference(r2) << endl;
     if (!r1.equalto(r2)) {
+      same = false;
       cerr << file1 << "::" << file2;
       cerr << ". Power Budgets Different " << r1.difference(r2) << endl;
-      // return 1;
     }
   }
 
@@ -83,13 +88,13 @@ int main(int argc, char** argv) {
 
     cout << "RadiationInput Difference = " << r1.difference(r2) << endl;
     if (!r1.equalto(r2)) {
+      same = false;
       cerr << file1 << "::" << file2;
       cerr << ". Radiation Patterns Different " << r1.difference(r2) << endl;
-      // return 1;
     }
   }
-  
-  return 0;
+
+  return same ? 0 : 1;
 }
 
 

src/nec_context_tb.cpp

@@ -506,3 +506,82 @@ TEST_CASE( "EX type-4 source on a patch centroid zeroes that patch", "[etmns][cu
     REQUIRE(std::real(c1) == Catch::Approx(6.9547e-05).margin(1e-7));
     REQUIRE(std::imag(c1) == Catch::Approx(-1.5977e-05).margin(1e-7));
 }
+
+// ---------------------------------------------------------------------------
+// Issue #139: Sommerfeld interpolation produced wire-order-dependent results.
+//
+// A 12 m monopole with an explicit radial screen 0.2 m above a Sommerfeld
+// ground (GN 2, epsr 13, sig 0.005) is built twice: once in natural wire
+// order and once with the wire cards reversed. The antenna input impedance,
+// accepted power and radiation pattern must not depend on the insertion
+// order.
+//
+// This deck also guards the gshank() Shanks-tail convergence floor: when it
+// was loosened (issue #131 item 3), the Sommerfeld tables degraded enough
+// that rom2() could never meet its Romberg test for this geometry and the
+// simulation ground on at the subdivision cap indefinitely. A regression
+// there shows up as this test hitting the ctest timeout.
+// ---------------------------------------------------------------------------
+namespace {
+	struct issue139_wire {
+		int tag, segments;
+		double x1, y1, z1, x2, y2, z2, radius;
+	};
+
+	void issue139_run(bool reverse_wires,
+		nec_complex& impedance, double& accepted_power, double& horizon_gain_db) {
+		issue139_wire wires[] = {
+			{ 1, 13, 0.0, 0.0, 0.2, 0.0, 0.0, 12.2, 0.025 },
+			{ 2, 20, 0.0, 0.0, 0.2, 20.0, 0.0, 0.2, 0.001 },
+			{ 3, 20, 0.0, 0.0, 0.2, 0.0, 20.0, 0.2, 0.001 },
+		};
+
+		nec_context nec;
+		nec.initialize();
+		c_geometry* geo = nec.get_geometry();
+		for (int i = 0; i < 3; i++) {
+			const issue139_wire& w = wires[reverse_wires ? 2 - i : i];
+			geo->wire(w.tag, w.segments, w.x1, w.y1, w.z1, w.x2, w.y2, w.z2,
+				w.radius, 1.0, 1.0);
+		}
+		nec.geometry_complete(1);
+		nec.gn_card(2, 0, 13.0, 0.005, 0.0, 0.0, 0.0, 0.0);
+		nec.ld_card(5, 0, 0, 0, 58.0e6, 0, 0);
+		nec.fr_card(0, 1, 2.7, 0.0);
+		nec.ex_card(EXCITATION_VOLTAGE, 1, 1, 0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+		nec.rp_card(0, 2, 1, 1, 0, 0, 0, 0.0, 0.0, 89.5, 0.0, 0.0, 0.0);
+
+		impedance = nec_complex(nec.get_impedance_real(0), nec.get_impedance_imag(0));
+		accepted_power = 0.0;
+		nec_antenna_input* ai = nec.get_input_parameters(0);
+		for (unsigned k = 0; k < ai->get_power().size(); k++)
+			accepted_power += ai->get_power()[k];
+
+		nec_radiation_pattern* rp = nec.get_radiation_pattern(0);
+		// theta = 0 is the symmetry null (gain floor sentinel); the 89.5
+		// degree sample (index 1) is the meaningful order-invariance probe.
+		horizon_gain_db = rp->get_gain_tot()[1];
+	}
+}
+
+TEST_CASE("Wire insertion order does not change Sommerfeld-ground results", "[issue139][sommerfeld]") {
+	nec_complex z_forward(0.0, 0.0), z_reverse(0.0, 0.0);
+	double p_forward = 0.0, p_reverse = 0.0;
+	double g_forward = 0.0, g_reverse = 0.0;
+
+	issue139_run(false, z_forward, p_forward, g_forward);
+	issue139_run(true, z_reverse, p_reverse, g_reverse);
+
+	// Reference values (cross-validated against the nec2c reference engine,
+	// which prints 1.3906E+01 -4.4584E+02 for this deck).
+	REQUIRE(std::real(z_forward) == Catch::Approx(13.906).margin(1e-2));
+	REQUIRE(std::imag(z_forward) == Catch::Approx(-445.84).margin(1e-1));
+
+	// Order invariance: the two runs share identical arithmetic up to the
+	// segment numbering, so the comparison can be tight.
+	REQUIRE(std::abs(z_forward - z_reverse)
+		<= 1e-8 * std::max(std::abs(z_forward), std::abs(z_reverse)));
+	REQUIRE(std::abs(p_forward - p_reverse)
+		<= 1e-8 * std::max(std::abs(p_forward), std::abs(p_reverse)));
+	REQUIRE(std::abs(g_forward - g_reverse) < 1e-6);
+}

testharness/data/short_dipole.nec

@@ -0,0 +1,8 @@
+CM Short (capacitive) 20 m dipole in free space
+CE
+GW 1 21 0 -4.5 0 0 4.5 0 0.001
+GE 0
+EX 0 1 11 0 1 0
+FR 0 1 0 0 14.2 0
+RP 0 1 3 1000 90 0 0 45
+EN

testharness/data/wire_order_normal.nec

@@ -0,0 +1,12 @@
+CM Issue 139 regression: monopole + radial screen over Sommerfeld ground (normal wire order)
+CE
+GW 1 13 0 0 0.2 0 0 12.2 0.025
+GW 2 20 0 0 0.2 20 0 0.2 0.001
+GW 3 20 0 0 0.2 0 20 0.2 0.001
+GE 1
+GN 2,0,0,0,13.,.005
+LD 5 0 0 0 58000000 0 0
+EX 0 1 1 0 1 0
+FR 0 1 0 0 2.7 0
+RP 0 2 1 1001 0 0 89.5 0
+EN

testharness/data/wire_order_reversed.nec

@@ -0,0 +1,12 @@
+CM Issue 139 regression: same antenna, reversed wire insertion order
+CE
+GW 3 20 0 0 0.2 0 20 0.2 0.001
+GW 2 20 0 0 0.2 20 0 0.2 0.001
+GW 1 13 0 0 0.2 0 0 12.2 0.025
+GE 1
+GN 2,0,0,0,13.,.005
+LD 5 0 0 0 58000000 0 0
+EX 0 1 1 0 1 0
+FR 0 1 0 0 2.7 0
+RP 0 2 1 1001 0 0 89.5 0
+EN

tests/CMakeLists.txt

@@ -99,4 +99,12 @@ add_test(NAME necpp_smoke_hertzian_dipole
     COMMAND ${CMAKE_COMMAND}
         -DNEC2PP=$<TARGET_FILE:nec2++>
         -DINPUT=${CMAKE_SOURCE_DIR}/testharness/data/herzian_dipole.nec
-        -P ${CMAKE_CURRENT_SOURCE_DIR}/smoke_test.cmake)
+        -P ${CMAKE_CURRENT_SOURCE_DIR}/smoke_test.cmake)
+
+# nec2diff behaviour regression (issue #138): negative-value comparison,
+# finite sign-flip differences, exit statuses and the argc guard.
+add_test(NAME necpp_nec2diff_issue138
+    COMMAND ${CMAKE_COMMAND}
+        -DNEC2DIFF=$<TARGET_FILE:nec2diff>
+        -DDATA_DIR=${CMAKE_CURRENT_SOURCE_DIR}/data/issue138
+        -P ${CMAKE_CURRENT_SOURCE_DIR}/nec2diff_test.cmake)