Coupling

The groundfield.coupling package collects every conductor-to-conductor and conductor-to-earth coupling relation that the solvers consume. The structure mirrors the physical decomposition: galvanic coupling through cluster nodes (handled by the solver builder, no module of its own), inductive coupling between distributed-conductor segments (inductance.py, ADR-0004), and the Carson earth-return correction that adds finite-conductivity effects to the inductive matrix (carson.py, ADR-0005).

Inductance — Neumann self and mutual integrals (ADR-0004)

Implements the Neumann partial-inductance assembly under a perfect-mirror earth assumption. Each pair of distributed-conductor segments contributes a self- or mutual-inductance entry to the branch-impedance block

\[ Z_b(\omega) \;=\; R \;+\; j\omega\, L_\text{Neumann}. \]

The thin-wire self-inductance uses Grover 1946 (closed form, plus optional internal-field \(\mu_0/(8\pi)\) contribution for the DC limit). Off-diagonal entries are evaluated by a closed-form parallel-segments fast path or, for arbitrary 3-D geometries, by a 16×16 Gauss–Legendre quadrature of the Neumann double-line integral. The image of every segment in the soil surface contributes one extra Neumann integral against the original — this is the perfect-mirror assumption and the starting point for the Carson refinement.

inductance

Self- and mutual-inductance helpers for distributed conductors.

This module implements the inductance physics described in ADR-0004. It provides three public functions:

• :func:thin_wire_self_inductance — closed-form self-inductance of a straight thin-wire segment of length \(\ell\) and radius \(a\) (Grover 1946): \(L_\text{self} = (\mu_0\,\ell)/(2\pi) \bigl[\ln(2\ell/a) - 1\bigr]\).
• :func:parallel_segments_mutual — closed-form mutual inductance of two parallel coaxial segments of equal length \(\ell\) at perpendicular distance \(d\).
• :func:neumann_mutual — generic Neumann double-line integral via two-point Gauss–Legendre quadrature, valid for any 3-D segment-pair geometry.

Earth-image contributions are handled by the caller: for a perfect-mirror earth the caller adds the integral against the mirror image of the source segment (see :func:perfect_mirror_self_pair_inductance and the assembly logic in :func:build_inductance_matrix).

References

• Grover, F. W. (1946). Inductance Calculations: Working Formulas and Tables. Dover (reprint 2004).
• Sunde, E. D. (1968). Earth Conduction Effects in Transmission Systems. Dover, ch. 7.
• Paul, C. R. (2010). Inductance: Loop and Partial. Wiley.

build_carson_correction_matrix

build_carson_correction_matrix(
    seg_endpoints: np.ndarray,
    wire_radii: np.ndarray,
    *,
    omega: float,
    sigma_earth: float
) -> np.ndarray

Assemble the dense Carson earth-return correction matrix.

Adds the per-segment-pair impedance correction \(\Delta Z^{(i,j)}_\text{Carson}(\omega)\) described in ADR-0005 on top of the perfect-mirror Neumann inductance matrix. The output has the same shape as :func:build_inductance_matrix but is complex and frequency-dependent — it is therefore rebuilt at every frequency by the solver.

Tiling convention (audit 2026-07-08, WP-C3)

Carson's correction is an infinite-wire, per-unit-length concept. To make the assembled matrix consistent — and, crucially, convergent under mesh refinement — the per-pair entries tile the per-unit-length values without double counting:

• Diagonal dZ[i, i]: self correction per metre × segment length. Summed over a refined straight wire this reproduces the full z'_self · ℓ exactly, independent of the segment count.
• Collinear pairs (same axis, negligible perpendicular offset): zero. The infinite-wire self value on the diagonal already contains the interaction with the rest of the same wire. The historic implementation additionally added a "mutual" term per collinear pair — with the axial midpoint distance misread as Carson's perpendicular wire separation — which made the assembled correction grow roughly linearly with the segment count (measured in the audit: no mesh convergence).
• Near-parallel pairs on distinct axes (|cosφ| > 0.99): mutual correction per metre at the perpendicular horizontal axis separation, weighted by the axial overlap length of the two segment intervals. Summing over the segments of the other wire recovers z'_mutual · ℓ_i exactly — refinement- invariant. Non-overlapping (e.g. staggered or end-to-end) pairs contribute zero: the long-range log tail is treated as part of the infinite-wire per-metre value, in line with the transmission-line locality of Carson's derivation.
• Oblique pairs (1e-9 < |cosφ| ≤ 0.99): the historic projected geometric-mean heuristic is retained (Carson's derivation only covers parallel wires; oblique pairs are rare in Manhattan-routed networks). Orthogonal pairs contribute zero by symmetry.

Earth-conductivity sign convention

In groundfield the \(z\)-axis points into the soil. A wire above ground has \(z < 0\) and height = -z > 0. A wire just below the surface has \(z > 0\) and we use height = z as the Sunde-equivalent depth (Carson's \(h\)). For the PEN cable at \(z = 0.6\,\mathrm{m}\) this produces the textbook "1 m below surface" Carson result.

Parameters:

Name	Type	Description	Default
seg_endpoints	ndarray	Array of shape (M, 2, 3) — same layout as :func:build_inductance_matrix.	required
wire_radii	ndarray	Per-branch wire radii in metres (currently only used as a regulariser when two segments would coincide; the Carson kernel itself does not depend on the radius — that is in the perfect-mirror \(\ln(2h/a)\) piece, which stays where it was in :func:build_inductance_matrix).	required
omega	float	Angular frequency in rad/s.	required
sigma_earth	float	Earth conductivity in S/m.	required

Returns:

Name	Type	Description
dZ	np.ndarray, shape (M, M), dtype complex	Symmetric Carson correction matrix in \(\Omega\). Entry dZ[i, j] is the integrated Carson contribution to the \((i, j)\) branch-impedance entry (per-unit-length Carson value times the geometric mean of the two segment lengths in the parallel-projection sense).

Source code in src/groundfield/coupling/inductance.py

def build_carson_correction_matrix(
    seg_endpoints: np.ndarray,        # shape (M, 2, 3)
    wire_radii: np.ndarray,           # shape (M,)
    *,
    omega: float,
    sigma_earth: float,
) -> np.ndarray:
    """Assemble the dense Carson earth-return correction matrix.

    Adds the per-segment-pair impedance correction
    $\\Delta Z^{(i,j)}_\\text{Carson}(\\omega)$ described in ADR-0005
    on top of the perfect-mirror Neumann inductance matrix. The
    output has the same shape as :func:`build_inductance_matrix`
    but is **complex** and **frequency-dependent** — it is therefore
    rebuilt at every frequency by the solver.

    Tiling convention (audit 2026-07-08, WP-C3)
    -------------------------------------------
    Carson's correction is an *infinite-wire, per-unit-length*
    concept. To make the assembled matrix consistent — and, crucially,
    **convergent under mesh refinement** — the per-pair entries tile
    the per-unit-length values without double counting:

    - **Diagonal** ``dZ[i, i]``: self correction per metre × segment
      length. Summed over a refined straight wire this reproduces the
      full ``z'_self · ℓ`` exactly, independent of the segment count.
    - **Collinear pairs** (same axis, negligible perpendicular
      offset): **zero.** The infinite-wire self value on the diagonal
      already contains the interaction with the rest of the same
      wire. The historic implementation additionally added a
      "mutual" term per collinear pair — with the axial midpoint
      distance misread as Carson's perpendicular wire separation —
      which made the assembled correction grow roughly linearly with
      the segment count (measured in the audit: no mesh convergence).
    - **Near-parallel pairs on distinct axes** (``|cosφ| > 0.99``):
      mutual correction per metre at the **perpendicular** horizontal
      axis separation, weighted by the **axial overlap length** of
      the two segment intervals. Summing over the segments of the
      other wire recovers ``z'_mutual · ℓ_i`` exactly — refinement-
      invariant. Non-overlapping (e.g. staggered or end-to-end)
      pairs contribute zero: the long-range log tail is treated as
      part of the infinite-wire per-metre value, in line with the
      transmission-line locality of Carson's derivation.
    - **Oblique pairs** (``1e-9 < |cosφ| ≤ 0.99``): the historic
      projected geometric-mean heuristic is retained (Carson's
      derivation only covers parallel wires; oblique pairs are rare
      in Manhattan-routed networks). Orthogonal pairs contribute
      zero by symmetry.

    Earth-conductivity sign convention
    ----------------------------------
    In ``groundfield`` the $z$-axis points *into* the soil. A wire
    above ground has $z < 0$ and ``height = -z > 0``. A wire just
    below the surface has $z > 0$ and we use ``height = z`` as the
    Sunde-equivalent depth (Carson's $h$). For the PEN cable at
    $z = 0.6\\,\\mathrm{m}$ this produces the textbook
    "1 m below surface" Carson result.

    Parameters
    ----------
    seg_endpoints
        Array of shape ``(M, 2, 3)`` — same layout as
        :func:`build_inductance_matrix`.
    wire_radii
        Per-branch wire radii in metres (currently only used as a
        regulariser when two segments would coincide; the Carson
        kernel itself does not depend on the radius — that is in
        the perfect-mirror $\\ln(2h/a)$ piece, which stays where
        it was in :func:`build_inductance_matrix`).
    omega
        Angular frequency in rad/s.
    sigma_earth
        Earth conductivity in S/m.

    Returns
    -------
    dZ : np.ndarray, shape (M, M), dtype complex
        Symmetric Carson correction matrix in $\\Omega$. Entry
        ``dZ[i, j]`` is the *integrated* Carson contribution to
        the $(i, j)$ branch-impedance entry (per-unit-length
        Carson value times the geometric mean of the two segment
        lengths in the parallel-projection sense).
    """
    from groundfield.coupling.carson import carson_mutual_correction

    if omega <= 0.0 or sigma_earth <= 0.0:
        return np.zeros(
            (seg_endpoints.shape[0], seg_endpoints.shape[0]),
            dtype=complex,
        )

    M = seg_endpoints.shape[0]
    dZ = np.zeros((M, M), dtype=complex)

    # Pre-compute mid-points, axis directions, lengths, heights.
    midpoints = 0.5 * (seg_endpoints[:, 0, :] + seg_endpoints[:, 1, :])
    axes = seg_endpoints[:, 1, :] - seg_endpoints[:, 0, :]
    lengths = np.linalg.norm(axes, axis=1)
    # Avoid divide-by-zero — caller is expected to filter degenerate
    # branches, but be defensive.
    safe_lengths = np.where(lengths > 0.0, lengths, 1.0)
    unit_axes = axes / safe_lengths[:, None]
    # Carson "height" h is the absolute value of the depth coordinate
    # (Sunde-equivalent). For a wire at the surface (z = 0) the
    # asymptote is logarithmic, regularised by the wire radius.
    heights = np.abs(midpoints[:, 2])
    heights = np.where(heights > 0.0, heights, np.maximum(wire_radii, 1e-3))

    for i in range(M):
        # Diagonal: self correction (parallel-to-self ⇒ θ = 0).
        a_self = 2.0 * heights[i] * math.sqrt(omega * MU_0 * sigma_earth)
        # Use the mutual function with d = 0 to share one code path
        # (it dispatches to θ = 0 internally). Per-unit-length value
        # times this branch's length.
        dz_per_m = carson_mutual_correction(
            omega=omega,
            height_i=heights[i],
            height_j=heights[i],
            horizontal_distance=0.0,
            sigma_earth=sigma_earth,
        )
        dZ[i, i] = dz_per_m * lengths[i]
        for j in range(i + 1, M):
            dot = float(unit_axes[i] @ unit_axes[j])
            if abs(dot) < 1e-9:
                continue  # purely orthogonal — zero by symmetry
            sign = 1.0 if dot > 0.0 else -1.0
            d_vec = midpoints[j] - midpoints[i]

            if abs(dot) > 0.99:
                # Near-parallel: transmission-line tiling (WP-C3).
                u = unit_axes[i]
                axial = float(d_vec @ u)
                perp = d_vec - axial * u
                d_perp_horiz = math.hypot(perp[0], perp[1])
                d_perp_vert = abs(perp[2])
                same_axis = (
                    d_perp_horiz < max(wire_radii[i], wire_radii[j], 1e-3)
                    and d_perp_vert < 1e-6
                )
                if same_axis:
                    # Same wire: earth term fully carried by the
                    # per-length diagonals — adding a "mutual" here
                    # double-counts (historic non-convergence bug).
                    continue
                # Axial overlap of the two segment intervals along u.
                t_i0 = float((seg_endpoints[i, 0] - midpoints[i]) @ u)
                t_i1 = float((seg_endpoints[i, 1] - midpoints[i]) @ u)
                lo_i, hi_i = min(t_i0, t_i1), max(t_i0, t_i1)
                t_j0 = float((seg_endpoints[j, 0] - midpoints[i]) @ u)
                t_j1 = float((seg_endpoints[j, 1] - midpoints[i]) @ u)
                lo_j, hi_j = min(t_j0, t_j1), max(t_j0, t_j1)
                overlap = min(hi_i, hi_j) - max(lo_i, lo_j)
                if overlap <= 0.0:
                    continue  # staggered/end-to-end: no local tiling
                dz_per_m = carson_mutual_correction(
                    omega=omega,
                    height_i=heights[i],
                    height_j=heights[j],
                    horizontal_distance=d_perp_horiz,
                    sigma_earth=sigma_earth,
                )
                dZ[i, j] = dZ[j, i] = sign * dz_per_m * overlap
                continue

            # Oblique pair: historic projected geometric-mean
            # heuristic (rare in Manhattan-routed networks).
            d_horiz = math.hypot(d_vec[0], d_vec[1])
            dz_per_m = carson_mutual_correction(
                omega=omega,
                height_i=heights[i],
                height_j=heights[j],
                horizontal_distance=d_horiz,
                sigma_earth=sigma_earth,
            )
            ell = math.sqrt(lengths[i] * lengths[j]) * abs(dot)
            dZ[i, j] = dZ[j, i] = sign * dz_per_m * ell
    return dZ

build_inductance_matrix

build_inductance_matrix(
    seg_endpoints: np.ndarray,
    wire_radii: np.ndarray,
    *,
    use_image: bool = True
) -> np.ndarray

Assemble the dense partial-inductance matrix over branches.

Vectorised implementation per ADR-0010 Tier 0b. Reproduces :func:_build_inductance_matrix_loop bit-exactly to floating- point precision but evaluates the off-diagonal entries one row at a time in batched NumPy calls — typical networks (~1000 segments) speed up by 1–2 orders of magnitude.

Parameters:

Name	Type	Description	Default
seg_endpoints	ndarray	Array of shape (M, 2, 3) with the start- and end-points of every distributed-conductor longitudinal-branch segment.	required
wire_radii	ndarray	Per-branch wire radii in metres, shape (M,).	required
use_image	bool	When True (default) the earth is treated as a perfect magnetic mirror: every branch's image (z → -z) is summed into both the diagonal (via :func:perfect_mirror_self_pair_inductance) and the off-diagonal (one extra Neumann integral against the mirrored partner).	True

Returns:

Name	Type	Description
L	(ndarray, shape(M, M))	Symmetric partial-inductance matrix in H. Entry L[i, j] is the partial mutual inductance between branch i and branch j (and the diagonal is the self-inductance + image).

Notes

The legacy loop-based implementation is kept as :func:_build_inductance_matrix_loop for regression testing. See :file:tests/test_inductance_vectorised.py for the bit-exact regression suite.

Source code in src/groundfield/coupling/inductance.py

def build_inductance_matrix(
    seg_endpoints: np.ndarray,        # shape (M, 2, 3): each branch's [start, end]
    wire_radii: np.ndarray,           # shape (M,)
    *,
    use_image: bool = True,
) -> np.ndarray:
    """Assemble the dense partial-inductance matrix over branches.

    Vectorised implementation per ADR-0010 Tier 0b. Reproduces
    :func:`_build_inductance_matrix_loop` bit-exactly to floating-
    point precision but evaluates the off-diagonal entries one row
    at a time in batched NumPy calls — typical networks
    (~1000 segments) speed up by 1–2 orders of magnitude.

    Parameters
    ----------
    seg_endpoints
        Array of shape ``(M, 2, 3)`` with the start- and end-points
        of every distributed-conductor longitudinal-branch segment.
    wire_radii
        Per-branch wire radii in metres, shape ``(M,)``.
    use_image
        When ``True`` (default) the earth is treated as a perfect
        magnetic mirror: every branch's image (z → -z) is summed
        into both the diagonal (via
        :func:`perfect_mirror_self_pair_inductance`) and the
        off-diagonal (one extra Neumann integral against the
        mirrored partner).

    Returns
    -------
    L : np.ndarray, shape (M, M)
        Symmetric partial-inductance matrix in H. Entry ``L[i, j]``
        is the partial mutual inductance between branch *i* and
        branch *j* (and the diagonal is the self-inductance + image).

    Notes
    -----
    The legacy loop-based implementation is kept as
    :func:`_build_inductance_matrix_loop` for regression testing.
    See :file:`tests/test_inductance_vectorised.py` for the
    bit-exact regression suite.
    """
    M = seg_endpoints.shape[0]
    if M == 0:
        return np.zeros((0, 0), dtype=float)

    L = np.zeros((M, M), dtype=float)

    # Pre-compute per-segment quantities that the loop would
    # otherwise rebuild on every call.
    p1 = seg_endpoints[:, 0, :].astype(float, copy=False)        # (M, 3)
    p2 = seg_endpoints[:, 1, :].astype(float, copy=False)        # (M, 3)
    da = p2 - p1                                                  # (M, 3)
    la = np.linalg.norm(da, axis=1)                               # (M,)
    if (la <= 0.0).any():
        raise ValueError("segment lengths must be positive")
    ua = da / la[:, None]                                         # (M, 3)

    # Mirrored endpoints (for the image contribution). Mirror by
    # flipping z; same shape and tangent direction up to a sign — the
    # Neumann integrand handles the sign via the dot product.
    if use_image:
        # Guard against the degenerate mirror geometry: a segment
        # lying (almost) in the soil-surface plane z = 0 coincides
        # with its own image, and the partial self-inductance
        # diverges logarithmically as the mirror distance 2|z|
        # approaches the wire radius. Without this check the
        # assembly silently produces self-inductances that are
        # orders of magnitude too large (observed: ~4.7 H instead
        # of ~55 µH for a 25 m segment), which in turn blocks the
        # branch current in the nodal solve. Conductors strung
        # between electrode heads at z = 0 are the typical trap.
        z_mid = 0.5 * np.abs(p1[:, 2] + p2[:, 2])
        degenerate = 2.0 * z_mid <= wire_radii
        if degenerate.any():
            k = int(np.argmax(degenerate))
            raise ValueError(
                f"build_inductance_matrix: segment {k} lies in the "
                f"soil-surface plane (|z_mid| = {z_mid[k]:.4g} m <= "
                f"wire_radius/2 = {0.5 * wire_radii[k]:.4g} m). Its "
                "perfect-mirror image coincides with the segment and "
                "the partial self-inductance diverges. Bury the "
                "conductor (e.g. z = 0.7 m) or disable the image "
                "term with use_image=False."
            )
        p1_img = p1.copy(); p1_img[:, 2] = -p1_img[:, 2]
        p2_img = p2.copy(); p2_img[:, 2] = -p2_img[:, 2]
        da_img = p2_img - p1_img
        la_img = np.linalg.norm(da_img, axis=1)
        ua_img = da_img / la_img[:, None]

    # Diagonal — self plus own-image. Cheap (M operations); leave it
    # as a Python loop so we can call the existing scalar helpers.
    for i in range(M):
        if use_image:
            L[i, i] = perfect_mirror_self_pair_inductance(
                p1[i], p2[i], wire_radii[i],
            )
        else:
            L[i, i] = thin_wire_self_inductance(la[i], wire_radii[i])

    # Off-diagonal — for each row i, compute L[i, i+1:M] as a batch.
    for i in range(M - 1):
        slc = slice(i + 1, M)
        K = M - i - 1

        ua_a = ua[i]                                              # (3,)
        la_a = float(la[i])
        da_a = da[i]                                              # (3,)
        p1_a = p1[i]                                              # (3,)

        ua_b = ua[slc]                                            # (K, 3)
        ub_b = ua_b   # alias for clarity below
        la_b = la[slc]                                            # (K,)
        db_b = da[slc]                                            # (K, 3)
        p1_b = p1[slc]                                            # (K, 3)

        # Pairwise dots between segment i and each j > i
        dots = ub_b @ ua_a                                        # (K,)

        # ---- Primary contribution: a vs. b ----
        # Closed form for parallel pairs, quadrature for the rest.
        M_par, par_mask = _parallel_filaments_mutual_batch(
            p1_a, ua_a, la_a, p1_b, ub_b, la_b, dots,
        )
        M_pri = np.where(par_mask, M_par, 0.0)
        npar_mask = ~par_mask
        if npar_mask.any():
            ks = np.where(npar_mask)[0]
            M_quad = _neumann_quadrature_batch(
                p1_a, da_a, la_a, ua_a,
                p1_b[ks], db_b[ks], la_b[ks], ub_b[ks],
                dots[ks],
            )
            M_pri[ks] = M_quad

        # ---- Image contribution: a vs. mirror(b) ----
        if use_image:
            ua_b_img = ua_img[slc]                                # (K, 3)
            la_b_img = la_img[slc]                                # (K,)
            db_b_img = (p2_img[slc] - p1_img[slc])                # (K, 3)
            p1_b_img = p1_img[slc]                                # (K, 3)
            dots_img = ua_b_img @ ua_a                            # (K,)

            M_par_img, par_mask_img = _parallel_filaments_mutual_batch(
                p1_a, ua_a, la_a,
                p1_b_img, ua_b_img, la_b_img, dots_img,
            )
            M_img = np.where(par_mask_img, M_par_img, 0.0)
            npar_mask_img = ~par_mask_img
            if npar_mask_img.any():
                ks = np.where(npar_mask_img)[0]
                M_quad_img = _neumann_quadrature_batch(
                    p1_a, da_a, la_a, ua_a,
                    p1_b_img[ks], db_b_img[ks], la_b_img[ks], ua_b_img[ks],
                    dots_img[ks],
                )
                M_img[ks] = M_quad_img
            M_pri = M_pri + M_img

        L[i, slc] = M_pri
        L[slc, i] = M_pri

    return L

neumann_mutual

neumann_mutual(
    p1_a: np.ndarray,
    p2_a: np.ndarray,
    p1_b: np.ndarray,
    p2_b: np.ndarray,
    *,
    min_distance: float = 1e-06,
    parallel_tol: float = 1e-09
) -> float

Neumann mutual-inductance integral between two straight segments.

Hybrid implementation:

1. If the two segments are parallel (or anti-parallel) within parallel_tol, the closed-form Grover expression in :func:_parallel_filaments_mutual is used — exact, no quadrature error. This fast path covers the bulk of the typical inductance assembly (parallel PEN strands, consecutive sub-segments along the same conductor).
2. Otherwise the Neumann double integral is evaluated by 16×16 Gauss–Legendre quadrature. Empirically accurate to ≲ 0.05 % across the typical geometry range; see ADR-0004 for the calibration data.

The kernel \(1/r\) is clamped at min_distance to suppress the integrable singularity for segments that touch — physically the diagonal uses :func:thin_wire_self_inductance instead, so this only affects pathological inputs.

Parameters:

Name	Type	Description	Default
p1_a	ndarray	Endpoints of segment a as (x, y, z) arrays in metres.	required
p2_a	ndarray	Endpoints of segment a as (x, y, z) arrays in metres.	required
p1_b	ndarray	Endpoints of segment b.	required
p2_b	ndarray	Endpoints of segment b.	required
min_distance	float	Numerical floor on \(|r_a - r_b|\) in metres.	1e-06
parallel_tol	float	Tolerance on \(||\hat{u}_a\cdot\hat{u}_b| - 1|\) that triggers the closed-form fast path.	1e-09

Returns:

Name	Type	Description
M	float	Partial mutual inductance in H.

Source code in src/groundfield/coupling/inductance.py

def neumann_mutual(
    p1_a: np.ndarray, p2_a: np.ndarray,
    p1_b: np.ndarray, p2_b: np.ndarray,
    *,
    min_distance: float = 1e-6,
    parallel_tol: float = 1e-9,
) -> float:
    """Neumann mutual-inductance integral between two straight segments.

    Hybrid implementation:

    1. If the two segments are **parallel** (or anti-parallel) within
       ``parallel_tol``, the closed-form Grover expression in
       :func:`_parallel_filaments_mutual` is used — exact, no
       quadrature error. This fast path covers the bulk of the typical
       inductance assembly (parallel PEN strands, consecutive
       sub-segments along the same conductor).
    2. Otherwise the Neumann double integral is evaluated by 16×16
       Gauss–Legendre quadrature. Empirically accurate to ≲ 0.05 %
       across the typical geometry range; see ADR-0004 for the
       calibration data.

    The kernel $1/r$ is clamped at ``min_distance`` to suppress the
    integrable singularity for segments that touch — physically the
    diagonal uses :func:`thin_wire_self_inductance` instead, so
    this only affects pathological inputs.

    Parameters
    ----------
    p1_a, p2_a
        Endpoints of segment *a* as ``(x, y, z)`` arrays in metres.
    p1_b, p2_b
        Endpoints of segment *b*.
    min_distance
        Numerical floor on $|r_a - r_b|$ in metres.
    parallel_tol
        Tolerance on $||\\hat{u}_a\\cdot\\hat{u}_b| - 1|$ that
        triggers the closed-form fast path.

    Returns
    -------
    M : float
        Partial mutual inductance in H.
    """
    p1_a = np.asarray(p1_a, dtype=float)
    p2_a = np.asarray(p2_a, dtype=float)
    p1_b = np.asarray(p1_b, dtype=float)
    p2_b = np.asarray(p2_b, dtype=float)

    # Fast path: parallel filaments with closed form.
    M_par = _parallel_filaments_mutual(
        p1_a, p2_a, p1_b, p2_b, parallel_tol=parallel_tol,
    )
    if M_par is not None:
        return M_par

    # Generic 3-D case: Gauss–Legendre quadrature.
    da = p2_a - p1_a
    db = p2_b - p1_b
    la = float(np.linalg.norm(da))
    lb = float(np.linalg.norm(db))
    if la <= 0.0 or lb <= 0.0:
        raise ValueError("segment lengths must be positive")
    ua = da / la
    ub = db / lb
    dot = float(ua @ ub)
    if dot == 0.0:
        return 0.0  # perpendicular — orthogonality kills the integral
    s_nodes = 0.5 * (_GL_NODES + 1.0)
    w_nodes = 0.5 * _GL_WEIGHTS
    pts_a = p1_a[None, :] + s_nodes[:, None] * da[None, :]
    pts_b = p1_b[None, :] + s_nodes[:, None] * db[None, :]
    diff = pts_a[:, None, :] - pts_b[None, :, :]
    r = np.linalg.norm(diff, axis=2)
    np.maximum(r, min_distance, out=r)
    inv_r = 1.0 / r
    integral = float(np.einsum("i,j,ij->", w_nodes, w_nodes, inv_r))
    return MU_0 / (4.0 * math.pi) * dot * la * lb * integral

parallel_segments_mutual

parallel_segments_mutual(
    length: float, distance: float
) -> float

Closed-form mutual inductance of two parallel coaxial segments.

For two segments of equal length \(\ell\) that share the same axis direction and are placed at perpendicular distance \(d\),

\[ M_\parallel(\ell, d) \;=\; \frac{\mu_0\,\ell}{2\pi} \Bigl[\ln\!\Bigl(\frac{\ell + \sqrt{\ell^2 + d^2}}{d}\Bigr) - \frac{\sqrt{\ell^2 + d^2} - d}{\ell}\Bigr]. \]

Used both as an internal closed-form fast path (when the segment pair geometry matches the assumption) and as a reference in the test suite.

Parameters:

Name	Type	Description	Default
length	float	Segment length \(\ell\) in m (both segments have the same).	required
distance	float	Perpendicular distance \(d\) between the two parallel axes in m. Must be > 0.	required

Returns:

Name	Type	Description
M	float	Partial mutual inductance in H.

Source code in src/groundfield/coupling/inductance.py

def parallel_segments_mutual(length: float, distance: float) -> float:
    """Closed-form mutual inductance of two parallel coaxial segments.

    For two segments of equal length $\\ell$ that share the same
    axis direction and are placed at perpendicular distance $d$,

    $$
    M_\\parallel(\\ell, d) \\;=\\; \\frac{\\mu_0\\,\\ell}{2\\pi}
    \\Bigl[\\ln\\!\\Bigl(\\frac{\\ell + \\sqrt{\\ell^2 + d^2}}{d}\\Bigr)
          - \\frac{\\sqrt{\\ell^2 + d^2} - d}{\\ell}\\Bigr].
    $$

    Used both as an internal closed-form fast path (when the segment
    pair geometry matches the assumption) and as a reference in the
    test suite.

    Parameters
    ----------
    length
        Segment length $\\ell$ in m (both segments have the same).
    distance
        Perpendicular distance $d$ between the two parallel axes
        in m. Must be > 0.

    Returns
    -------
    M : float
        Partial mutual inductance in H.
    """
    if length <= 0.0 or distance <= 0.0:
        raise ValueError("length and distance must be positive")
    s = math.hypot(length, distance)
    return (MU_0 * length) / (2.0 * math.pi) * (
        math.log((length + s) / distance) - (s - distance) / length
    )

perfect_mirror_self_pair_inductance

perfect_mirror_self_pair_inductance(
    p1: np.ndarray, p2: np.ndarray, wire_radius: float
) -> float

Self-inductance contribution of a segment plus its earth image.

For a single segment the partial self-inductance against itself uses the thin-wire formula. The image segment (mirrored at the soil surface) is treated as an external segment whose Neumann integral with the original is added to the total.

Parameters:

Name	Type	Description	Default
p1	ndarray	Endpoints in metres (with \(z > 0\) pointing into the soil).	required
p2	ndarray	Endpoints in metres (with \(z > 0\) pointing into the soil).	required
wire_radius	float	Wire radius in m (only used for the thin-wire self term).	required

Returns:

Name	Type	Description
L	float	Self-inductance plus image contribution, in henries.

Source code in src/groundfield/coupling/inductance.py

def perfect_mirror_self_pair_inductance(
    p1: np.ndarray, p2: np.ndarray, wire_radius: float
) -> float:
    """Self-inductance contribution of a segment plus its earth image.

    For a single segment the *partial* self-inductance against itself
    uses the thin-wire formula. The image segment (mirrored at the
    soil surface) is treated as an external segment whose Neumann
    integral with the original is added to the total.

    Parameters
    ----------
    p1, p2
        Endpoints in metres (with $z > 0$ pointing into the soil).
    wire_radius
        Wire radius in m (only used for the thin-wire self term).

    Returns
    -------
    L : float
        Self-inductance plus image contribution, in henries.
    """
    p1 = np.asarray(p1, dtype=float)
    p2 = np.asarray(p2, dtype=float)
    L_self = thin_wire_self_inductance(
        float(np.linalg.norm(p2 - p1)), wire_radius,
    )
    p1_img = _mirror(p1)
    p2_img = _mirror(p2)
    # Image segment runs in the opposite direction (current image of a
    # vertical filament is anti-parallel; for a horizontal filament it
    # is parallel). The general rule is: image of (p1 → p2) is
    # (p1' → p2'), with the *same* tangent direction up to the
    # sign of ẑ. For Neumann the dot product ǔ_a·ǔ_b takes care of
    # the sign automatically — we just feed the mirrored endpoints.
    L_image = neumann_mutual(p1, p2, p1_img, p2_img)
    return L_self + L_image

thin_wire_self_inductance

thin_wire_self_inductance(
    length: float,
    wire_radius: float,
    *,
    include_internal: bool = True
) -> float

Closed-form self-inductance of a straight thin-wire segment.

Decomposed into the external Grover 1946 thin-wire term and the internal contribution from the magnetic field inside the conductor (uniform DC current distribution, non-magnetic material):

\[ L_\text{ext} \;=\; \frac{\mu_0\,\ell}{2\pi} \Bigl[\ln\!\Bigl(\frac{2\ell}{a}\Bigr) - 1\Bigr], \qquad L_\text{int} \;=\; \frac{\mu_0\,\ell}{8\pi}. \]

Adding the internal term reproduces the standard Oeding/Oswald loop-inductance formula (Gl. 9.13c) when used as the diagonal contribution of a two-wire loop. For high-frequency studies where the current flows on the surface only (skin effect at \(f \gg f_\text{skin}\)), pass include_internal=False and the external term alone is returned.

Parameters:

Name	Type	Description	Default
length	float	Segment length \(\ell\) in metres.	required
wire_radius	float	Wire radius \(a\) in metres.	required
include_internal	bool	True (default) — add the internal-field contribution \(\mu_0 \ell / (8\pi)\). Appropriate for the \(f < 1\,\mathrm{kHz}\) regime with mostly uniform current distribution. Set False for the pure external-field thin-wire formula.	True

Returns:

Name	Type	Description
L	float	Partial self-inductance in henries.

Source code in src/groundfield/coupling/inductance.py

def thin_wire_self_inductance(
    length: float,
    wire_radius: float,
    *,
    include_internal: bool = True,
) -> float:
    """Closed-form self-inductance of a straight thin-wire segment.

    Decomposed into the **external** Grover 1946 thin-wire term and
    the **internal** contribution from the magnetic field inside the
    conductor (uniform DC current distribution, non-magnetic
    material):

    $$
    L_\\text{ext} \\;=\\; \\frac{\\mu_0\\,\\ell}{2\\pi}
    \\Bigl[\\ln\\!\\Bigl(\\frac{2\\ell}{a}\\Bigr) - 1\\Bigr],
    \\qquad
    L_\\text{int} \\;=\\; \\frac{\\mu_0\\,\\ell}{8\\pi}.
    $$

    Adding the internal term reproduces the standard Oeding/Oswald
    loop-inductance formula (Gl. 9.13c) when used as the diagonal
    contribution of a two-wire loop. For high-frequency studies
    where the current flows on the surface only (skin effect at
    $f \\gg f_\\text{skin}$), pass ``include_internal=False`` and
    the external term alone is returned.

    Parameters
    ----------
    length
        Segment length $\\ell$ in metres.
    wire_radius
        Wire radius $a$ in metres.
    include_internal
        ``True`` (default) — add the internal-field contribution
        $\\mu_0 \\ell / (8\\pi)$. Appropriate for the
        $f < 1\\,\\mathrm{kHz}$ regime with mostly uniform
        current distribution. Set ``False`` for the pure
        external-field thin-wire formula.

    Returns
    -------
    L : float
        Partial self-inductance in henries.
    """
    if length <= 0.0 or wire_radius <= 0.0:
        raise ValueError("length and wire_radius must be positive")
    L_ext = (MU_0 * length) / (2.0 * math.pi) * (
        math.log(2.0 * length / wire_radius) - 1.0
    )
    if not include_internal:
        return L_ext
    L_int = (MU_0 * length) / (8.0 * math.pi)
    return L_ext + L_int

Carson earth-return correction (ADR-0005)

Adds Carson 1926's finite-conductivity correction \(\Delta Z_\text{Carson}(\omega)\) on top of the perfect-mirror inductance. The branch-impedance block becomes

\[ Z_b(\omega) \;=\; R \;+\; j\omega\, L_\text{Neumann} \;+\; \Delta Z_\text{Carson}(\omega, \sigma_\text{earth}, h_i, h_j, d_{ij}). \]

The correction is evaluated as

\[ \Delta Z_\text{Carson}(\omega) \;=\; \frac{\omega\,\mu_0}{\pi}\, \bigl[P(a, \theta) \,+\, j\,Q(a, \theta)\bigr], \qquad a \;=\; D\,\sqrt{\omega\,\mu_0\,\sigma_\text{earth}} \;=\; \frac{D\sqrt{2}}{\delta(\omega)}, \]

with \(D = 2h_i\) (\(\theta = 0\)) for the self contribution and \(D = \sqrt{(h_i+h_j)^2 + d_{ij}^2}\), \(\theta = \arctan(d_{ij}/(h_i+h_j))\) for the mutual contribution. \(\delta(\omega) = \sqrt{2 / (\omega\mu_0\sigma_\text{earth})}\) is the electromagnetic skin depth in soil — the natural length scale at which the perfect-mirror approximation starts to break down.

Three regimes

Following Carson 1926 §III the implementation switches between three numerical regimes:

Regime	Range of \(a\)	Method
Small	\(a \le 0.25\)	Closed-form leading-term expansion (Carson eqs. 34/35)
Intermediate	\(0.25 < a \le 5\)	64-point Gauss–Legendre quadrature of Carson eq. 29
Asymptotic	\(a > 5\)	Inverse-power expansion (Carson eqs. 36/37)

Validity and limitations

• Homogeneous soil — the Carson series is exact (within the quasi-static / sub-kHz assumption Carson himself states) when \(\sigma_\text{earth}\) is uniform.
• Layered soil — the implementation falls back to \(\sigma = 1/\rho_1\) of the upper layer with a runtime UserWarning. For a rigorous result switch to backend="mom_sommerfeld", which uses the full Pollaczek kernel.
• Frequency — derived for \(\omega \ll 1 / \mu_0\sigma\); for typical (sub-kHz, \(\rho_\text{earth} \in [50, 5000]\,\Omega\,\mathrm{m}\)) the assumption is comfortably satisfied.
• Geometry — Carson's derivation assumes parallel wires above a plane homogeneous half-space. Non-parallel segment pairs are handled by projection onto the parallel component (orthogonal components contribute zero by Neumann symmetry).

carson

Carson 1926 correction for the earth-return path.

This module implements the magnetic-image correction described in ADR-0005. It complements the perfect-mirror Neumann assembly from ADR-0004 (in :mod:groundfield.coupling.inductance) with the finite-conductivity contribution that Carson 1926 derived for a homogeneous, semi-infinite, conductive half-space.

Mathematical background

Carson 1926 (Bell Syst. Tech. J. 5(4)) writes the earth-return correction to the per-unit-length series impedance of a long straight wire above a homogeneous earth as

\[ \Delta Z_\text{Carson}(\omega) \;=\; \frac{\omega\,\mu_0}{\pi}\, \bigl[P(a, \theta) \,+\, j\,Q(a, \theta)\bigr], \]

with the dimensionless Carson parameter

\[ a \;=\; D\,\sqrt{\omega\,\mu_0\,\sigma_\text{earth}} \;=\; \frac{D\sqrt{2}}{\delta(\omega)}, \qquad \delta(\omega) \;=\; \sqrt{2 \,/\, (\omega\mu_0\sigma_\text{earth})} \;[\text{skin depth in soil}], \]

where \(D = 2h_i\) for the self-impedance correction (\(\theta = 0\)) and \(D = \sqrt{(h_i + h_j)^2 + d_{ij}^2}\), \(\theta = \arctan(d_{ij} / (h_i + h_j))\) for the mutual-impedance correction between two parallel wires at heights \(h_i, h_j\) with horizontal separation \(d_{ij}\).

The functions \(P, Q\) are the real and imaginary parts of Carson's infinite integral \(J(p, q)\) (Carson eq. 29, with \(p = a\cos\theta\), \(q = a\sin\theta\)):

\[ J(p, q) \;=\; \int_0^{\infty}\!\bigl(\sqrt{\mu^2 + j}-\mu\bigr)\, e^{-p\mu}\cos(q\mu)\,d\mu \;=\; P(a,\theta) + j\,Q(a,\theta). \]

Three evaluation regimes

We follow Carson's own discussion in section III of the original paper:

1. Small \(a\) (\(a \le 0.25\)) — Carson eqs. 34/35, leading-term form. Closed form, only \(\sin / \cos / \ln\).
2. Intermediate \(a\) (\(0.25 < a \le 5\)) — direct numerical quadrature of Carson's \(J(p, q)\) via Gauss–Legendre on a truncated interval. The original Carson 1926 series in that range is technically convergent but its recurrence is numerically delicate; quadrature is robust and converges to machine precision in \(\le 64\) nodes.
3. Large \(a\) (\(a > 5\)) — Carson eqs. 36/37, asymptotic expansion in inverse powers of \(a\).

Every regime boundary is smoke-tested for continuity at \(\le 10^{-6}\) — see :mod:tests.test_carson_coupling.

Implementation notes

• All formulas use SI units throughout. Carson's CGS pre-factor \(4\omega\) becomes \(\omega\mu_0 / \pi\) in SI.
• \(\sigma_\text{earth}\) is an explicit argument; the caller is responsible for selecting the right resistivity (homogeneous \(1/\rho\), or upper-layer \(1/\rho_1\) with a warning).
• \(\omega = 0\) or \(\sigma = 0\) short-circuit to 0+0j: the prefactor \(\omega\mu_0/\pi\) vanishes, so the (logarithmically diverging) \(Q\)-asymptote is harmless.
• The complex-depth Deri/Semlyen approximation is provided as an internal sanity check (:func:deri_semlyen_correction); it is not the production code path.

References

• Carson, J. R. (1926). Wave propagation in overhead wires with ground return. Bell Syst. Tech. J. 5(4), 539–554.
• Deri, A.; Tevan, G.; Semlyen, A.; Castanheira, A. (1981). The complex ground return plane. IEEE Trans. PAS 100(8), 3686–3693.
• Tleis, N. D. (2008). Power Systems Modelling and Fault Analysis, Newnes, ch. 3.

carson_mutual_correction

carson_mutual_correction(
    omega: float,
    height_i: float,
    height_j: float,
    horizontal_distance: float,
    sigma_earth: float,
) -> complex

Carson earth-return correction between two parallel horizontal wires.

Evaluates the per-unit-length mutual-impedance correction

\[ \Delta Z_\text{mutual}(\omega) \;=\; \frac{\omega\mu_0}{\pi} \bigl[P(a_m, \theta_m) \,+\, j Q(a_m, \theta_m)\bigr], \]

with

\[ a_m \;=\; D_m \sqrt{\omega\mu_0\sigma_\text{earth}}, \quad D_m \;=\; \sqrt{(h_i + h_j)^2 + d^2}, \quad \theta_m \;=\; \arctan(d / (h_i + h_j)), \]

after Carson eq. 31. The parallel-wire assumption is the same as for ADR-0004's Neumann fast path; for non-parallel segments the caller must split into projection components.

Parameters:

Name	Type	Description	Default
omega	float	Angular frequency in rad/s.	required
height_i	float	Heights above earth surface in metres (both positive).	required
height_j	float	Heights above earth surface in metres (both positive).	required
horizontal_distance	float	Horizontal separation \(d\) between the two wires in metres.	required
sigma_earth	float	Earth conductivity in S/m.	required

Returns:

Name	Type	Description
Z	complex	Per-unit-length mutual correction in \(\Omega/\text{m}\).

Source code in src/groundfield/coupling/carson.py

def carson_mutual_correction(
    omega: float,
    height_i: float,
    height_j: float,
    horizontal_distance: float,
    sigma_earth: float,
) -> complex:
    """Carson earth-return correction between two parallel horizontal wires.

    Evaluates the per-unit-length mutual-impedance correction

    $$
    \\Delta Z_\\text{mutual}(\\omega) \\;=\\; \\frac{\\omega\\mu_0}{\\pi}
    \\bigl[P(a_m, \\theta_m) \\,+\\, j Q(a_m, \\theta_m)\\bigr],
    $$

    with

    $$
    a_m \\;=\\; D_m \\sqrt{\\omega\\mu_0\\sigma_\\text{earth}}, \\quad
    D_m \\;=\\; \\sqrt{(h_i + h_j)^2 + d^2}, \\quad
    \\theta_m \\;=\\; \\arctan(d / (h_i + h_j)),
    $$

    after Carson eq. 31. The parallel-wire assumption is the same
    as for ADR-0004's Neumann fast path; for non-parallel segments
    the caller must split into projection components.

    Parameters
    ----------
    omega
        Angular frequency in rad/s.
    height_i, height_j
        Heights above earth surface in metres (both positive).
    horizontal_distance
        Horizontal separation $d$ between the two wires in metres.
    sigma_earth
        Earth conductivity in S/m.

    Returns
    -------
    Z : complex
        Per-unit-length mutual correction in $\\Omega/\\text{m}$.
    """
    if height_i <= 0.0 or height_j <= 0.0:
        raise ValueError("heights must be positive")
    if horizontal_distance < 0.0:
        raise ValueError("horizontal_distance must be non-negative")
    if omega <= 0.0 or sigma_earth <= 0.0:
        return 0.0 + 0.0j
    h_sum = height_i + height_j
    D = math.hypot(h_sum, horizontal_distance)
    theta = math.atan2(horizontal_distance, h_sum)
    a = carson_parameter(D, omega, sigma_earth)
    P, Q = carson_p_q(a, theta)
    pref = omega * MU_0 / math.pi
    return complex(pref * P, pref * Q)

carson_p_q

carson_p_q(a: float, theta: float) -> tuple[float, float]

Evaluate Carson's \(P(a, \theta)\), \(Q(a, \theta)\).

Dispatches to the appropriate regime based on the magnitude of \(a\):

• \(a \le 0.25\): closed-form leading-term expansion (Carson eqs. 34/35).
• \(0.25 < a \le 5\): direct numerical quadrature of Carson's \(J(p, q)\) (Gauss–Legendre for \(\theta \approx 0\), QUADPACK Fourier integrator otherwise — audit 2026-07-08, WP-B2).
• \(a > 5\): asymptotic expansion (Carson eqs. 36/37), except for \(\theta > 1.4\) rad where the expansion degrades and the quadrature path is used instead.

The regimes are continuous at the boundaries to within the tolerances documented in ADR-0005 §5/§6.

Parameters:

Name	Type	Description	Default
a	float	Dimensionless Carson parameter \(a = D \sqrt{\omega \mu_0 \sigma_\text{earth}}\). Must be \(\ge 0\).	required
theta	float	Angle of the image-distance vector to the vertical, in radians. \(\theta = 0\) for the self-impedance correction, \(\theta = \arctan(d/(h_i+h_j))\) for the mutual.	required

Returns:

Type	Description
P, Q : tuple[float, float]	Real and imaginary parts of Carson's integral \(J\).

Raises:

Type	Description
ValueError	If a is negative.

Source code in src/groundfield/coupling/carson.py

def carson_p_q(a: float, theta: float) -> tuple[float, float]:
    """Evaluate Carson's $P(a, \\theta)$, $Q(a, \\theta)$.

    Dispatches to the appropriate regime based on the magnitude of
    $a$:

    - $a \\le 0.25$: closed-form leading-term expansion (Carson
      eqs. 34/35).
    - $0.25 < a \\le 5$: direct numerical quadrature of Carson's
      $J(p, q)$ (Gauss–Legendre for $\\theta \\approx 0$, QUADPACK
      Fourier integrator otherwise — audit 2026-07-08, WP-B2).
    - $a > 5$: asymptotic expansion (Carson eqs. 36/37), except for
      $\\theta > 1.4$ rad where the expansion degrades and the
      quadrature path is used instead.

    The regimes are continuous at the boundaries to within the
    tolerances documented in ADR-0005 §5/§6.

    Parameters
    ----------
    a : float
        Dimensionless Carson parameter $a = D \\sqrt{\\omega
        \\mu_0 \\sigma_\\text{earth}}$. Must be $\\ge 0$.
    theta : float
        Angle of the image-distance vector to the vertical, in
        radians. $\\theta = 0$ for the self-impedance correction,
        $\\theta = \\arctan(d/(h_i+h_j))$ for the mutual.

    Returns
    -------
    P, Q : tuple[float, float]
        Real and imaginary parts of Carson's integral $J$.

    Raises
    ------
    ValueError
        If ``a`` is negative.
    """
    if a < 0.0:
        raise ValueError(f"Carson parameter a must be non-negative, got {a}")
    if a == 0.0:
        return math.pi / 8.0, math.inf
    if a <= _REGIME_SMALL_MAX:
        return _p_q_small(a, theta)
    if a <= _REGIME_LARGE_MIN or theta > _ASYMPTOTE_THETA_MAX:
        return _p_q_quadrature(a, theta)
    return _p_q_large(a, theta)

carson_parameter

carson_parameter(
    distance: float, omega: float, sigma_earth: float
) -> float

Dimensionless Carson parameter \(a = D \sqrt{\omega\mu_0\sigma}\).

Parameters:

Name	Type	Description	Default
distance	float	Geometric distance \(D\) in metres (\(2h\) for self-impedance, \(\sqrt{(h_i+h_j)^2 + d^2}\) for mutual).	required
omega	float	Angular frequency in rad/s.	required
sigma_earth	float	Earth conductivity in S/m.	required

Returns:

Name	Type	Description
a	float	Dimensionless Carson parameter \(a\).

Source code in src/groundfield/coupling/carson.py

def carson_parameter(distance: float, omega: float, sigma_earth: float) -> float:
    """Dimensionless Carson parameter $a = D \\sqrt{\\omega\\mu_0\\sigma}$.

    Parameters
    ----------
    distance
        Geometric distance $D$ in metres
        ($2h$ for self-impedance, $\\sqrt{(h_i+h_j)^2 + d^2}$
        for mutual).
    omega
        Angular frequency in rad/s.
    sigma_earth
        Earth conductivity in S/m.

    Returns
    -------
    a : float
        Dimensionless Carson parameter $a$.
    """
    if distance < 0.0:
        raise ValueError("distance must be non-negative")
    if omega <= 0.0 or sigma_earth <= 0.0:
        return 0.0
    return distance * math.sqrt(omega * MU_0 * sigma_earth)

carson_self_correction

carson_self_correction(
    omega: float, height: float, sigma_earth: float
) -> complex

Carson earth-return correction for a single horizontal wire.

Evaluates the per-unit-length impedance correction

\[ \Delta Z_\text{self}(\omega) \;=\; \frac{\omega \mu_0}{\pi} \bigl[P(a_s, 0) + j Q(a_s, 0)\bigr], \]

with \(a_s = 2h\sqrt{\omega\mu_0\sigma_\text{earth}}\) and \(\theta = 0\) (Carson eq. 30, ADR-0005).

Parameters:

Name	Type	Description	Default
omega	float	Angular frequency \(\omega = 2\pi f\) in rad/s.	required
height	float	Height of the wire above the earth surface in metres (positive — for buried wires use the Sunde-equivalent height).	required
sigma_earth	float	Earth conductivity in S/m (\(\sigma = 1/\rho\) for a homogeneous earth).	required

Returns:

Name	Type	Description
Z	complex	Per-unit-length earth-return correction in \(\Omega/\text{m}\).

Raises:

Type	Description
ValueError	If height is non-positive.

Source code in src/groundfield/coupling/carson.py

def carson_self_correction(
    omega: float,
    height: float,
    sigma_earth: float,
) -> complex:
    """Carson earth-return correction for a single horizontal wire.

    Evaluates the per-unit-length impedance correction

    $$
    \\Delta Z_\\text{self}(\\omega) \\;=\\; \\frac{\\omega \\mu_0}{\\pi}
    \\bigl[P(a_s, 0) + j Q(a_s, 0)\\bigr],
    $$

    with $a_s = 2h\\sqrt{\\omega\\mu_0\\sigma_\\text{earth}}$ and
    $\\theta = 0$ (Carson eq. 30, ADR-0005).

    Parameters
    ----------
    omega
        Angular frequency $\\omega = 2\\pi f$ in rad/s.
    height
        Height of the wire above the earth surface in metres
        (positive — for buried wires use the Sunde-equivalent
        height).
    sigma_earth
        Earth conductivity in S/m
        ($\\sigma = 1/\\rho$ for a homogeneous earth).

    Returns
    -------
    Z : complex
        Per-unit-length earth-return correction in $\\Omega/\\text{m}$.

    Raises
    ------
    ValueError
        If ``height`` is non-positive.
    """
    if height <= 0.0:
        raise ValueError("height must be positive (use |z| of the wire)")
    if omega <= 0.0 or sigma_earth <= 0.0:
        return 0.0 + 0.0j
    a = carson_parameter(2.0 * height, omega, sigma_earth)
    P, Q = carson_p_q(a, 0.0)
    pref = omega * MU_0 / math.pi
    return complex(pref * P, pref * Q)

deri_semlyen_correction

deri_semlyen_correction(
    omega: float,
    height_i: float,
    height_j: float,
    horizontal_distance: float,
    sigma_earth: float,
) -> complex

Deri/Semlyen 1981 complex-depth approximation.

Replaces the Carson integral by the closed-form expression

\[ \Delta Z_\text{Deri-Semlyen} \;=\; \frac{j\omega\mu_0}{2\pi} \ln\!\Bigl(\frac{D'}{D}\Bigr), \]

with \(D = \sqrt{(h_i - h_j)^2 + d^2}\) the direct distance, \(D' = \sqrt{(h_i + h_j + 2p)^2 + d^2}\) the distance to a complex-depth image, and the complex penetration depth

\[ p \;=\; 1 \big/ \sqrt{j\omega\mu_0\sigma_\text{earth}}. \]

The Deri/Semlyen approximation is not the production code path. It is provided as an alternative independent estimator that the test suite can compare against the Carson series — agreement within \(\approx 5\,\%\) over the typical parameter range confirms that neither implementation contains a sign or pre-factor bug.

Parameters:

Name	Type	Description	Default
omega	float	Same meaning as in :func:carson_mutual_correction.	required
height_i	float	Same meaning as in :func:carson_mutual_correction.	required
height_j	float	Same meaning as in :func:carson_mutual_correction.	required
horizontal_distance	float	Same meaning as in :func:carson_mutual_correction.	required
sigma_earth	float	Same meaning as in :func:carson_mutual_correction.	required

Returns:

Name	Type	Description
Z	complex	Per-unit-length earth-return correction in \(\Omega/\text{m}\) according to Deri/Semlyen 1981.

Notes

The "self" version is recovered by setting height_i = height_j and horizontal_distance = 0: the direct-distance term degenerates to the wire radius, which the caller must substitute manually (the formula does not include a wire-radius regularisation by itself).

Source code in src/groundfield/coupling/carson.py

def deri_semlyen_correction(
    omega: float,
    height_i: float,
    height_j: float,
    horizontal_distance: float,
    sigma_earth: float,
) -> complex:
    """Deri/Semlyen 1981 complex-depth approximation.

    Replaces the Carson integral by the closed-form expression

    $$
    \\Delta Z_\\text{Deri-Semlyen} \\;=\\; \\frac{j\\omega\\mu_0}{2\\pi}
    \\ln\\!\\Bigl(\\frac{D'}{D}\\Bigr),
    $$

    with $D = \\sqrt{(h_i - h_j)^2 + d^2}$ the direct distance,
    $D' = \\sqrt{(h_i + h_j + 2p)^2 + d^2}$ the distance to a
    complex-depth image, and the complex penetration depth

    $$
    p \\;=\\; 1 \\big/ \\sqrt{j\\omega\\mu_0\\sigma_\\text{earth}}.
    $$

    The Deri/Semlyen approximation is **not** the production
    code path. It is provided as an alternative independent
    estimator that the test suite can compare against the Carson
    series — agreement within $\\approx 5\\,\\%$ over the typical
    parameter range confirms that neither implementation contains
    a sign or pre-factor bug.

    Parameters
    ----------
    omega, height_i, height_j, horizontal_distance, sigma_earth
        Same meaning as in :func:`carson_mutual_correction`.

    Returns
    -------
    Z : complex
        Per-unit-length earth-return correction in $\\Omega/\\text{m}$
        according to Deri/Semlyen 1981.

    Notes
    -----
    The "self" version is recovered by setting
    ``height_i = height_j`` and ``horizontal_distance = 0``: the
    direct-distance term degenerates to the wire radius, which the
    caller must substitute manually (the formula does not include
    a wire-radius regularisation by itself).
    """
    if omega <= 0.0 or sigma_earth <= 0.0:
        return 0.0 + 0.0j
    p = 1.0 / np.sqrt(1j * omega * MU_0 * sigma_earth)
    D = math.hypot(height_i - height_j, horizontal_distance)
    if D == 0.0:
        D = 1e-9
    D_prime = np.sqrt((height_i + height_j + 2.0 * p) ** 2 + horizontal_distance ** 2)
    return complex(1j * omega * MU_0 / (2.0 * math.pi) * np.log(D_prime / D))

skin_depth

skin_depth(omega: float, sigma_earth: float) -> float

Electromagnetic skin depth in soil.

Returns \(\delta = \sqrt{2 / (\omega \mu_0 \sigma)}\) in metres. Diverges at \(\omega = 0\) (purely conductive earth- return path); the caller must handle that case by skipping the Carson correction altogether.

Parameters:

Name	Type	Description	Default
omega	float	Angular frequency \(\omega = 2\pi f\) in rad/s.	required
sigma_earth	float	Earth conductivity \(\sigma\) in S/m.	required

Returns:

Name	Type	Description
delta	float	Skin depth in metres. +inf when omega or sigma_earth is zero.

Source code in src/groundfield/coupling/carson.py

def skin_depth(omega: float, sigma_earth: float) -> float:
    """Electromagnetic skin depth in soil.

    Returns $\\delta = \\sqrt{2 / (\\omega \\mu_0 \\sigma)}$ in
    metres. Diverges at $\\omega = 0$ (purely conductive earth-
    return path); the caller must handle that case by skipping the
    Carson correction altogether.

    Parameters
    ----------
    omega
        Angular frequency $\\omega = 2\\pi f$ in rad/s.
    sigma_earth
        Earth conductivity $\\sigma$ in S/m.

    Returns
    -------
    delta : float
        Skin depth in metres. ``+inf`` when ``omega`` or
        ``sigma_earth`` is zero.
    """
    if omega <= 0.0 or sigma_earth <= 0.0:
        return math.inf
    return math.sqrt(2.0 / (omega * MU_0 * sigma_earth))

Sommerfeld geometric earth-return Green's function (ADR-0006)

The rigorous formulation of the earth-return inductive coupling. Integrates the σ-dependent vector-potential Green's function

\[ G_\text{mag}(\vec{r}, \vec{r}';\,\omega,\sigma_e) \;=\; \frac{1}{R} \;+\; \int_0^{\infty}\! \Gamma_\text{mag}(\lambda)\, e^{-\lambda(z+z')}\,J_0(\lambda\rho)\,d\lambda \]

over the actual segment-pair geometry, where the reflection coefficient \(\Gamma_\text{mag}\) encodes the homogeneous or layered-earth structure. Reduces to the perfect-mirror case at \(\sigma\to\infty\), to free space at \(\sigma\to 0\), and converges to ADR-0005's per-meter Carson asymptote for long parallel wires over homogeneous earth. Unlike Carson, it correctly handles short wires, non-parallel geometries, and layered earth without approximation (Pollaczek/Wait kernel).

When to use which

Problem class	Recommended option
Reference computations, layered earth, mixed wire lengths	"sommerfeld"
Long parallel PEN over homogeneous earth, fast scoping	"carson_series"
Pure DC studies, perfect-mirror reference	"perfect_mirror"

sommerfeld_inductance

Geometric Sommerfeld earth-return Green function (ADR-0006).

This module implements the rigorous geometric formulation of the inductive earth-return coupling described in ADR-0006. Compared to ADR-0005's Carson per-meter scaling, it integrates the actual vector-potential Green's function over the segment-pair geometry, which correctly handles short wires, non-parallel arrangements, and layered earth.

Mathematical background

(Rewritten in the 2026-07-08 audit, WP-C1: the historic kernel 1/R + ∫Γ·e^{-λ(z+z')}J0 dropped the in-soil propagation of the primary term — precisely where the buried-wire earth-return resistance lives — and applied the soil-side reflection sign to overhead conductors. The implemented kernel is now the genuine quasi-static Pollaczek form, dispatched by which side of the soil surface the two segments are on.)

For a horizontal current element in a conducting half-space (\(z > 0\) soil, \(z < 0\) air, \(\gamma^2 = j\omega\mu_0\sigma_1\), \(u_1 = \sqrt{\lambda^2 + \gamma^2}\)), the quasi-static vector-potential Green's functions are:

Buried source, buried observer (depths \(h, h' > 0\)):

\[ G_\text{bb} \;=\; \frac{e^{-\gamma R}}{R} \;+\; \int_0^{\infty}\! \frac{\lambda}{u_1}\,\Gamma(\lambda)\, e^{-u_1 (h + h')}\,J_0(\lambda\rho)\,d\lambda, \]

Overhead source and observer (heights \(H, H' > 0\)):

\[ G_\text{oo} \;=\; \frac{1}{R} \;-\; \int_0^{\infty}\! \Gamma(\lambda)\, e^{-\lambda (H + H')}\,J_0(\lambda\rho)\,d\lambda \qquad\text{(Carson's geometry)}, \]

Mixed (buried \(h\), overhead \(H\)):

\[ G_\text{bo} \;=\; \int_0^{\infty}\! \frac{2\lambda}{u_1 + \lambda}\, e^{-u_1 h - \lambda H}\,J_0(\lambda\rho)\,d\lambda, \]

with the reflection coefficient

• homogeneous earth (Pillar A): \(\Gamma(\lambda) = (u_1 - \lambda)/(u_1 + \lambda)\),
• \(n\)-layer earth (Pillar B): the recursive Wait-style reflection coefficient (Wait 1972 §3, Tleis 2008 §3.5); the \(\lambda/u_1\) weight and the exponents use the top-layer \(u_1\) (the conductors live in layer 1).

The functions in this module return the correction beyond the ADR-0004 additive perfect-mirror baseline (\(1/R + 1/R'\)), so the solver assembles Z_b = jω·(L_Neumann+mirror) + ΔZ_Sommerfeld. The mirror deduplication \(-1/R'\) is carried in closed form; only the smooth reflected/transmitted remainder is integrated numerically.

Limit checks (built into the test suite)

• \(\sigma_e\to 0\): every kernel collapses to free space — the correction tends to \(-1/R'\) and exactly cancels the ADR-0004 mirror. (At DC the soil is magnetically transparent: no image.)
• \(\sigma_e\to\infty\), overhead pair: \(-\Gamma \to -1\), the net image becomes the anti-parallel PEC image \(-1/R'\) (Carson's baseline).
• \(\sigma_e\to\infty\), buried pair: \(e^{-\gamma R} \to 0\) — the conductor is screened by the surrounding medium and the total external coupling tends to zero.
• Long parallel buried wires: the per-metre limit reproduces Pollaczek's mutual $\frac{j\omega\mu_0}{2\pi}\bigl[K_0(\gamma d)
• K_0(\gamma D') + J_P\bigr]$, whose low-frequency real part is Carson's \(\omega\mu_0/8\) (with the positive sign — the historic kernel produced \(-\omega\mu_0/8\)).

References

• Stratton, J. A. (1941). Electromagnetic Theory, McGraw-Hill, §9-10 — derivation of the half-space vector potential.
• Sommerfeld, A. (1909). Über die Ausbreitung der Wellen in der drahtlosen Telegraphie. Ann. Phys. 28(4), 665–736.
• Wait, J. R. (1972). Electromagnetic Waves in Stratified Media, Pergamon. Ch. 3.
• Tleis, N. D. (2008). Power Systems Modelling and Fault Analysis, Newnes. Ch. 3.

LayeredEarth dataclass

LayeredEarth(
    rhos: tuple[float, ...], thicknesses: tuple[float, ...]
)

Frozen layered-earth configuration for the Sommerfeld kernel.

Attributes:

Name	Type	Description
rhos	tuple[float, ...]	Resistivities \(\rho_1, \dots, \rho_n\) of the layers in \(\Omega\,\mathrm{m}\). The last entry is the semi-infinite bottom layer.
thicknesses	tuple[float, ...]	Thicknesses \(h_1, \dots, h_{n-1}\) in metres. The bottom layer has no thickness (semi-infinite). For len(thicknesses) == len(rhos) - 1.

Notes

Numerical precision contract. The reflection-coefficient evaluators (:func:reflection_coefficient_homogeneous, :func:reflection_coefficient_layered) and every consumer of this dataclass operate in IEEE-754 double precision (FP64). Future hardware-accelerated backends (e.g. an MLX path on Apple silicon) must honour the same precision or the cross-backend cross-check in tests/test_layered_green.py::test_cross_backend_precision will fail. np.complex128 is the default at every entry point and is preserved through the Sommerfeld quadrature; do not silently down-cast to FP32 in a derived backend.

build_sommerfeld_correction_matrix

build_sommerfeld_correction_matrix(
    seg_endpoints: np.ndarray,
    wire_radii: np.ndarray,
    *,
    omega: float,
    earth: LayeredEarth
) -> np.ndarray

Assemble the dense Sommerfeld earth-return correction matrix.

The output is the σ-dependent addition to the perfect-mirror Neumann inductance matrix from :func:groundfield.coupling.inductance.build_inductance_matrix. The two should be added (after multiplying \(L_\text{Neumann}\) by \(j\omega\)):

.. code-block:: python

Z_b = jω · L_Neumann + dZ_Sommerfeld

For earth.n_layers == 1 this is the homogeneous-earth Sommerfeld kernel; for n >= 2 the layered Pollaczek-Wait kernel is used.

Parameters:

Name	Type	Description	Default
seg_endpoints	ndarray	Array of shape (M, 2, 3) with the start- and end-points of every distributed-conductor longitudinal-branch segment.	required
wire_radii	ndarray	Per-branch wire radii. Currently unused (the radius is already in the perfect-mirror diagonal handled by build_inductance_matrix); reserved for future use when the diagonal needs a wire-radius regularisation in the Sommerfeld self-pair integral.	required
omega	float	Angular frequency in rad/s.	required
earth	LayeredEarth	Layered-earth configuration.	required

Returns:

Name	Type	Description
dZ	np.ndarray, shape (M, M), dtype complex	Symmetric Sommerfeld correction matrix in \(\Omega\).

Source code in src/groundfield/coupling/sommerfeld_inductance.py

def build_sommerfeld_correction_matrix(
    seg_endpoints: np.ndarray,        # shape (M, 2, 3)
    wire_radii: np.ndarray,           # shape (M,) — currently unused but kept for API parity
    *,
    omega: float,
    earth: LayeredEarth,
) -> np.ndarray:
    """Assemble the dense Sommerfeld earth-return correction matrix.

    The output is the σ-dependent addition to the perfect-mirror
    Neumann inductance matrix from
    :func:`groundfield.coupling.inductance.build_inductance_matrix`.
    The two should be added (after multiplying $L_\\text{Neumann}$
    by $j\\omega$):

    .. code-block:: python

        Z_b = jω · L_Neumann + dZ_Sommerfeld

    For ``earth.n_layers == 1`` this is the homogeneous-earth
    Sommerfeld kernel; for ``n >= 2`` the layered Pollaczek-Wait
    kernel is used.

    Parameters
    ----------
    seg_endpoints
        Array of shape ``(M, 2, 3)`` with the start- and end-points
        of every distributed-conductor longitudinal-branch segment.
    wire_radii
        Per-branch wire radii. Currently unused (the radius is
        already in the perfect-mirror diagonal handled by
        ``build_inductance_matrix``); reserved for future use when
        the diagonal needs a wire-radius regularisation in the
        Sommerfeld self-pair integral.
    omega
        Angular frequency in rad/s.
    earth
        Layered-earth configuration.

    Returns
    -------
    dZ : np.ndarray, shape (M, M), dtype complex
        Symmetric Sommerfeld correction matrix in $\\Omega$.
    """
    if omega <= 0.0:
        return np.zeros(
            (seg_endpoints.shape[0], seg_endpoints.shape[0]),
            dtype=complex,
        )
    M = seg_endpoints.shape[0]
    dZ = np.zeros((M, M), dtype=complex)
    for i in range(M):
        p1_i = seg_endpoints[i, 0]
        p2_i = seg_endpoints[i, 1]
        for j in range(i, M):
            p1_j = seg_endpoints[j, 0]
            p2_j = seg_endpoints[j, 1]
            val = sommerfeld_pair_integral_layered(
                p1_i, p2_i, p1_j, p2_j,
                omega=omega, earth=earth,
            )
            dZ[i, j] = dZ[j, i] = val
    return dZ

earth_return_correction_homogeneous

earth_return_correction_homogeneous(
    *,
    rho: float,
    z_i: float,
    z_j: float,
    omega: float,
    sigma_earth: float
) -> complex

σ-dependent earth-return correction beyond the additive mirror.

Point-kernel form of the buried–buried Pollaczek correction (audit 2026-07-08, WP-C1; see the module docstring):

\[ \Delta G \;=\; \frac{e^{-\gamma R} - 1}{R} \;-\; \frac{1}{R'} \;+\; \int_0^\infty \frac{\lambda}{u_1}\, \Gamma(\lambda)\, e^{-u_1 (z + z')}\, J_0(\lambda\rho)\,d\lambda, \]

with \(R = \sqrt{\rho^2 + (z - z')^2}\), \(R' = \sqrt{\rho^2 + (z + z')^2}\). This is the σ-dependent piece to be added to the ADR-0004 result (which carries \(1/R + 1/R'\)): the \(-1/R'\) removes the additive mirror, the attenuated direct term carries the buried-path earth-return resistance, and the \(\lambda/u_1\)-weighted reflection is the genuine soil-side boundary response.

Limit checks:

• \(\sigma_e \to 0\): \(\gamma \to 0\), \(\Gamma \to 0\) — correction \(\to -1/R'\) → cancels the ADR-0004 image, total → free space \(1/R\) (at DC the soil is magnetically transparent: no image). ✓
• \(\sigma_e \to \infty\): \(e^{-\gamma R} \to 0\), \(\lambda/u_1 \to 0\) — correction \(\to -1/R - 1/R'\), total → 0 (buried conductor screened by the surrounding medium). ✓ (The historic kernel instead returned the additive mirror in this limit, which is the electrostatic — not the magnetic — image; see the audit report.)

Parameters:

Name	Type	Description	Default
rho	float	Horizontal distance between source and field point in m. Must be ≥ 0.	required
z_i	float	Depths (positive into soil) in m. Must be > 0 for the integral to converge; for wires at the surface (Sunde-equivalent depth = 0) use a small regularisation \(z = \max(|\text{depth}|, r)\) where \(r\) is the wire radius.	required
z_j	float	Depths (positive into soil) in m. Must be > 0 for the integral to converge; for wires at the surface (Sunde-equivalent depth = 0) use a small regularisation \(z = \max(|\text{depth}|, r)\) where \(r\) is the wire radius.	required
omega	float	Angular frequency in rad/s.	required
sigma_earth	float	Earth conductivity in S/m.	required

Returns:

Name	Type	Description
correction	complex	The kernel above, dimensionless (the calling :func:build_sommerfeld_correction_matrix multiplies by \(\mu_0 / (4\pi)\) × ... × line integrations).

Source code in src/groundfield/coupling/sommerfeld_inductance.py

def earth_return_correction_homogeneous(
    *,
    rho: float, z_i: float, z_j: float,
    omega: float, sigma_earth: float,
) -> complex:
    """σ-dependent earth-return correction beyond the additive mirror.

    Point-kernel form of the **buried–buried** Pollaczek correction
    (audit 2026-07-08, WP-C1; see the module docstring):

    $$
    \\Delta G \\;=\\; \\frac{e^{-\\gamma R} - 1}{R}
    \\;-\\; \\frac{1}{R'}
    \\;+\\; \\int_0^\\infty \\frac{\\lambda}{u_1}\\,
    \\Gamma(\\lambda)\\, e^{-u_1 (z + z')}\\,
    J_0(\\lambda\\rho)\\,d\\lambda,
    $$

    with $R = \\sqrt{\\rho^2 + (z - z')^2}$,
    $R' = \\sqrt{\\rho^2 + (z + z')^2}$. This is the σ-dependent
    piece *to be added to the ADR-0004 result* (which carries
    $1/R + 1/R'$): the $-1/R'$ removes the additive mirror, the
    attenuated direct term carries the buried-path earth-return
    resistance, and the $\\lambda/u_1$-weighted reflection is the
    genuine soil-side boundary response.

    Limit checks:

    - $\\sigma_e \\to 0$: $\\gamma \\to 0$, $\\Gamma \\to 0$ —
      correction $\\to -1/R'$ → cancels the ADR-0004 image, total
      → free space $1/R$ (at DC the soil is magnetically
      transparent: no image). ✓
    - $\\sigma_e \\to \\infty$: $e^{-\\gamma R} \\to 0$,
      $\\lambda/u_1 \\to 0$ — correction $\\to -1/R - 1/R'$, total
      → 0 (buried conductor screened by the surrounding medium). ✓
      *(The historic kernel instead returned the additive mirror in
      this limit, which is the electrostatic — not the magnetic —
      image; see the audit report.)*

    Parameters
    ----------
    rho
        Horizontal distance between source and field point in m.
        Must be ≥ 0.
    z_i, z_j
        Depths (positive into soil) in m. Must be > 0 for the
        integral to converge; for wires *at* the surface
        (Sunde-equivalent depth = 0) use a small regularisation
        $z = \\max(|\\text{depth}|, r)$ where $r$ is the wire
        radius.
    omega
        Angular frequency in rad/s.
    sigma_earth
        Earth conductivity in S/m.

    Returns
    -------
    correction : complex
        The kernel above, dimensionless (the calling
        :func:`build_sommerfeld_correction_matrix` multiplies by
        $\\mu_0 / (4\\pi)$ × ... × line integrations).
    """
    if omega <= 0.0 or sigma_earth <= 0.0:
        return 0.0 + 0.0j
    z_sum = z_i + z_j
    if z_sum <= 0.0:
        z_sum = 1e-3
    gamma2 = 1j * omega * MU_0 * sigma_earth
    gamma = complex(np.sqrt(gamma2))
    R_dir = math.sqrt(rho * rho + (z_i - z_j) ** 2)
    R_mir = math.sqrt(rho * rho + z_sum * z_sum)
    if R_dir < 1e-9:
        direct = -gamma
    else:
        direct = (np.exp(-gamma * R_dir) - 1.0) / R_dir
    lambdas, weights = _build_lambda_grid(
        z_sum=z_sum, rho_max=max(rho, 1e-6),
        omega=omega, sigma_top=sigma_earth,
    )
    Gamma = reflection_coefficient_homogeneous(
        lambdas, omega=omega, sigma_earth=sigma_earth,
    )
    u1 = np.sqrt(lambdas * lambdas + gamma2)
    decay = np.exp(-u1 * z_sum)
    bessel = j0(lambdas * rho)
    reflected = complex(
        np.sum(weights * (lambdas / u1) * Gamma * decay * bessel)
    )
    return complex(direct) - 1.0 / max(R_mir, 1e-9) + reflected

earth_return_correction_layered

earth_return_correction_layered(
    *,
    rho: float,
    z_i: float,
    z_j: float,
    omega: float,
    earth: LayeredEarth
) -> complex

Layered-earth analogue of :func:earth_return_correction_homogeneous.

Uses :func:reflection_coefficient_layered for the reflected term; the \(\lambda/u_1\) weight, the exponents and the direct attenuation use the top-layer parameters (conductors live in layer 1). For earth.n_layers == 1 it short-circuits to the single-layer formula.

Source code in src/groundfield/coupling/sommerfeld_inductance.py

def earth_return_correction_layered(
    *,
    rho: float, z_i: float, z_j: float,
    omega: float, earth: LayeredEarth,
) -> complex:
    """Layered-earth analogue of :func:`earth_return_correction_homogeneous`.

    Uses :func:`reflection_coefficient_layered` for the reflected
    term; the $\\lambda/u_1$ weight, the exponents and the direct
    attenuation use the top-layer parameters (conductors live in
    layer 1). For ``earth.n_layers == 1`` it short-circuits to the
    single-layer formula.
    """
    if omega <= 0.0:
        return 0.0 + 0.0j
    if earth.n_layers == 1:
        return earth_return_correction_homogeneous(
            rho=rho, z_i=z_i, z_j=z_j,
            omega=omega, sigma_earth=1.0 / earth.rhos[0],
        )
    z_sum = z_i + z_j
    if z_sum <= 0.0:
        z_sum = 1e-3
    sigma_top = 1.0 / earth.rhos[0]
    gamma2 = 1j * omega * MU_0 * sigma_top
    gamma = complex(np.sqrt(gamma2))
    R_dir = math.sqrt(rho * rho + (z_i - z_j) ** 2)
    R_mir = math.sqrt(rho * rho + z_sum * z_sum)
    if R_dir < 1e-9:
        direct = -gamma
    else:
        direct = (np.exp(-gamma * R_dir) - 1.0) / R_dir
    lambdas, weights = _build_lambda_grid(
        z_sum=z_sum, rho_max=max(rho, 1e-6),
        omega=omega, sigma_top=sigma_top,
    )
    Gamma = reflection_coefficient_layered(lambdas, omega=omega, earth=earth)
    u1 = np.sqrt(lambdas * lambdas + gamma2)
    decay = np.exp(-u1 * z_sum)
    bessel = j0(lambdas * rho)
    reflected = complex(
        np.sum(weights * (lambdas / u1) * Gamma * decay * bessel)
    )
    return complex(direct) - 1.0 / max(R_mir, 1e-9) + reflected

reflection_coefficient_homogeneous

reflection_coefficient_homogeneous(
    lambdas: np.ndarray, *, omega: float, sigma_earth: float
) -> np.ndarray

Magnetic reflection coefficient for a homogeneous half-space.

.. math::

\Gamma_\text{mag}^{(1)}(\lambda) \;=\;
\frac{u_e - \lambda}{u_e + \lambda},
\qquad u_e \;=\; \sqrt{\lambda^2 + j\omega\mu_0\sigma_e}.

Parameters:

Name	Type	Description	Default
lambdas	ndarray	Spectral variable, shape (N_lam,). Must be non-negative.	required
omega	float	Angular frequency in rad/s.	required
sigma_earth	float	Earth conductivity in S/m.	required

Returns:

Name	Type	Description
Gamma	(ndarray, complex, shape(N_lam))	Reflection coefficient at each \(\lambda\).

Source code in src/groundfield/coupling/sommerfeld_inductance.py

def reflection_coefficient_homogeneous(
    lambdas: np.ndarray, *, omega: float, sigma_earth: float,
) -> np.ndarray:
    """Magnetic reflection coefficient for a homogeneous half-space.

    .. math::

        \\Gamma_\\text{mag}^{(1)}(\\lambda) \\;=\\;
        \\frac{u_e - \\lambda}{u_e + \\lambda},
        \\qquad u_e \\;=\\; \\sqrt{\\lambda^2 + j\\omega\\mu_0\\sigma_e}.

    Parameters
    ----------
    lambdas : np.ndarray
        Spectral variable, shape ``(N_lam,)``. Must be non-negative.
    omega
        Angular frequency in rad/s.
    sigma_earth
        Earth conductivity in S/m.

    Returns
    -------
    Gamma : np.ndarray, complex, shape (N_lam,)
        Reflection coefficient at each $\\lambda$.
    """
    if omega <= 0.0 or sigma_earth <= 0.0:
        return np.zeros_like(lambdas, dtype=complex)
    u_e = np.sqrt(lambdas * lambdas + 1j * omega * MU_0 * sigma_earth)
    return (u_e - lambdas) / (u_e + lambdas)

reflection_coefficient_layered

reflection_coefficient_layered(
    lambdas: np.ndarray,
    *,
    omega: float,
    earth: LayeredEarth
) -> np.ndarray

Magnetic reflection coefficient for an \(n\)-layer earth.

Implements the recursive Tagg/Sunde-Wait formula

.. math::

\Gamma_k(\lambda) \;=\;
\frac{u_k - u_{k+1} - (u_k + u_{k+1})\,\Gamma_{k+1}\,e^{-2 u_k h_k}}
     {u_k + u_{k+1} + (u_k - u_{k+1})\,\Gamma_{k+1}\,e^{-2 u_k h_k}},

starting from \(\Gamma_n = 0\) (semi-infinite bottom layer) and walking up to layer 1. The top-layer reflection is \(\Gamma_\text{mag}^{(n)}(\lambda) = (u_e - \lambda)/(u_e+\lambda)\) with \(u_e \to u_1\) in the homogeneous limit, but for \(n>1\) the recursion modifies \(u_1\) effectively. The formulation below is the standard one in Tleis 2008 §3.5.

For \(n=1\) this collapses to :func:reflection_coefficient_homogeneous.

Parameters:

Name	Type	Description	Default
lambdas	(ndarray, shape(N_lam))		required
omega	float	Angular frequency in rad/s.	required
earth	LayeredEarth	Layered-earth configuration.	required

Returns:

Name	Type	Description
Gamma	(ndarray, complex)

Source code in src/groundfield/coupling/sommerfeld_inductance.py

def reflection_coefficient_layered(
    lambdas: np.ndarray, *, omega: float, earth: LayeredEarth,
) -> np.ndarray:
    """Magnetic reflection coefficient for an $n$-layer earth.

    Implements the recursive Tagg/Sunde-Wait formula

    .. math::

        \\Gamma_k(\\lambda) \\;=\\;
        \\frac{u_k - u_{k+1} - (u_k + u_{k+1})\\,\\Gamma_{k+1}\\,e^{-2 u_k h_k}}
             {u_k + u_{k+1} + (u_k - u_{k+1})\\,\\Gamma_{k+1}\\,e^{-2 u_k h_k}},

    starting from $\\Gamma_n = 0$ (semi-infinite bottom layer)
    and walking up to layer 1. The top-layer reflection is
    $\\Gamma_\\text{mag}^{(n)}(\\lambda) = (u_e - \\lambda)/(u_e+\\lambda)$
    with $u_e \\to u_1$ in the homogeneous limit, but for $n>1$
    the recursion modifies $u_1$ effectively. The formulation
    below is the standard one in Tleis 2008 §3.5.

    For $n=1$ this collapses to
    :func:`reflection_coefficient_homogeneous`.

    Parameters
    ----------
    lambdas : np.ndarray, shape (N_lam,)
    omega
        Angular frequency in rad/s.
    earth : LayeredEarth
        Layered-earth configuration.

    Returns
    -------
    Gamma : np.ndarray, complex
    """
    if omega <= 0.0:
        return np.zeros_like(lambdas, dtype=complex)
    n = earth.n_layers
    if n == 1:
        return reflection_coefficient_homogeneous(
            lambdas, omega=omega, sigma_earth=1.0 / earth.rhos[0],
        )
    # Standard recursive multilayer surface reflection (WP-D2 fix,
    # 2026-07-20). Vertical wavenumbers with ``u[0]`` the air (quasi-
    # static, u = lambda) and ``u[k]`` layer k (k = 1..n):
    u = [lambdas.astype(complex)]
    for rho_k in earth.rhos:
        u.append(np.sqrt(lambdas * lambdas + 1j * omega * MU_0 / rho_k))
    # Interface reflection at boundary k (between medium k and k+1), in the
    # magnetic convention ``Gamma_hom = (u_e - lambda)/(u_e + lambda)``:
    #   r_k = (u_{k+1} - u_k) / (u_{k+1} + u_k),  k = 0..n-1.
    # Cascade from the deepest interface (top of the semi-infinite bottom
    # layer) up to the air surface, adding the round-trip phase
    # ``exp(-2 u_{k+1} h_{k+1})`` across each intermediate layer:
    #   Gamma <- (r_k + Gamma·phase) / (1 + r_k·Gamma·phase).
    # This correctly reduces to the homogeneous rho_1 half-space as
    # h_1 -> inf and to the rho_2 half-space as h_1 -> 0. (The buried-in-
    # layer-1 *correction* no longer uses this Gamma for the two-layer
    # buried case — see ``_two_layer_reflected_weight`` — but it is the
    # correct surface reflection for overhead conductors over layered
    # earth and for the n>=3 approximation.)
    gamma = (u[n] - u[n - 1]) / (u[n] + u[n - 1])  # r_{n-1}, deepest interface
    for k in range(n - 2, -1, -1):
        r_k = (u[k + 1] - u[k]) / (u[k + 1] + u[k])
        phase = np.exp(-2.0 * u[k + 1] * earth.thicknesses[k])
        gamma = (r_k + gamma * phase) / (1.0 + r_k * gamma * phase)
    return gamma

sommerfeld_pair_integral_homogeneous

sommerfeld_pair_integral_homogeneous(
    p1_a: np.ndarray,
    p2_a: np.ndarray,
    p1_b: np.ndarray,
    p2_b: np.ndarray,
    *,
    omega: float,
    sigma_earth: float
) -> complex

Integrate the σ-dependent magnetic Green function over a segment pair.

Computes

\[ \Delta Z^{(i,j)}_\text{Sommerfeld} \;=\; \frac{j\omega\mu_0}{4\pi} \int_{C_i}\!\!\int_{C_j} (\hat{l}_i\cdot\hat{l}_j)\, \Delta G_\text{mag}(\vec{r}_i, \vec{r}_j;\,\omega,\sigma_e)\, dl_i\,dl_j, \]

where \(\Delta G_\text{mag}\) is the σ-dependent earth-return correction beyond the ADR-0004 additive-mirror baseline (vanishes for σ → 0 together with the mirror; see the module docstring for the case-dispatched Pollaczek kernels). Used by :func:build_sommerfeld_correction_matrix.

The integration is 16×16 Gauss–Legendre over the two segment parameterisations. The Sommerfeld inner integral is not re-evaluated at every outer node — the spectral kernel is computed once per segment pair on a shared λ-grid and reused (valid because the segments are assumed approximately horizontal so the depth sum varies negligibly along them; see ADR-0006 numerical notes). The closed-form pieces (mirror deduplication \(-1/R'\); buried–buried direct attenuation \((e^{-\gamma R} - 1)/R\)) are evaluated exactly per outer node.

Returns:

Name	Type	Description
Z	complex	Per-pair earth-return correction in \(\Omega\) (already includes the \(j\omega\mu_0/(4\pi)\) pre-factor).

Source code in src/groundfield/coupling/sommerfeld_inductance.py

def sommerfeld_pair_integral_homogeneous(
    p1_a: np.ndarray, p2_a: np.ndarray,
    p1_b: np.ndarray, p2_b: np.ndarray,
    *,
    omega: float, sigma_earth: float,
) -> complex:
    """Integrate the σ-dependent magnetic Green function over a segment pair.

    Computes

    $$
    \\Delta Z^{(i,j)}_\\text{Sommerfeld} \\;=\\; \\frac{j\\omega\\mu_0}{4\\pi}
    \\int_{C_i}\\!\\!\\int_{C_j}
    (\\hat{l}_i\\cdot\\hat{l}_j)\\,
    \\Delta G_\\text{mag}(\\vec{r}_i, \\vec{r}_j;\\,\\omega,\\sigma_e)\\,
    dl_i\\,dl_j,
    $$

    where $\\Delta G_\\text{mag}$ is the σ-dependent earth-return
    correction beyond the ADR-0004 additive-mirror baseline
    (vanishes for σ → 0 together with the mirror; see the module
    docstring for the case-dispatched Pollaczek kernels). Used by
    :func:`build_sommerfeld_correction_matrix`.

    The integration is 16×16 Gauss–Legendre over the two segment
    parameterisations. The Sommerfeld inner integral is **not**
    re-evaluated at every outer node — the spectral kernel is
    computed once per segment pair on a shared λ-grid and reused
    (valid because the segments are assumed approximately horizontal
    so the depth sum varies negligibly along them; see ADR-0006
    numerical notes). The closed-form pieces (mirror deduplication
    $-1/R'$; buried–buried direct attenuation
    $(e^{-\\gamma R} - 1)/R$) are evaluated exactly per outer node.

    Returns
    -------
    Z : complex
        Per-pair earth-return correction in $\\Omega$ (already
        includes the $j\\omega\\mu_0/(4\\pi)$ pre-factor).
    """
    if omega <= 0.0 or sigma_earth <= 0.0:
        return 0.0 + 0.0j

    def _gamma(lambdas: np.ndarray) -> np.ndarray:
        return reflection_coefficient_homogeneous(
            lambdas, omega=omega, sigma_earth=sigma_earth,
        )

    return _pollaczek_pair_integral(
        p1_a, p2_a, p1_b, p2_b,
        omega=omega, sigma_top=sigma_earth, gamma_provider=_gamma,
    )

sommerfeld_pair_integral_layered

sommerfeld_pair_integral_layered(
    p1_a: np.ndarray,
    p2_a: np.ndarray,
    p1_b: np.ndarray,
    p2_b: np.ndarray,
    *,
    omega: float,
    earth: LayeredEarth
) -> complex

Layered-earth analogue of :func:sommerfeld_pair_integral_homogeneous.

The reflected term uses the layered reflection coefficient; the \(\lambda/u_1\) weight, the \(e^{-u_1(\cdot)}\) exponents and the buried direct-term attenuation use the top-layer \(u_1\) / \(\gamma_1\) (the conductors live in layer 1 — deeper layers act through the reflection coefficient only).

Source code in src/groundfield/coupling/sommerfeld_inductance.py

def sommerfeld_pair_integral_layered(
    p1_a: np.ndarray, p2_a: np.ndarray,
    p1_b: np.ndarray, p2_b: np.ndarray,
    *,
    omega: float, earth: LayeredEarth,
) -> complex:
    """Layered-earth analogue of :func:`sommerfeld_pair_integral_homogeneous`.

    The reflected term uses the layered reflection coefficient; the
    $\\lambda/u_1$ weight, the $e^{-u_1(\\cdot)}$ exponents and the
    buried direct-term attenuation use the **top-layer** $u_1$ /
    $\\gamma_1$ (the conductors live in layer 1 — deeper layers act
    through the reflection coefficient only).
    """
    if earth.n_layers == 1:
        return sommerfeld_pair_integral_homogeneous(
            p1_a, p2_a, p1_b, p2_b,
            omega=omega, sigma_earth=1.0 / earth.rhos[0],
        )
    if omega <= 0.0:
        return 0.0 + 0.0j
    sigma_top = 1.0 / earth.rhos[0]

    def _gamma(lambdas: np.ndarray) -> np.ndarray:
        return reflection_coefficient_layered(
            lambdas, omega=omega, earth=earth,
        )

    return _pollaczek_pair_integral(
        p1_a, p2_a, p1_b, p2_b,
        omega=omega, sigma_top=sigma_top, gamma_provider=_gamma,
        earth=earth,
    )

Cross-layer Green's function (ADR-0007 Phase B)

groundfield.coupling.layered_green solves the two-layer matching problem for the electric Green's function — used by image_2layer, mom_sommerfeld, cim, and bem whenever a source / observer pair straddles the upper-layer interface (driven rods, deep meshes, foundation electrodes that cross \(z = h_1\)).

Entry points:

• two_layer_spectral_kernel — the kernel \(\widetilde{G}(\lambda; z, z_s)\) in spectral space,
• two_layer_real_space_kernel — its real-space counterpart \(G(s, z, z_s)\) obtained by Sommerfeld inversion,
• two_layer_layered_correction_real_space / two_layer_layered_correction_group — the layered correction \(\Delta G = G_\text{2-layer} - G_\text{homog}\) for one distance or for a whole reaction block at one \((z, z_s)\) pair,
• two_layer_probe_matrix — the post-processing potential matrix \(\varphi = G\,I\) for arbitrary probe/source layer combinations.

Together with the cross-layer-aware self-action factory (_two_layer_self_kernel_factory) they lift the long-standing z_max < h_1 precondition for n_layers == 2. For n_layers >= 3 the layered backends still emit a documented UserWarning — the n-layer extension is on the roadmap.

Numerical strategy: analytic primary term, quadrature for the rest

The spectral kernel splits into a primary term — the particular solution \(\rho_k/(2\lambda)\,e^{-\lambda|z-z_s|}\) in the layer that contains the source — and a reflected remainder generated by the free surface at \(z = 0\) and the interface at \(z = h_1\). The two behave very differently in \(\lambda\) and are treated differently:

\[ G(s, z, z_s) \;=\; \underbrace{\frac{\rho_k}{2\,\sqrt{s^2 + (z-z_s)^2}} \;+\; \sum_{d_j < d_\text{cut}} \frac{c_j}{2\sqrt{s^2+d_j^2}}} _{\text{closed form}} \;+\; \underbrace{\int_0^{\lambda_\text{max}}\! \Big[\Phi_\text{refl}(\lambda) - \!\!\sum_{d_j < d_\text{cut}}\!\! \tfrac{c_j}{2\lambda}e^{-\lambda d_j}\Big] J_0(\lambda s)\,\lambda\,d\lambda} _{\text{split-grid quadrature}} . \]

The primary term is the only part that does not decay in \(\lambda\) at all, and its Hankel inverse is the exact \(1/R\) field. Since 0.15.0 it is evaluated in closed form.

The reflected remainder decays like \(e^{-\lambda d_\text{img}}\) with \(d_\text{img}\) the smallest image distance — which is not by itself a licence to truncate, because every image distance can degenerate: \(z + z_s \to 0\) for a pair at the air/soil surface, \(|2h_1 - z - z_s| \to 0\) for a pair at the interface, and \(z_l - z_u \to 0\) for a pair straddling it. Up to 0.15.0.dev \(\lambda_\text{max}\) was lambda_max_factor \(/\min(h_1, s+z+z_s)\), which has no relation to \(d_\text{img}\), and those three regimes were wrong by −2.7e-2, +4.1 % and −39.5 % respectively with no warning — the last one growing as the mesh was refined towards \(h_1\).

Since 0.15.0 the truncation is guaranteed instead of assumed. The reflected kernel is written as its exact image superposition \(\sum_j \frac{c_j}{2\lambda}e^{-\lambda d_j}\) (the Sunde/Tagg series, for all four layer combinations); every term with \(\lambda_\text{max} d_j < 30\) is subtracted from the integrand and added back through its exact Hankel inverse \(c_j/(2R_j)\), so the discarded tail is bounded by \(\sum_j |c_j| e^{-\lambda_\text{max} d_j}/(2 d_j) \le e^{-30} \sum_j |c_j| / (2 d_j)\) by construction. Measured against an independently derived closed form, surface pairs, interface pairs, pairs a micrometre below the interface, pairs straddling it and radii from 0 to 20 km all agree to ~3e-10.

Up to 0.14.1 the primary term went through the quadrature as well, which left an \(O((\lambda_\text{max}R)^{-1/2})\) truncation residue oscillating in \(s\) with period \(2\pi/\lambda_\text{max}\): −4.5 % at \(s = 1\) m and ±4 % swings around \(s = 10\) m for an equal-depth pair in \(h_1 = 6.1\) m soil. The homogeneous limit \(\rho_2 = \rho_1 \Rightarrow \rho/2\,(1/r + 1/r')\) is now reproduced to ~1e-10 instead of ~1e-2 (ADR-0007 validation item 1 asks for 1e-9).

At \(s = 0\) with \(z = z_s\) every \(1/R\) pole is clamped at \(R = 1/\lambda_\text{clamp}\), \(\lambda_\text{clamp} =\) lambda_max_factor \(/\min(h_1, s+z+z_s)\) — so the clamp depends only on the knob and the depths, never on the largest radius that happens to share a group. The value stays finite and continuous but is not a physical potential; the caller regularises at the wire radius.

Depths must satisfy \(z, z_s \ge 0\) (positive into the soil). Negative depths raise ValueError; up to 0.15.0.dev they returned nan, which propagated silently into grounding_impedance().

Oscillation resolution and SommerfeldResolutionWarning

The linear region of the \(\lambda\) grid places one 16-node Gauss panel per \(J_0(\lambda s)\) period, so the required node count grows with the oscillation content \(n_\text{osc} = \lambda_\text{max}\,s / 2\pi\). The policy has three tiers:

Regime	Resolution
\(16\,n_\text{osc} \le 2^{16}\) nodes	one panel per period (16 nodes / oscillation)
beyond the soft budget	reduced to the floor of 8 nodes / oscillation
\(8\,n_\text{osc} > 2^{20}\) nodes	hard allocation budget — SommerfeldResolutionWarning

Eight nodes per oscillation is the same adequacy criterion the single-panel fast path uses at small radius. Before 0.15.0 the panel count was capped at 4096 with no adequacy check, so raising lambda_max_factor from 200 to 6000 turned a layered correction of −0.2311 into +0.0349 — a sign flip with no diagnostic.

What lambda_max_factor does since 0.15.0

Because the truncated tail is bounded by the closed-form image split above rather than by making \(\lambda_\text{max}\) large, \(\lambda_\text{max}\) is now chosen by the module and lambda_max_factor is only an upper bound on it plus the \(1/R\) regularisation wavenumber. The choice is:

Step	Rule
requested	\(\lambda_\text{user} =\) lambda_max_factor \(/\min(h_1, s_\text{min}+z+z_s)\)
cost cap	at most \(1024 \cdot\) lambda_max_factor\(/200\) panels, never more than 4096 (the historic allocation)
decay cap	\(120 / d_\text{rest}\) — no point integrating past \(e^{-120}\)
repair	raised back to \(30/d_\text{rest}\) if the image-term budget ran out (then also warns)

Consequences, both intended:

• Raising the knob normally does nothing at all. Measured: the full kernel of a surface pair is identical to ten digits for any factor from 30 to 30 000. That is the resolution of F14 — the knob can no longer make the answer worse, because the answer no longer depends on it.
• Lowering it is cheap and equally exact, because work moves into the closed-form image sum. A convergence check should therefore sweep downwards as well, and see nothing move.
• Because \(\lambda_\text{max}\) no longer grows with the knob or with \(s\), SommerfeldResolutionWarning is no longer reachable from a public entry point; it remains as a guard behind the image-term budget, and is tested through it.

The cost cap is also what fixes the multi-kilometre performance of the remote-injection workload: at 0.15.0.dev a 4 km group cost 42x a 90 m one (the node count grew like \(\lambda_\text{max}\,s\) with nothing bounding it). Measured for a 4000-pair group, \(z = z_s = 0.7\) m, \(h_1 = 5\) m, \(\rho = 1000/30\), default knob:

\(s_\text{max}\)	v0.14.1	0.15.0.dev	0.15.0
90 m	94 ms	131 ms	59 ms
400 m	186 ms	418 ms	49 ms
1 km	190 ms	1 107 ms	66 ms
4 km	199 ms	11 216 ms	71 ms
20 km	196 ms	11 641 ms	75 ms

Note

coupling.sommerfeld_inductance._build_lambda_grid (the magnetic counterpart, ADR-0006) still carries the historic unchecked 4096-panel cap with 8-node panels.

Large-group interpolation

Reaction blocks share one \((z, z_s)\) pair, so the kernel there depends on the horizontal distance only, and what is left for the quadrature after the closed-form split is analytic in \(t = \ln(s + c)\) with \(c = d_\text{rest}\), the decay length of that remainder. Groups above 64 distance pairs are therefore sampled on a log-spaced \(s\)-grid — sized from the group's actual dynamic range at ~12 nodes per e-fold — and interpolated with a cubic spline, turning a block from \(O(n_\text{pairs} \cdot n_\lambda)\) into \(O(n_\text{nodes} \cdot n_\lambda)\). Interpolation is skipped whenever it would need at least as many nodes as the group has pairs. The closed-form primary and image terms are always evaluated per pair, outside the interpolation.

Measured agreement with the exact per-pair contraction is ~1e-6 (worst case 1.1e-5, at a surface pair as \(s \to 0\) where the scalar branch is the less accurate of the two), against ~1e-2 for the fixed 48-node linear interpolation used up to 0.14.1; for two_layer_probe_matrix 8e-7 against a documented 1e-4. Because the node count no longer depends on a fixed constant, the same matrix entry no longer changes with the number of pairs in its block — that discontinuity was ~3 % at the 64-pair threshold and is now below 1.1e-5 for every group size from 60 to 129 pairs.

Do not restate this as ~1e-8

Source comments claiming "~1e-8 relative" and "four orders of margin" stood in layered_green.py at 0.15.0.dev; they were 2–3 orders optimistic and contradicted this page. The numbers above are measured.

layered_green

Cross-layer scalar Green's function for layered soil (ADR-0007).

This module implements the electric (scalar-potential) Green's function for a 2-layer soil where source and observer can be in either layer. It complements the existing image-series solver in :mod:groundfield.solver.image_2layer (which assumes both points in the upper layer) and the Sommerfeld kernel in :mod:groundfield.solver.mom_sommerfeld (same restriction).

Mathematical background

For a 2-layer soil with upper-layer resistivity \(\rho_1\) and thickness \(h_1\), lower-layer resistivity \(\rho_2\) semi-infinite, the scalar Green's function for a unit point source at depth \(z_s\) observed at depth \(z\) takes the Hankel-transform form

\[ \varphi(s, z, z_s) \;=\; \int_0^{\infty}\! \Phi(\lambda, z, z_s)\,J_0(\lambda s)\,\lambda\,d\lambda, \]

where the spectral kernel \(\Phi(\lambda, z, z_s)\) is determined by the boundary-value problem

• \(\partial\Phi/\partial z = 0\) at \(z = 0\) (free surface),
• \(\Phi\) continuous at \(z = h_1\),
• \(\sigma_k\,\partial\Phi/\partial z\) continuous at \(z = h_1\).

In each layer \(\Phi\) is a sum of \(e^{\pm\lambda z}\) plus a particular solution \(e^{-\lambda|z-z_s|}/(2\lambda\sigma_k)\) in the layer that contains the source. The four layer-pair cases — \(\Phi_{uu}, \Phi_{ul}, \Phi_{lu}, \Phi_{ll}\) — correspond to the four (source-layer, observer-layer) combinations and are related by reciprocity (\(\Phi_{ul} = \Phi_{lu}\)).

Implementation

The 4×4 boundary-value system is solved numerically once per \((\lambda, z_s)\) pair (4 unknowns: amplitudes in each layer for each direction). Per call the cost is dominated by \(J_0\) evaluation and the linear solve at each \(\lambda\) node — fully vectorisable with numpy.

The implementation focuses on the 2-layer case (:class:groundfield.soil.models.TwoLayerSoil); the n-layer extension follows the same recursion and is deferred to ADR-0007 Phase C.

Numerical strategy

The Hankel inversion is split into an analytic and a numerical part, because the two halves of \(\Phi\) behave completely differently in \(\lambda\):

1. Primary term — the particular solution \(\rho_k/(2\lambda)\,e^{-\lambda|z-z_s|}\) in the layer that contains the source. It does not decay when \(|z - z_s| \to 0\), and its Hankel inverse is known exactly, \(\rho_k/(2R)\) with \(R = \sqrt{s^2 + (z-z_s)^2}\). It is therefore evaluated in closed form (:func:_primary_term_real_space). Integrating it numerically instead — as versions up to 0.14.1 did — truncates at \(\lambda_\text{max}\) and leaves an \(O((\lambda_\text{max}R)^{-1/2})\) residue that oscillates in \(s\) with period \(2\pi/\lambda_\text{max}\): up to 4.5 % on the total kernel, and non-monotone in lambda_max_factor.
2. Reflected remainder — everything generated by the free surface and the interface. It is a superposition of image terms \(c_j/(2\lambda)\,e^{-\lambda d_j}\) (:func:_image_generation gives the exact coefficients and distances), so it decays like \(e^{-\lambda d_\text{img}}\) with \(d_\text{img} = \min_j d_j\) the smallest image distance — but only in proportion to that distance. Every image distance can degenerate to zero (\(z + z_s \to 0\) for a surface pair, \(|2h_1 - z - z_s| \to 0\) for a pair at the interface, \(z_l - z_u \to 0\) for a cross-layer pair straddling it), and then the remainder does not decay at all and truncating it is not safe. Versions up to 0.15.0.dev sized \(\lambda_\text{max}\) from \(\min(h_1, s + z + z_s)\), which has no relation to \(d_\text{img}\), and were consequently wrong by 0.7 % … 40 % on such pairs with no warning (review pass 9 re-audit). Since 0.15.0 the truncation is guaranteed: every image term with \(\lambda_\text{max} d_j < 30\) is subtracted from the spectral integrand and added back through its exact Hankel inverse \(c_j/(2\sqrt{s^2 + d_j^2})\) (:func:_reflected_image_terms), so the numerically integrated rest provably decays like \(e^{-\lambda d_\text{rest}}\) with \(\lambda_\text{max} d_\text{rest} \ge 30\) and the discarded tail is bounded by \(e^{-30} \approx 10^{-13}\). The rest goes through the split-grid quadrature (:func:_hankel_lambda_grid): logarithmic nodes below the first \(J_0\) oscillation, then one Gauss panel per oscillation period, with a documented resolution floor of eight nodes per period and a :class:SommerfeldResolutionWarning when the allocation budget cannot meet it.

Because the tail bound above is analytic, \(\lambda_\text{max}\) no longer has to be large — it only has to be consistent with the set of subtracted image terms. That is what makes the cost of a remote pair bounded: the node count of the oscillation-resolved grid grows like \(\lambda_\text{max}\,s\), so the module caps \(\lambda_\text{max}\) at :data:_HANKEL_PANEL_TARGET panels (scaled by lambda_max_factor / 200) and subtracts as many image terms as that cap requires (:func:_plan_truncation). A 4 km pair costs the same as a 90 m one, and lambda_max_factor becomes a genuine convergence knob: raising it moves work from the closed form into the quadrature and must not change the answer.

Because the numerically integrated rest is analytic in \(t = \ln(s + c)\) — and, after the image split above, has no scale shorter than \(d_\text{rest}\), which is what \(c\) is set to — reaction blocks with many distance pairs at one \((z, z_s)\) are evaluated on a range-sized log grid and cubic-spline interpolated instead of pair by pair (~1e-6 relative, worst case 3e-5; see :func:two_layer_layered_correction_group).

Limit checks (in the test suite)

• \(\rho_2 = \rho_1\): all four pair kernels reduce to the homogeneous Green's function \(\rho/(4\pi r) + \rho/(4\pi r')\).
• \(\rho_2 \to \infty\) (PEC bottom): for source in upper layer, \(G_{ul} \to 0\) (no penetration into the lower layer); \(G_{uu} \to G_{uu}^{K \to 1}\) — image structure with \(K \to 1\).
• \(\rho_2 \to 0\) (sink): \(G_{uu} \to G_{uu}^{K \to -1}\), spreading resistance drops dramatically.
• Source/observer in upper layer: bit-exact match to the Tagg/Sunde image series in :func:groundfield.solver.image_2layer.solve_image_2layer.

References

• Sunde, E. D. (1968). Earth Conduction Effects in Transmission Systems, Dover, Ch. 3.
• Tagg, G. F. (1964). Earth Resistances, Newnes.
• Wait, J. R. (1972). Electromagnetic Waves in Stratified Media, Pergamon.

SommerfeldResolutionWarning

Bases: UserWarning

The Sommerfeld quadrature could not resolve the \(J_0\) oscillations.

The linear region of the Hankel grid places one Gauss panel per \(J_0(\lambda s)\) period, so the node count grows like \(\lambda_\text{max}\,s\). A hard node budget (:data:_HANKEL_MAX_LINEAR_NODES) bounds the transient allocation; when that budget forces fewer than :data:_HANKEL_MIN_NODES_PER_OSCILLATION nodes per period the quadrature is no longer converged and this warning is raised instead of returning a silently wrong number (a coarse grid on an oscillatory integrand does not merely lose digits — it can flip the sign of the result).

The remedy is always to lower lambda_max_factor or to reduce the horizontal distance range of the group. Lowering the factor does not lose accuracy: the truncated tail is subtracted in closed form (:func:_reflected_image_terms), so a smaller \(\lambda_\text{max}\) only moves work from the quadrature into the exact image sum. Silence with warnings.simplefilter("ignore", SommerfeldResolutionWarning) only if the affected entries are known to be negligible.

SommerfeldTailTruncationWarning

Bases: UserWarning

The \(\lambda\)-truncated tail of the reflected kernel is not bounded.

Sibling of :class:SommerfeldResolutionWarning: that one reports a grid too coarse to follow \(J_0(\lambda s)\), this one a grid too short to have integrated the reflected remainder.

The reflected part of \(\Phi\) is a superposition of image terms \(c_j/(2\lambda)\,e^{-\lambda d_j}\), so truncating the Hankel integral at \(\lambda_\text{max}\) discards \(\sum_j c_j/2 \int_{\lambda_\text{max}}^\infty e^{-\lambda d_j} J_0(\lambda s)\,d\lambda\), bounded by \(\sum_j |c_j| e^{-\lambda_\text{max} d_j}/(2 d_j)\). The bound is only small if \(\lambda_\text{max} d_j\) is large for every retained term, and \(d_j\) has nothing to do with the historic $\lambda_\text{max} = $ lambda_max_factor \(/\min(h_1, s + z + z_s)\): a pair at the air/soil surface (\(z + z_s \to 0\)), at the layer interface (\(|2h_1 - z - z_s| \to 0\)) or straddling it has \(d_j \to 0\) and a tail of order one. Review pass 9 measured −39.5 % on such a pair with no diagnostic; the failure grows as a mesh is refined towards the interface, which is what makes it release-blocking rather than cosmetic.

Every term with $\lambda_\text{max} d_j < $ :data:_HANKEL_TAIL_DECAY_TARGET is therefore taken out of the quadrature and inverted in closed form. This warning fires when that could not be completed — in practice only when the term budget :data:_HANKEL_IMAGE_TERM_BUDGET is exhausted, which needs a very thin top layer, \(|K| \to 1\) and a multi-kilometre span at once. The result is then not bounded by the documented tolerance: lower lambda_max_factor (fewer image terms are needed), reduce the distance range, or raise the budget.

two_layer_layered_correction_group

two_layer_layered_correction_group(
    s_values: np.ndarray,
    z: float,
    z_s: float,
    *,
    rho_1: float,
    rho_2: float,
    h_1: float,
    rho_baseline: float | None = None,
    lambda_max_factor: float = 200.0,
    n_log: int = 32,
    n_lin: int = 96
) -> np.ndarray

Vectorised form of :func:two_layer_layered_correction_real_space.

Evaluates the layered correction for many horizontal distances at one (z, z_s) pair with a single spectral-amplitude solve: the amplitudes depend on (λ, z, z_s) only, so the per-pair cost of the historic scalar loop (one BVP solve per matrix entry — the dominant cost of the ADR-0007 cross-layer reaction-matrix assembly after the WP-B1 oscillation-resolved grids) collapses to one BVP solve per (z, z_s) group plus a J0 contraction.

The shared λ-grid is built for the most demanding radius of the group: the lambda_max_factor bound and the \(1/R\) regularisation from the smallest s (near-field content), oscillation resolution, log-region cap and panel-cost budget from the largest (:func:_plan_truncation).

Accuracy (review pass 9, finding F15)

Groups above :data:_HANKEL_INTERP_THRESHOLD pairs are sampled on a log-spaced s-grid and interpolated. Three changes make that branch agree with the exact per-pair contraction to ~1e-6 relative (worst case measured 1.1e-5, over ten geometries and every group size from 60 to 129 pairs; previously up to 1 % for a 400 m span and 17 % when rho_baseline did not match the source layer):

1. the residual primary term \((\rho_\text{src} - \rho_\text{baseline})/(2\lambda)\, e^{-\lambda|z-z_s|}\) — which is not zero unless the caller matches rho_baseline to the source layer, and which carries a \(1/R\) singularity plus a \(2\pi/\lambda_\text{max}\)-periodic truncation ripple in s — is taken out of the quadrature and added back in closed form, leaving a genuinely analytic integrand;
2. the node count is sized from the actual dynamic range of the group and the interpolant is a cubic spline, so the accuracy no longer depends on the network extent;
3. every reflected image term short enough for the λ-truncation to matter gets the same treatment as the primary term — closed form, per pair, outside the interpolation (:func:_reflected_image_terms). This is what removes the 0.7 % … 40 % errors of surface, interface and interface-straddling pairs, and it also smooths what remains: the interpolated function then has no length scale shorter than d_rest, which is what the log-coordinate offset is set to.

Because interpolation is skipped whenever it would need at least as many nodes as the group has pairs, the exact and the interpolated branch also no longer differ observably at the threshold (the historic ~3 % jump between a 64-pair and a 65-pair block is gone; worst deviation over every size from 60 to 129 pairs is 1.1e-5).

Source code in src/groundfield/coupling/layered_green.py

def two_layer_layered_correction_group(
    s_values: np.ndarray,
    z: float,
    z_s: float,
    *,
    rho_1: float,
    rho_2: float,
    h_1: float,
    rho_baseline: float | None = None,
    lambda_max_factor: float = 200.0,
    n_log: int = 32,
    n_lin: int = 96,
) -> np.ndarray:
    """Vectorised form of :func:`two_layer_layered_correction_real_space`.

    Evaluates the layered correction for **many horizontal distances
    at one (z, z_s) pair** with a single spectral-amplitude solve:
    the amplitudes depend on (λ, z, z_s) only, so the per-pair cost
    of the historic scalar loop (one BVP solve per matrix entry —
    the dominant cost of the ADR-0007 cross-layer reaction-matrix
    assembly after the WP-B1 oscillation-resolved grids) collapses
    to one BVP solve per (z, z_s) group plus a ``J0`` contraction.

    The shared λ-grid is built for the most demanding radius of the
    group: the ``lambda_max_factor`` bound and the $1/R$ regularisation
    from the smallest ``s`` (near-field content), oscillation
    resolution, log-region cap and panel-cost budget from the largest
    (:func:`_plan_truncation`).

    Accuracy (review pass 9, finding F15)
    -------------------------------------
    Groups above :data:`_HANKEL_INTERP_THRESHOLD` pairs are sampled on
    a log-spaced ``s``-grid and interpolated. Three changes make that
    branch agree with the exact per-pair contraction to ~1e-6 relative
    (worst case measured 1.1e-5, over ten geometries and every group
    size from 60 to 129 pairs; previously up to 1 % for a 400 m span
    and 17 % when ``rho_baseline`` did not match the source layer):

    1. the residual primary term
       $(\\rho_\\text{src} - \\rho_\\text{baseline})/(2\\lambda)\\,
       e^{-\\lambda|z-z_s|}$ — which is *not* zero unless the caller
       matches ``rho_baseline`` to the source layer, and which carries a
       $1/R$ singularity plus a $2\\pi/\\lambda_\\text{max}$-periodic
       truncation ripple in ``s`` — is taken out of the quadrature and
       added back in closed form, leaving a genuinely analytic
       integrand;
    2. the node count is sized from the actual dynamic range of the
       group and the interpolant is a cubic spline, so the accuracy no
       longer depends on the network extent;
    3. every reflected image term short enough for the λ-truncation to
       matter gets the same treatment as the primary term — closed form,
       per pair, outside the interpolation
       (:func:`_reflected_image_terms`). This is what removes the
       0.7 % … 40 % errors of surface, interface and interface-straddling
       pairs, and it also *smooths* what remains: the interpolated
       function then has no length scale shorter than ``d_rest``, which
       is what the log-coordinate offset is set to.

    Because interpolation is skipped whenever it would need at least as
    many nodes as the group has pairs, the exact and the interpolated
    branch also no longer differ observably at the threshold (the
    historic ~3 % jump between a 64-pair and a 65-pair block is gone;
    worst deviation over every size from 60 to 129 pairs is 1.1e-5).
    """
    if rho_baseline is None:
        rho_baseline = rho_1
    _validate_depths(z, z_s)
    s_values = np.asarray(s_values, dtype=float)
    if abs(rho_2 - rho_1) < 1e-12 * max(rho_1, 1.0) and \
            abs(rho_baseline - rho_1) < 1e-12 * max(rho_1, 1.0):
        return np.zeros(s_values.shape, dtype=float)
    s_min = float(s_values.min())
    s_max = float(s_values.max())
    plan = _plan_truncation(
        z=z, z_s=z_s, rho_1=rho_1, rho_2=rho_2, h_1=h_1,
        rho_baseline=rho_baseline, s_near=s_min, s_far=s_max,
        lambda_max_factor=lambda_max_factor, warn_stacklevel=3,
    )
    lambdas, weights = _hankel_lambda_grid(
        lambda_max=plan.lambda_max, char_length=plan.char_length, s=s_max,
        n_log=n_log, n_lin=n_lin, warn_stacklevel=3,
    )
    A, B, _, C = _spectral_amplitudes(lambdas, z_s, rho_1, rho_2, h_1)
    A_h, B_h, _, C_h = _spectral_amplitudes(
        lambdas, z_s, rho_baseline, rho_baseline, h_1,
    )
    Phi_lay = _phi_from_stable_amplitudes(
        lambdas, A, B, C, z, z_s, rho_1, rho_2, h_1,
        include_primary=False,
    )
    Phi_hom = _phi_from_stable_amplitudes(
        lambdas, A_h, B_h, C_h, z, z_s, rho_baseline, rho_baseline, h_1,
        include_primary=False,
    )
    w_phi = (
        Phi_lay - Phi_hom - _image_spectral_sum(
            lambdas, plan.coeffs, plan.dists,
        )
    ) * weights * lambdas
    # The residual primary term. Zero — bit-exactly — whenever the
    # caller follows the ADR-0007 recommendation and matches
    # ``rho_baseline`` to the resistivity of the source layer.
    d_rho = _primary_rho(z, z_s, rho_1, rho_2, h_1) - _primary_rho(
        z, z_s, rho_baseline, rho_baseline, h_1,
    )
    c = _log_grid_offset(z, z_s, h_1, d_rest=plan.d_rest)
    out = _hankel_contract_smooth(lambdas, s_values, w_phi, c=c)
    out = out + _image_real_space_sum(
        plan.coeffs, plan.dists, s_values, plan.lambda_clamp,
    )
    if d_rho == 0.0:
        return out
    return out + _primary_term_real_space(
        d_rho, s_values, z - z_s, plan.lambda_clamp,
    )

two_layer_layered_correction_real_space

two_layer_layered_correction_real_space(
    s: float,
    z: float,
    z_s: float,
    *,
    rho_1: float,
    rho_2: float,
    h_1: float,
    rho_baseline: float | None = None,
    lambda_max_factor: float = 200.0,
    n_log: int = 32,
    n_lin: int = 96
) -> float

Real-space layered correction over the homogeneous-rho_1 baseline.

Computes $$ \Delta G(s, z, z_s) \;=\; G_{\text{2-layer}}(s, z, z_s; \rho_1, \rho_2, h_1) \;-\; G_{\text{homog}}(s, z, z_s; \rho_1) $$ by subtracting the spectral kernels before the Hankel integration. The singular source term \(\rho_1/(2\lambda)\,e^{-\lambda|z-z_s|}\) is identical in both spectral kernels and cancels exactly, leaving a smooth, exponentially decaying integrand. Two big advantages:

1. In the homogeneous limit \(\rho_2 = \rho_1\) the spectral difference is identically zero at every \(\lambda\), so the integrator returns 0 to machine precision (no spurious diagonal contribution from finite-Sommerfeld-quadrature error at \(z = z_s\)).
2. The kernel difference $(a_1 - a_1^{\text{hom}})e^{\lambda z}
3. (b_1 - b_1^{\text{hom}})e^{-\lambda z}$ has no \(1/(2\lambda)\) piece, so the small-\(\lambda\) region is finite and easy to integrate.

The caller adds this correction on top of the homogeneous potential computed with the existing :func:groundfield.solver.image._self_corrected_kernel (which correctly handles the line-self diagonal).

The optional rho_baseline parameter selects which homog soil to subtract:

• None (default): subtract homog with \(\rho_1\). Suitable when the calling solver builds \(\phi_\text{hom}\) uniformly with \(\rho_1\).
• explicit float: subtract homog with the given resistivity. Use rho_2 for source segments in the lower layer when the calling solver matches its \(\phi_\text{hom}\) baseline to the source layer (recommended for production-grade solver, ADR-0007 §"Phase A.1").

Returns 0 in the homogeneous limit by construction; otherwise a finite real value with the same units as :func:two_layer_real_space_kernel.

Source code in src/groundfield/coupling/layered_green.py

def two_layer_layered_correction_real_space(
    s: float,
    z: float,
    z_s: float,
    *,
    rho_1: float,
    rho_2: float,
    h_1: float,
    rho_baseline: float | None = None,
    lambda_max_factor: float = 200.0,
    n_log: int = 32,
    n_lin: int = 96,
) -> float:
    """Real-space layered *correction* over the homogeneous-rho_1 baseline.

    Computes
    $$
    \\Delta G(s, z, z_s) \\;=\\;
        G_{\\text{2-layer}}(s, z, z_s; \\rho_1, \\rho_2, h_1)
        \\;-\\; G_{\\text{homog}}(s, z, z_s; \\rho_1)
    $$
    by subtracting the spectral kernels **before** the Hankel
    integration. The singular source term
    $\\rho_1/(2\\lambda)\\,e^{-\\lambda|z-z_s|}$ is identical in
    both spectral kernels and **cancels exactly**, leaving a
    smooth, exponentially decaying integrand. Two big advantages:

    1. In the homogeneous limit $\\rho_2 = \\rho_1$ the spectral
       difference is *identically zero* at every $\\lambda$, so the
       integrator returns 0 to machine precision (no spurious
       diagonal contribution from finite-Sommerfeld-quadrature
       error at $z = z_s$).
    2. The kernel difference $(a_1 - a_1^{\\text{hom}})e^{\\lambda z}
       + (b_1 - b_1^{\\text{hom}})e^{-\\lambda z}$ has no
       $1/(2\\lambda)$ piece, so the small-$\\lambda$ region is
       finite and easy to integrate.

    The caller adds this correction on top of the homogeneous
    potential computed with the existing
    :func:`groundfield.solver.image._self_corrected_kernel` (which
    correctly handles the line-self diagonal).

    The optional ``rho_baseline`` parameter selects which homog
    soil to subtract:

    - ``None`` (default): subtract homog with $\\rho_1$. Suitable
      when the calling solver builds $\\phi_\\text{hom}$ uniformly
      with $\\rho_1$.
    - explicit float: subtract homog with the given resistivity.
      Use ``rho_2`` for source segments in the lower layer when
      the calling solver matches its $\\phi_\\text{hom}$ baseline
      to the source layer (recommended for production-grade solver,
      ADR-0007 §"Phase A.1").

    Returns 0 in the homogeneous limit by construction; otherwise
    a finite real value with the same units as
    :func:`two_layer_real_space_kernel`.
    """
    if rho_baseline is None:
        rho_baseline = rho_1
    if abs(rho_2 - rho_1) < 1e-12 * max(rho_1, 1.0) and \
            abs(rho_baseline - rho_1) < 1e-12 * max(rho_1, 1.0):
        return 0.0

    _validate_depths(z, z_s)
    plan = _plan_truncation(
        z=z, z_s=z_s, rho_1=rho_1, rho_2=rho_2, h_1=h_1,
        rho_baseline=rho_baseline, s_near=s, s_far=s,
        lambda_max_factor=lambda_max_factor, warn_stacklevel=3,
    )
    lambdas, weights = _hankel_lambda_grid(
        lambda_max=plan.lambda_max, char_length=plan.char_length, s=s,
        n_log=n_log, n_lin=n_lin, warn_stacklevel=3,
    )

    # Layered amplitudes
    A, B, _, C = _spectral_amplitudes(
        lambdas, z_s, rho_1, rho_2, h_1,
    )
    # Baseline (homog) amplitudes — both layers at rho_baseline.
    A_h, B_h, _, C_h = _spectral_amplitudes(
        lambdas, z_s, rho_baseline, rho_baseline, h_1,
    )

    # Evaluate the *reflected* Phi for both and subtract. When the
    # source is in the same layer as the baseline-rho region, the
    # singular particular term cancels identically and the residual is
    # bounded. When the baseline differs from rho_at_source (e.g.,
    # rho_baseline=rho_2 with source in upper at z_s<h_1, where the
    # layered Phi has rho_1/(2λ)·exp() but the baseline-homog Phi has
    # rho_2/(2λ)·exp()), the residual carries a
    # (rho_at_source - rho_baseline)/(2R) divergence. That piece is
    # handled in closed form below rather than by quadrature — it is
    # the only non-decaying part of the spectral difference, so leaving
    # it in the integral produced an O((λ_max·R)^{-1/2}) ripple in s.
    # Callers should still select rho_baseline = rho at the source
    # layer (ADR-0007 §"Phase A.1"), which makes the term vanish.
    Phi_lay = _phi_from_stable_amplitudes(
        lambdas, A, B, C, z, z_s, rho_1, rho_2, h_1,
        include_primary=False,
    )
    Phi_hom = _phi_from_stable_amplitudes(
        lambdas, A_h, B_h, C_h, z, z_s, rho_baseline, rho_baseline, h_1,
        include_primary=False,
    )
    delta_Phi = Phi_lay - Phi_hom - _image_spectral_sum(
        lambdas, plan.coeffs, plan.dists,
    )

    integrand = delta_Phi * j0(lambdas * s) * lambdas
    total = float(np.sum(weights * integrand))
    total += float(
        _image_real_space_sum(plan.coeffs, plan.dists, s, plan.lambda_clamp)
    )

    d_rho = _primary_rho(z, z_s, rho_1, rho_2, h_1) - _primary_rho(
        z, z_s, rho_baseline, rho_baseline, h_1,
    )
    if d_rho == 0.0:
        return total
    return total + float(
        _primary_term_real_space(d_rho, s, z - z_s, plan.lambda_clamp)
    )

two_layer_probe_matrix

two_layer_probe_matrix(
    probe_points: np.ndarray,
    source_points: np.ndarray,
    *,
    rho_1: float,
    rho_2: float,
    h_1: float,
    min_distance: float = 0.001,
    lambda_max_factor: float = 200.0
) -> np.ndarray

Potential matrix phi = G @ I for arbitrary layer combinations.

Returns the (M, N) matrix with G[m, n] = potential at probe_points[m] per unit current injected at source_points[n], evaluated with the full two-layer spectral kernel (:func:two_layer_real_space_kernel). Valid for every probe/source layer combination — upper–upper, cross-layer and lower–lower — unlike the Tagg/Sunde uu image series that the fast paths use (audit 2026-07-08, WP-D1: the uu series applied to layer-2 points was measured +7 % at z = 7 m and +25 % at z = 12 m for K = +0.818, h_1 = 5 m).

Cost: one spectral-amplitude solve per (z, z_s) depth group plus a J0 contraction (see "Cost control" below) — reserve for the sub-blocks that the uu series cannot handle.

Notes

The near-singularity is regularised like the image-series paths: the horizontal distance is clamped at min_distance whenever the vertical separation is below min_distance as well.

Cost control: probe/source pairs are grouped by their (z, z_s) depth combination (rounded to 1 µm). Within a group the primary \(1/R\) field is evaluated in closed form at every pair, and the reflected remainder — analytic in \(t = \ln(s + c)\) — is what the interpolation sees. Groups larger than :data:_HANKEL_INTERP_THRESHOLD pairs are sampled on a log-spaced s-grid whose node count is derived from the group's dynamic range (:func:_interp_node_count) and interpolated with a cubic spline, keeping surface-grid evaluations tractable.

Accuracy (review pass 9, finding F15)

The historic fixed 48-node linear interpolation of the total kernel was documented at "relative error ≲ 1e-4" but measured 2.1e-3, because (a) the node count ignored the dynamic range of s and (b) the total kernel contains the \(1/R\) primary term whose λ-truncated quadrature ripples in s with period \(2\pi/\lambda_\text{max}\) — a ripple no coarse interpolant can follow. With the primary term and the degenerate image terms in closed form and range-sized cubic nodes, the interpolated and the exact per-pair branch agree to 8e-7 relative (measured over surface, equal-depth-in-layer-2 and cross-layer probe/source combinations, 400 probes over a 90 m square), i.e. the documented 1e-4 holds with two orders of margin — not the "four orders" and "~1e-8" that stood here at 0.15.0.dev, which were 2–3 orders optimistic and contradicted both the module header and docs/api/coupling.md.

Source code in src/groundfield/coupling/layered_green.py

def two_layer_probe_matrix(
    probe_points: np.ndarray,     # (M, 3)
    source_points: np.ndarray,    # (N, 3)
    *,
    rho_1: float,
    rho_2: float,
    h_1: float,
    min_distance: float = 1e-3,
    lambda_max_factor: float = 200.0,
) -> np.ndarray:
    """Potential matrix ``phi = G @ I`` for arbitrary layer combinations.

    Returns the ``(M, N)`` matrix with ``G[m, n]`` = potential at
    ``probe_points[m]`` per unit current injected at
    ``source_points[n]``, evaluated with the full two-layer spectral
    kernel (:func:`two_layer_real_space_kernel`). Valid for **every**
    probe/source layer combination — upper–upper, cross-layer and
    lower–lower — unlike the Tagg/Sunde uu image series that the
    fast paths use (audit 2026-07-08, WP-D1: the uu series applied to
    layer-2 points was measured +7 % at z = 7 m and +25 % at z = 12 m
    for K = +0.818, h_1 = 5 m).

    Cost: one spectral-amplitude solve per (z, z_s) depth group plus a
    ``J0`` contraction (see "Cost control" below) — reserve for the
    sub-blocks that the uu series cannot handle.

    Notes
    -----
    The near-singularity is regularised like the image-series paths:
    the horizontal distance is clamped at ``min_distance`` whenever
    the vertical separation is below ``min_distance`` as well.

    Cost control: probe/source pairs are grouped by their (z, z_s)
    depth combination (rounded to 1 µm). Within a group the primary
    $1/R$ field is evaluated in closed form at every pair, and the
    *reflected* remainder — analytic in $t = \\ln(s + c)$ — is what the
    interpolation sees. Groups larger than
    :data:`_HANKEL_INTERP_THRESHOLD` pairs are sampled on a log-spaced
    ``s``-grid whose node count is derived from the group's dynamic
    range (:func:`_interp_node_count`) and interpolated with a cubic
    spline, keeping surface-grid evaluations tractable.

    Accuracy (review pass 9, finding F15)
    -------------------------------------
    The historic fixed 48-node linear interpolation of the *total*
    kernel was documented at "relative error ≲ 1e-4" but measured
    2.1e-3, because (a) the node count ignored the dynamic range of
    ``s`` and (b) the total kernel contains the $1/R$ primary term
    whose λ-truncated quadrature ripples in ``s`` with period
    $2\\pi/\\lambda_\\text{max}$ — a ripple no coarse interpolant can
    follow. With the primary term and the degenerate image terms in
    closed form and range-sized cubic nodes, the interpolated and the
    exact per-pair branch agree to **8e-7** relative (measured over
    surface, equal-depth-in-layer-2 and cross-layer probe/source
    combinations, 400 probes over a 90 m square), i.e. the documented
    1e-4 holds with two orders of margin — not the "four orders" and
    "~1e-8" that stood here at 0.15.0.dev, which were 2–3 orders
    optimistic and contradicted both the module header and
    ``docs/api/coupling.md``.
    """
    M = probe_points.shape[0]
    N = source_points.shape[0]
    if M and N:
        _validate_depths(
            float(np.min(probe_points[:, 2])),
            float(np.min(source_points[:, 2])),
            labels=("min(probe_points[:, 2])", "min(source_points[:, 2])"),
        )
    G = np.zeros((M, N), dtype=float)

    def _group_values(
        s_flat: np.ndarray, zm: float, zn: float,
    ) -> np.ndarray:
        """Vectorised kernel over many s at one (z, z_s) pair.

        The spectral amplitudes depend on (λ-grid, z, z_s) only, so
        they are solved **once per group** on a shared λ-grid built
        for the most demanding radius (λ_max from the smallest s,
        oscillation resolution and log-region cap from the largest);
        all s-values then reduce to one ``J0``-matrix contraction of
        the reflected part plus a closed-form primary term.
        """
        plan = _plan_truncation(
            z=zm, z_s=zn, rho_1=rho_1, rho_2=rho_2, h_1=h_1,
            rho_baseline=None, s_near=float(s_flat.min()),
            s_far=float(s_flat.max()),
            lambda_max_factor=lambda_max_factor,
            # 4 frames: _plan_truncation -> _group_values ->
            # two_layer_probe_matrix -> user.
            warn_stacklevel=4,
        )
        lambdas, weights = _hankel_lambda_grid(
            lambda_max=plan.lambda_max, char_length=plan.char_length,
            s=float(s_flat.max()), n_log=32, n_lin=96, warn_stacklevel=4,
        )
        A, B, _, C = _spectral_amplitudes(lambdas, zn, rho_1, rho_2, h_1)
        Phi_refl = _phi_from_stable_amplitudes(
            lambdas, A, B, C, zm, zn, rho_1, rho_2, h_1,
            include_primary=False,
        ) - _image_spectral_sum(lambdas, plan.coeffs, plan.dists)
        w_phi = Phi_refl * weights * lambdas
        c = _log_grid_offset(
            zm, zn, h_1, floor=min_distance, d_rest=plan.d_rest,
        )
        vals = _hankel_contract_smooth(lambdas, s_flat, w_phi, c=c)
        vals = vals + _image_real_space_sum(
            plan.coeffs, plan.dists, s_flat, plan.lambda_clamp,
        )
        rho_primary = _primary_rho(zm, zn, rho_1, rho_2, h_1)
        if rho_primary == 0.0:
            return vals
        return vals + _primary_term_real_space(
            rho_primary, s_flat, zm - zn, plan.lambda_clamp,
        )

    # Group by the (z, z_s) pair — the kernel then depends on s only.
    z_keys = np.round(probe_points[:, 2] / 1e-6).astype(np.int64)
    zs_keys = np.round(source_points[:, 2] / 1e-6).astype(np.int64)
    for zk in np.unique(z_keys):
        rows = np.where(z_keys == zk)[0]
        zm = float(probe_points[rows[0], 2])
        for zsk in np.unique(zs_keys):
            cols = np.where(zs_keys == zsk)[0]
            zn = float(source_points[cols[0], 2])
            dxy = (
                probe_points[rows][:, None, 0:2]
                - source_points[cols][None, :, 0:2]
            )
            s_grid = np.sqrt(np.einsum("mnk,mnk->mn", dxy, dxy))
            if abs(zm - zn) < min_distance:
                np.maximum(s_grid, min_distance, out=s_grid)
            vals = _group_values(s_grid.ravel(), zm, zn)
            G[np.ix_(rows, cols)] = vals.reshape(s_grid.shape)
    # Physical prefactor: phi = G_phys / (2π) per unit current (the
    # spectral kernel carries the rho/2 convention — see
    # solver/mom_sommerfeld.py for the bookkeeping).
    return G / (2.0 * math.pi)

two_layer_real_space_kernel

two_layer_real_space_kernel(
    s: float,
    z: float,
    z_s: float,
    *,
    rho_1: float,
    rho_2: float,
    h_1: float,
    lambda_max_factor: float = 200.0,
    n_log: int = 32,
    n_lin: int = 96
) -> float

Cross-layer real-space scalar Green's function \(\varphi(s, z, z_s)\).

Computes \(\varphi(s, z, z_s) / I = \int_0^\infty \Phi(\lambda, z, z_s)\,J_0(\lambda s)\,\lambda\,d\lambda\) for a unit point source at \((0, 0, z_s)\) observed at \((s, 0, z)\). The result already includes the \(1/(4\pi)\)-style prefactor that is conventional in the existing mom_sommerfeld kernel representation.

Numerical strategy: the primary (source-layer particular) term is evaluated in closed form as \(\rho_k/(2R)\) with \(R = \sqrt{s^2 + (z-z_s)^2}\) (:func:_primary_term_real_space), and so is every reflected image term whose distance is too short for the truncation at \(\lambda_\text{max}\) to be negligible (:func:_reflected_image_terms, :func:_plan_truncation); only the provably fast-decaying rest goes through the split-grid Sommerfeld quadrature (logarithmic for small \(\lambda\), oscillation-resolved for large \(\lambda\)) consistent with :mod:groundfield.coupling.sommerfeld_inductance (ADR-0006).

Measured against an independently derived closed form (image expansion summed to convergence): ~3e-10 relative for surface pairs, interface pairs, pairs a millimetre below the interface, cross-layer pairs straddling it, and radii from 0 to 20 km — and identical to ten digits for lambda_max_factor anywhere in 30 … 30 000. Up to 0.15.0.dev the same geometries were wrong by 0.7 % … 40 % with no warning, and non-monotone in the knob.

Splitting the primary term off matters because it is the only part of \(\Phi\) that does not decay in \(\lambda\): integrating it numerically up to \(\lambda_\text{max}\) left an \(O((\lambda_\text{max}R)^{-1/2})\) oscillatory residue — measured −4.5 % at \(s = 1\) m and ±4 % swings with period \(2\pi/\lambda_\text{max} \approx 0.03\) m for equal-depth pairs, and non-monotone in lambda_max_factor. In the homogeneous limit \(\rho_2 = \rho_1\) the closed form now reproduces \(\rho/2\,(1/r + 1/r')\) to ~1e-13 instead of ~1e-2.

For s = 0 and z = z_s every \(1/R\) pole is clamped at \(R = 1/\lambda_\text{clamp}\) with \(\lambda_\text{clamp} =\) lambda_max_factor \(/\min(h_1, s + z + z_s)\) — so the return value stays finite and continuous and depends only on the knob and the depths, but it is not a physical potential: the caller must regularise at the wire radius.

Parameters:

Name	Type	Description	Default
s	float	Cylindrical radius in m, \(\ge 0\).	required
z	float	Observer / source depths in m, \(\ge 0\) (positive into the soil). Negative depths raise :class:ValueError — they used to return nan.	required
z_s	float	Observer / source depths in m, \(\ge 0\) (positive into the soil). Negative depths raise :class:ValueError — they used to return nan.	required
rho_1	float	Two-layer soil parameters.	required
rho_2	float	Two-layer soil parameters.	required
h_1	float	Two-layer soil parameters.	required
lambda_max_factor	float	Upper bound on the quadrature truncation, as a multiple of \(1/\bar h\) with \(\bar h = \min(h_1, s + z + z_s)\), and the regularisation wavenumber of the \(1/R\) poles. It is not the truncation actually used: since 0.15.0 the truncated tail is bounded by a closed-form image split rather than by making \(\lambda_\text{max}\) large, so \(\lambda_\text{max}\) is lowered to whatever that split allows (:func:_plan_truncation). Consequence: raising this knob is normally a no-op and cannot change the answer, which is the whole point — review pass 9 found the result flipping sign as it was raised. Lowering it below ~30/(smallest image distance) still moves work into the closed form and must likewise leave the answer unchanged; that is the convergence check to run.	200.0
n_log	int	Number of nodes in the logarithmic / linear part of the quadrature grid.	32
n_lin	int	Number of nodes in the logarithmic / linear part of the quadrature grid.	32

Returns:

Name	Type	Description
G	float	Real-space scalar Green's function (dimensionless after dividing by \(\rho\)).

Source code in src/groundfield/coupling/layered_green.py

def two_layer_real_space_kernel(
    s: float,
    z: float,
    z_s: float,
    *,
    rho_1: float,
    rho_2: float,
    h_1: float,
    lambda_max_factor: float = 200.0,
    n_log: int = 32,
    n_lin: int = 96,
) -> float:
    """Cross-layer real-space scalar Green's function $\\varphi(s, z, z_s)$.

    Computes $\\varphi(s, z, z_s) / I = \\int_0^\\infty
    \\Phi(\\lambda, z, z_s)\\,J_0(\\lambda s)\\,\\lambda\\,d\\lambda$
    for a unit point source at $(0, 0, z_s)$ observed at $(s, 0, z)$.
    The result already includes the $1/(4\\pi)$-style prefactor that
    is conventional in the existing `mom_sommerfeld` kernel
    representation.

    Numerical strategy: the primary (source-layer particular) term is
    evaluated in **closed form** as $\\rho_k/(2R)$ with
    $R = \\sqrt{s^2 + (z-z_s)^2}$
    (:func:`_primary_term_real_space`), and so is every reflected image
    term whose distance is too short for the truncation at
    $\\lambda_\\text{max}$ to be negligible
    (:func:`_reflected_image_terms`, :func:`_plan_truncation`); only the
    provably fast-decaying rest goes through the split-grid Sommerfeld
    quadrature (logarithmic for small $\\lambda$, oscillation-resolved
    for large $\\lambda$) consistent with
    :mod:`groundfield.coupling.sommerfeld_inductance` (ADR-0006).

    Measured against an independently derived closed form (image
    expansion summed to convergence): ~3e-10 relative for surface
    pairs, interface pairs, pairs a millimetre below the interface,
    cross-layer pairs straddling it, and radii from 0 to 20 km — and
    *identical to ten digits* for ``lambda_max_factor`` anywhere in
    30 … 30 000. Up to 0.15.0.dev the same geometries were wrong by
    0.7 % … 40 % with no warning, and non-monotone in the knob.

    Splitting the primary term off matters because it is the only part
    of $\\Phi$ that does not decay in $\\lambda$: integrating it
    numerically up to $\\lambda_\\text{max}$ left an
    $O((\\lambda_\\text{max}R)^{-1/2})$ oscillatory residue — measured
    −4.5 % at $s = 1$ m and ±4 % swings with period
    $2\\pi/\\lambda_\\text{max} \\approx 0.03$ m for equal-depth pairs,
    and *non-monotone* in ``lambda_max_factor``. In the homogeneous
    limit $\\rho_2 = \\rho_1$ the closed form now reproduces
    $\\rho/2\\,(1/r + 1/r')$ to ~1e-13 instead of ~1e-2.

    For ``s = 0`` and ``z = z_s`` every $1/R$ pole is clamped at
    $R = 1/\\lambda_\\text{clamp}$ with $\\lambda_\\text{clamp} =$
    ``lambda_max_factor`` $/\\min(h_1, s + z + z_s)$ — so the return
    value stays finite and continuous and depends only on the knob and
    the depths, but it is *not* a physical potential: the caller must
    regularise at the wire radius.

    Parameters
    ----------
    s : float
        Cylindrical radius in m, $\\ge 0$.
    z, z_s : float
        Observer / source depths in m, $\\ge 0$ (positive into the
        soil). Negative depths raise :class:`ValueError` — they used to
        return ``nan``.
    rho_1, rho_2, h_1 : float
        Two-layer soil parameters.
    lambda_max_factor
        *Upper bound* on the quadrature truncation, as a multiple of
        $1/\\bar h$ with $\\bar h = \\min(h_1, s + z + z_s)$, and the
        regularisation wavenumber of the $1/R$ poles. It is **not** the
        truncation actually used: since 0.15.0 the truncated tail is
        bounded by a closed-form image split rather than by making
        $\\lambda_\\text{max}$ large, so $\\lambda_\\text{max}$ is
        lowered to whatever that split allows
        (:func:`_plan_truncation`). Consequence: raising this knob is
        normally a no-op and *cannot* change the answer, which is the
        whole point — review pass 9 found the result flipping sign as it
        was raised. Lowering it below ~30/(smallest image distance)
        still moves work into the closed form and must likewise leave
        the answer unchanged; that is the convergence check to run.
    n_log, n_lin
        Number of nodes in the logarithmic / linear part of the
        quadrature grid.

    Returns
    -------
    G : float
        Real-space scalar Green's function (dimensionless after
        dividing by $\\rho$).
    """
    _validate_depths(z, z_s)
    plan = _plan_truncation(
        z=z, z_s=z_s, rho_1=rho_1, rho_2=rho_2, h_1=h_1,
        rho_baseline=None, s_near=s, s_far=s,
        lambda_max_factor=lambda_max_factor, warn_stacklevel=3,
    )
    lambdas, weights = _hankel_lambda_grid(
        lambda_max=plan.lambda_max, char_length=plan.char_length, s=s,
        n_log=n_log, n_lin=n_lin, warn_stacklevel=3,
    )

    A, B, _, C = _spectral_amplitudes(lambdas, z_s, rho_1, rho_2, h_1)
    Phi_refl = _phi_from_stable_amplitudes(
        lambdas, A, B, C, z, z_s, rho_1, rho_2, h_1,
        include_primary=False,
    ) - _image_spectral_sum(lambdas, plan.coeffs, plan.dists)
    integrand = Phi_refl * j0(lambdas * s) * lambdas
    total = float(np.sum(weights * integrand))
    total += float(
        _image_real_space_sum(plan.coeffs, plan.dists, s, plan.lambda_clamp)
    )

    rho_primary = _primary_rho(z, z_s, rho_1, rho_2, h_1)
    if rho_primary == 0.0:
        return total
    return total + float(
        _primary_term_real_space(rho_primary, s, z - z_s, plan.lambda_clamp)
    )

two_layer_spectral_kernel

two_layer_spectral_kernel(
    lambdas: np.ndarray,
    z: float,
    z_s: float,
    *,
    rho_1: float,
    rho_2: float,
    h_1: float
) -> np.ndarray

Spectral Green's function \(\Phi(\lambda, z, z_s)\).

Returns the value of the spectral kernel at every \(\lambda\) in lambdas for given source depth \(z_s\) and observer depth \(z\). Both depths can be in either layer.

Parameters:

Name	Type	Description	Default
lambdas	ndarray	1-D array of spectral variable values, \(\lambda > 0\) in units of m⁻¹.	required
z	float	Observer / source depth in metres (\(z \ge 0\), positive into the soil).	required
z_s	float	Observer / source depth in metres (\(z \ge 0\), positive into the soil).	required
rho_1	float	Resistivities of upper / lower layer in \(\Omega\,\mathrm{m}\).	required
rho_2	float	Resistivities of upper / lower layer in \(\Omega\,\mathrm{m}\).	required
h_1	float	Upper-layer thickness in m.	required

Returns:

Name	Type	Description
Phi	(ndarray, shape(n_lambda))

Source code in src/groundfield/coupling/layered_green.py

def two_layer_spectral_kernel(
    lambdas: np.ndarray,
    z: float,
    z_s: float,
    *,
    rho_1: float,
    rho_2: float,
    h_1: float,
) -> np.ndarray:
    """Spectral Green's function $\\Phi(\\lambda, z, z_s)$.

    Returns the value of the spectral kernel at every $\\lambda$
    in ``lambdas`` for given source depth $z_s$ and observer depth
    $z$. Both depths can be in either layer.

    Parameters
    ----------
    lambdas
        1-D array of spectral variable values, $\\lambda > 0$ in
        units of m⁻¹.
    z, z_s
        Observer / source depth in metres ($z \\ge 0$, positive into
        the soil).
    rho_1, rho_2
        Resistivities of upper / lower layer in $\\Omega\\,\\mathrm{m}$.
    h_1
        Upper-layer thickness in m.

    Returns
    -------
    Phi : np.ndarray, shape (n_lambda,)
    """
    A, B, _, C = _spectral_amplitudes(
        lambdas, z_s, rho_1, rho_2, h_1,
    )
    return _phi_from_stable_amplitudes(
        lambdas, A, B, C, z, z_s, rho_1, rho_2, h_1,
    )

Earth-conductivity / earth-layer resolvers

These helpers normalise any soil model to a homogeneous \(\sigma_\text{earth}\) (Carson) or a layered structure (Sommerfeld). They are typically invoked by the engine builders, not by user code.

from groundfield.coupling import (
    resolve_earth_conductivity, resolve_earth_layers,
)

sigma = resolve_earth_conductivity(world.soil)        # Carson
layers = resolve_earth_layers(world.soil)             # Sommerfeld

resolve_earth_conductivity

resolve_earth_conductivity(soil) -> float

Resolve earth conductivity \(\sigma_\text{earth}\) from a soil model.

Used by the Carson earth-return correction (ADR-0005) to get a single scalar conductivity from the (potentially layered) soil model. The mapping is:

• :class:~groundfield.soil.models.HomogeneousSoil → \(\sigma = 1/\rho\), exact.
• :class:~groundfield.soil.models.TwoLayerSoil → \(\sigma = 1/\rho_1\) (upper layer), with a runtime warning.
• :class:~groundfield.soil.models.MultiLayerSoil → \(\sigma = 1/\rho_1\) (top layer), with a runtime warning.

The layered cases are an approximation; ADR-0005 references the Pollaczek/Sommerfeld kernel (via mom_sommerfeld) as the rigorous alternative.

Parameters:

Name	Type	Description	Default
soil		Any concrete soil model from :mod:groundfield.soil.models.	required

Returns:

Name	Type	Description
sigma_earth	float	Earth conductivity in S/m.

Raises:

Type	Description
TypeError	If the soil model has no resistivity field that can be interpreted as an upper-layer value.

Source code in src/groundfield/coupling/__init__.py

def resolve_earth_conductivity(soil) -> float:
    """Resolve earth conductivity $\\sigma_\\text{earth}$ from a soil model.

    Used by the Carson earth-return correction (ADR-0005) to get a
    single scalar conductivity from the (potentially layered) soil
    model. The mapping is:

    - :class:`~groundfield.soil.models.HomogeneousSoil` →
      $\\sigma = 1/\\rho$, exact.
    - :class:`~groundfield.soil.models.TwoLayerSoil` →
      $\\sigma = 1/\\rho_1$ (upper layer), with a runtime warning.
    - :class:`~groundfield.soil.models.MultiLayerSoil` →
      $\\sigma = 1/\\rho_1$ (top layer), with a runtime warning.

    The layered cases are an *approximation*; ADR-0005 references the
    Pollaczek/Sommerfeld kernel (via ``mom_sommerfeld``) as the
    rigorous alternative.

    Parameters
    ----------
    soil
        Any concrete soil model from
        :mod:`groundfield.soil.models`.

    Returns
    -------
    sigma_earth : float
        Earth conductivity in S/m.

    Raises
    ------
    TypeError
        If the soil model has no resistivity field that can be
        interpreted as an upper-layer value.
    """
    from groundfield.soil.models import (
        HomogeneousSoil,
        MultiLayerSoil,
        TwoLayerSoil,
    )

    if isinstance(soil, HomogeneousSoil):
        return 1.0 / float(soil.resistivity)
    if isinstance(soil, TwoLayerSoil):
        warnings.warn(
            "Carson series uses upper-layer rho_1 for a TwoLayerSoil — "
            "this is an approximation. For a rigorous result on layered "
            "soils use the Pollaczek/Sommerfeld kernel via "
            "backend='mom_sommerfeld'. See ADR-0005 §'Layered earth'.",
            UserWarning,
            stacklevel=2,
        )
        return 1.0 / float(soil.rho_1)
    if isinstance(soil, MultiLayerSoil):
        warnings.warn(
            "Carson series uses top-layer rho for a MultiLayerSoil — "
            "this is an approximation. For a rigorous result on layered "
            "soils use the Pollaczek/Sommerfeld kernel via "
            "backend='mom_sommerfeld'. See ADR-0005 §'Layered earth'.",
            UserWarning,
            stacklevel=2,
        )
        return 1.0 / float(soil.layers[0].resistivity)
    raise TypeError(
        f"Unsupported soil type for Carson series: {type(soil).__name__}"
    )

resolve_earth_layers

resolve_earth_layers(soil) -> 'LayeredEarth'

Resolve a soil model into a :class:LayeredEarth configuration.

Used by the Sommerfeld earth-return kernel (ADR-0006). Unlike :func:resolve_earth_conductivity, this preserves the full layered structure — no warning is emitted, and layered soils are handled natively by the Pollaczek/Wait kernel.

Parameters:

Name	Type	Description	Default
soil		Any concrete soil model from :mod:groundfield.soil.models.	required

Returns:

Type	Description
LayeredEarth	Frozen layered-earth configuration (resistivities and layer thicknesses) consumable by :func:groundfield.coupling.sommerfeld_inductance.build_sommerfeld_correction_matrix.

Source code in src/groundfield/coupling/__init__.py

def resolve_earth_layers(soil) -> "LayeredEarth":
    """Resolve a soil model into a :class:`LayeredEarth` configuration.

    Used by the Sommerfeld earth-return kernel (ADR-0006). Unlike
    :func:`resolve_earth_conductivity`, this preserves the full
    layered structure — no warning is emitted, and layered soils
    are handled natively by the Pollaczek/Wait kernel.

    Parameters
    ----------
    soil
        Any concrete soil model from
        :mod:`groundfield.soil.models`.

    Returns
    -------
    LayeredEarth
        Frozen layered-earth configuration (resistivities and layer
        thicknesses) consumable by
        :func:`groundfield.coupling.sommerfeld_inductance.build_sommerfeld_correction_matrix`.
    """
    from groundfield.coupling.sommerfeld_inductance import LayeredEarth
    from groundfield.soil.models import (
        HomogeneousSoil,
        MultiLayerSoil,
        TwoLayerSoil,
    )

    if isinstance(soil, HomogeneousSoil):
        return LayeredEarth(
            rhos=(float(soil.resistivity),), thicknesses=(),
        )
    if isinstance(soil, TwoLayerSoil):
        return LayeredEarth(
            rhos=(float(soil.rho_1), float(soil.rho_2)),
            thicknesses=(float(soil.h_1),),
        )
    if isinstance(soil, MultiLayerSoil):
        rhos = tuple(float(layer.resistivity) for layer in soil.layers)
        thicknesses = tuple(
            float(layer.thickness) for layer in soil.layers[:-1]
        )
        return LayeredEarth(rhos=rhos, thicknesses=thicknesses)
    raise TypeError(
        f"Unsupported soil type for Sommerfeld kernel: "
        f"{type(soil).__name__}"
    )