From b765112316ee3f89958cf040acbec0f8c2f7e555 Mon Sep 17 00:00:00 2001 From: Eric Wheeler Date: Tue, 11 Aug 2026 17:49:52 -0700 Subject: [PATCH 1/3] fix: bind ROM2 subdivision cap to NEC-2's NM constant, not structure size nec_context::rom2 derived its minimum-interval epsilon and step-halving bound from the structure's segment-plus-patch count instead of the routine-local NM=65536 constant NEC-2 declares in SUBROUTINE ROM2. Small models over ground therefore received the coarsest Sommerfeld integration, exactly where per-segment ground interaction is strongest, and could exhaust the bound before the relative-error test converged, accumulating an unconverged Richardson extrapolation into the field sum. - declare nma = 65536 in both rom2 copies (#if 1 and #else), alongside the reference quartet each already carries - consume nma in place of m_geometry->n_plus_m() at the epsilon derivation and both step-halving bound checks in each copy, restoring the sixteen halvings NEC-2 permits regardless of geometry size - No control flow changes; hfk, intx, and rom1 already used the correct constant and are unaffected. Signed-off-by: Eric Wheeler --- src/nec_context.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/nec_context.cpp b/src/nec_context.cpp index e6fb98e..ddb7ffb 100644 --- a/src/nec_context.cpp +++ b/src/nec_context.cpp @@ -5638,7 +5638,7 @@ void nec_context::rom2( nec_float a, nec_float b, complex_array& sum, nec_float WRITE(3,18) STOP */ - int nts = 4, nx = 1; + int nts = 4, nx = 1, nma = 65536; const int n = 9; nec_float rx = 1e-4; nec_float dz = 0.0; @@ -5666,7 +5666,7 @@ void nec_context::rom2( nec_float a, nec_float b, complex_array& sum, nec_float complex_array t01(9), t10(9), t20(9), t11(9); /*! tolerance for hitting upper limit */ - const nec_float ep = _s/(1.0e4 * m_geometry->n_plus_m()); + const nec_float ep = _s/(1.0e4 * nma); /*! upper limit */ const nec_float zend = ze - ep; @@ -5846,7 +5846,7 @@ C GO TO 5 */ nt=0; - if ( ns >= m_geometry->n_plus_m() ) + if ( ns >= nma ) { nec_error_mode em(m_output); m_output.string("ROM2 -- STEP SIZE LIMITED AT Z = "); @@ -5894,7 +5894,7 @@ void nec_context::rom2( nec_float a, nec_float b, complex_array& sum, nec_float bool recalculate_fields = true; - int nts = 4, nx = 1, n = 9; + int nts = 4, nx = 1, n = 9, nma = 65536; /*! subinterval size */ nec_float dz=0.0; @@ -5918,7 +5918,7 @@ void nec_context::rom2( nec_float a, nec_float b, complex_array& sum, nec_float } /*! tolerance for hitting upper limit */ - const nec_float ep = _s/(1.0e4 * m_geometry->n_plus_m()); + const nec_float ep = _s/(1.0e4 * nma); /*! upper limit */ const nec_float zend = ze - ep; @@ -6015,7 +6015,7 @@ void nec_context::rom2( nec_float a, nec_float b, complex_array& sum, nec_float if ( tr > rx) { nt=0; - if ( ns <= m_geometry->n_plus_m() ) + if ( ns <= nma ) { // halve step size ns = ns*2; From d0ee6a6bc11ef9b7d81e95eda052d372ed077e68 Mon Sep 17 00:00:00 2001 From: Eric Wheeler Date: Tue, 11 Aug 2026 17:50:51 -0700 Subject: [PATCH 2/3] c_ggrid: restore INTRP guard branch order for cached interpolation c_ggrid::interpolate cached the 4x4 grid region coefficients around the previous request and reused them when the next point fell in the same region, but its opening guard had the two branches exchanged relative to the NEC-2 reference. A point at or above the cached origin recomputed unconditionally, while a point below the origin reused stale coefficients whenever the leftover indices happened to satisfy the proximity test, so interpolation returned different values for the same argument pair depending on request order. This perturbed every Sommerfeld ground-reflected field the interpolation grid supplies to the impedance matrix, flipping the feed resistance sign at low heights over ground where it is a near-cancellation of large reactive terms. - introduce a jump flag set only when the point lies below the cached region origin, matching SUBROUTINE INTRP: indices are computed solely in the complementary branch - combine jump with the four index proximity tests into one disjunction that decides recomputation, replacing the inverted two-branch logic Signed-off-by: Eric Wheeler --- src/c_ggrid.cpp | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/c_ggrid.cpp b/src/c_ggrid.cpp index 46b3a42..ceb8669 100644 --- a/src/c_ggrid.cpp +++ b/src/c_ggrid.cpp @@ -40,19 +40,23 @@ void c_ggrid::interpolate( nec_float x, nec_float y, nec_complex *f1, { static const int nda[3] = {11,17,9}, ndpa[3] = {110, 85, 72}; - bool recalculate = true; + /* A point below the cached region origin carries no valid index, so it + forces the coefficients to be rebuilt. */ + bool jump = false; - if( (x >= m_ip_xs) && (y >= m_ip_ys) ) { + if( (x < m_ip_xs) || (y < m_ip_ys) ) { + jump = true; + } else { m_ip_ix = static_cast((x - m_ip_xs) / m_ip_dx) + 1; m_ip_iy = static_cast((y - m_ip_ys) / m_ip_dy) + 1; - } else { - /* if point lies in same 4 by 4 point region */ - /* as previous point, old values are reused. */ - if ( ((m_ip_ix >= m_ip_ixeg) && (m_ip_iy >= m_ip_iyeg)) && - ((std::abs(m_ip_ix - m_ip_ixs) < 2) && (std::abs(m_ip_iy - m_ip_iys) < 2)) ) - recalculate = false; } + /* if point lies in same 4 by 4 point region */ + /* as previous point, old values are reused. */ + bool recalculate = jump || + (m_ip_ix < m_ip_ixeg) || (m_ip_iy < m_ip_iyeg) || + (std::abs(m_ip_ix - m_ip_ixs) >= 2) || (std::abs(m_ip_iy - m_ip_iys) >= 2); + if (true == recalculate) { /* determine correct grid and grid region */ int igr; From a89a6f9450b75087539f842294b5228ea37535b2 Mon Sep 17 00:00:00 2001 From: Eric Wheeler Date: Tue, 11 Aug 2026 18:01:44 -0700 Subject: [PATCH 3/3] fix: use squared-radial-distance norm2 for efld thin-wire displacement nec_context::efld displaces the observation point off the wire axis by the wire radius when the field point lies near the source segment during Sommerfeld ground evaluation. The displacement divides by the squared radial distance, but the code computed norm (a square root) instead of norm2 (the squared sum), taking the root twice and placing the field point at the wrong radius on every thin-wire ground evaluation. - replace norm(rhox, rhoy, rhoz) with norm2(rhox, rhoy, rhoz) in efld to match the squared-sum quantity the reference nec2c and nec2dx.f divide by, fixing the thin-wire radial displacement over Sommerfeld ground Signed-off-by: Eric Wheeler --- src/nec_context.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/nec_context.cpp b/src/nec_context.cpp index ddb7ffb..08e4d9e 100644 --- a/src/nec_context.cpp +++ b/src/nec_context.cpp @@ -2975,7 +2975,7 @@ void nec_context::efld( nec_float xi, nec_float yi, nec_float zi, nec_float ai, rhox= sabj* zij- salpr* yij; rhoy= salpr* xij- cabj* zij; rhoz= cabj* yij- sabj* xij; - rh = norm(rhox, rhoy, rhoz); // rhox* rhox+ rhoy* rhoy+ rhoz* rhoz; + rh = norm2(rhox, rhoy, rhoz); if ( rh <= 1.e-10) { xo= xi- ai* ysn;