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\)).

Two 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.

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.

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.

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.

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: λ_max from the smallest s (near-field content), oscillation resolution and log-region cap from the largest.

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: ``λ_max`` from the smallest ``s`` (near-field content),
    oscillation resolution and log-region cap from the largest.
    """
    if rho_baseline is None:
        rho_baseline = rho_1
    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())
    char_min = max(min(h_1, s_min + z + z_s + 1e-9), 1e-3)
    lambdas, weights = _hankel_lambda_grid(
        char_length=char_min, s=s_max,
        lambda_max_factor=lambda_max_factor, n_log=n_log, n_lin=n_lin,
    )
    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,
    )
    Phi_hom = _phi_from_stable_amplitudes(
        lambdas, A_h, B_h, C_h, z, z_s, rho_baseline, rho_baseline, h_1,
    )
    w_phi = (Phi_lay - Phi_hom) * weights * lambdas
    if s_values.size <= _HANKEL_INTERP_THRESHOLD:
        return _hankel_contract(lambdas, s_values, w_phi)
    # Large group (a whole-grid reaction block at one depth pair): the
    # correction is smooth and monotone in ``s`` (the singular 1/r part
    # cancels in ``Phi_lay - Phi_hom``), so sample it on a log-spaced
    # ``s``-grid and interpolate — same policy as ``two_layer_probe_matrix``.
    c = max(abs(z - z_s), 1e-3)
    t_lo = math.log(s_min + c)
    t_hi = math.log(max(s_max + c, s_min + c + 1e-9))
    t_nodes = np.linspace(t_lo, t_hi, _HANKEL_INTERP_NODES)
    s_nodes = np.maximum(np.exp(t_nodes) - c, 0.0)
    g_nodes = _hankel_contract(lambdas, s_nodes, w_phi)
    flat = np.interp(np.log(s_values.ravel() + c), t_nodes, g_nodes)
    return flat.reshape(s_values.shape)

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

    char_length = max(min(h_1, s + z + z_s + 1e-9), 1e-3)
    lambdas, weights = _hankel_lambda_grid(
        char_length=char_length, s=s,
        lambda_max_factor=lambda_max_factor, n_log=n_log, n_lin=n_lin,
    )

    # 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 full Phi for both and subtract. When the source
    # is in the same layer as the baseline-rho region, the
    # singular particular term cancels 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)/(4πr) divergence — but the
    # caller is responsible for selecting rho_baseline = rho at
    # source layer to avoid this.
    Phi_lay = _phi_from_stable_amplitudes(
        lambdas, A, B, C, z, z_s, rho_1, rho_2, h_1,
    )
    Phi_hom = _phi_from_stable_amplitudes(
        lambdas, A_h, B_h, C_h, z, z_s, rho_baseline, rho_baseline, h_1,
    )
    delta_Phi = Phi_lay - Phi_hom

    integrand = delta_Phi * j0(lambdas * s) * lambdas
    return float(np.sum(weights * integrand))

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 scalar Sommerfeld quadrature per (probe, source) pair — 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 kernel is a smooth, monotone function of the horizontal distance s alone; groups larger than 24 pairs are evaluated on 24 log-spaced s-nodes and interpolated (relative error ≲ 1e-4 — far below the physical validity of the post-processing evaluation), keeping surface-grid evaluations tractable.

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 scalar Sommerfeld quadrature per (probe, source) pair —
    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 kernel is
    a smooth, monotone function of the horizontal distance ``s``
    alone; groups larger than 24 pairs are evaluated on 24 log-spaced
    ``s``-nodes and interpolated (relative error ≲ 1e-4 — far below
    the physical validity of the post-processing evaluation), keeping
    surface-grid evaluations tractable.
    """
    M = probe_points.shape[0]
    N = source_points.shape[0]
    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.
        """
        s_min = float(s_flat.min())
        s_max = float(s_flat.max())
        char_min = max(min(h_1, s_min + zm + zn + 1e-9), 1e-3)
        lambdas, weights = _hankel_lambda_grid(
            char_length=char_min, s=s_max,
            lambda_max_factor=lambda_max_factor, n_log=32, n_lin=96,
        )
        Phi = two_layer_spectral_kernel(
            lambdas, zm, zn, rho_1=rho_1, rho_2=rho_2, h_1=h_1,
        )
        w_phi = Phi * weights * lambdas
        if s_flat.size <= 64:
            return j0(lambdas[None, :] * s_flat[:, None]) @ w_phi
        # Large groups (surface grids): evaluate on 48 log-spaced
        # s-nodes and interpolate — the Hankel integral is smooth and
        # monotone in s.
        c = max(abs(zm - zn), min_distance)
        t_lo = math.log(s_min + c)
        t_hi = math.log(max(s_max + c, s_min + c + 1e-9))
        t_nodes = np.linspace(t_lo, t_hi, 48)
        s_nodes = np.maximum(np.exp(t_nodes) - c, 0.0)
        g_nodes = j0(lambdas[None, :] * s_nodes[:, None]) @ w_phi
        return np.interp(np.log(s_flat + c), t_nodes, g_nodes)

    # 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: split-grid Sommerfeld quadrature (logarithmic for small \(\lambda\), uniform-with-Bessel-resolution for large \(\lambda\)) consistent with :mod:groundfield.coupling.sommerfeld_inductance (ADR-0006).

For s = 0 (coaxial limit) the kernel diverges logarithmically when z = z_s; 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.	required
z_s	float	Observer / source depths in m.	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 of the quadrature, given as a multiple of \(1 / \bar h\) with \(\bar h = \min(h_1, s + z + z_s)\).	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: split-grid Sommerfeld quadrature
    (logarithmic for small $\\lambda$, uniform-with-Bessel-resolution
    for large $\\lambda$) consistent with
    :mod:`groundfield.coupling.sommerfeld_inductance` (ADR-0006).

    For ``s = 0`` (coaxial limit) the kernel diverges logarithmically
    when ``z = z_s``; 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.
    rho_1, rho_2, h_1 : float
        Two-layer soil parameters.
    lambda_max_factor
        Upper bound of the quadrature, given as a multiple of
        $1 / \\bar h$ with $\\bar h = \\min(h_1, s + z + z_s)$.
    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$).
    """
    char_length = max(min(h_1, s + z + z_s + 1e-9), 1e-3)
    lambdas, weights = _hankel_lambda_grid(
        char_length=char_length, s=s,
        lambda_max_factor=lambda_max_factor, n_log=n_log, n_lin=n_lin,
    )

    Phi = two_layer_spectral_kernel(
        lambdas, z, z_s, rho_1=rho_1, rho_2=rho_2, h_1=h_1,
    )
    integrand = Phi * j0(lambdas * s) * lambdas
    return float(np.sum(weights * integrand))

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__}"
    )