Solver

The solver layer holds the numerical core. Available backends:

Backend	Soil model	Method	Status
image	HomogeneousSoil	image-charge sum (point sources + line self-action)	implemented
image_2layer	TwoLayerSoil	Tagg/Sunde image-charge series	implemented
image_nlayer	HomogeneousSoil, TwoLayerSoil, MultiLayerSoil	image-charge dispatcher (delegates to image for n=1, to image_2layer for n=2; raises for n ≥ 3)	implemented
cim	any layered	Complex Image Method (matrix-pencil fit of \(\Gamma_1(\lambda)\))	implemented
mom	HomogeneousSoil or TwoLayerSoil	Galerkin Method-of-Moments on segment level (independent resolution scheme over the same Green's-function kernels)	implemented
mom_sommerfeld	any layered	Galerkin MoM with direct Sommerfeld quadrature (reference engine)	implemented
bem	any layered	Boundary-element collocation with the CIM kernel	implemented
fem	any layered	Axisymmetric volume PDE with equivalent-hemisphere reduction	implemented

Mathematical / physical model

For a point current source \(I\) at \(z_s\) in the upper layer of a horizontally stratified half-space (insulating soil surface at \(z = 0\)), every backend evaluates the same quasi-static Sommerfeld representation of the potential:

\[ \varphi(s, z) \;=\; \frac{\rho_1\, I}{4\pi} \int_0^{\infty} \bigl[ e^{-\lambda |z - z_s|} + \Gamma_1(\lambda)\, e^{-\lambda (z + z_s)} \bigr]\, J_0(\lambda s)\, d\lambda, \]

with the upward-looking reflection \(\Gamma_1(\lambda)\) built recursively from the bottom up (groundfield.solver._layered.reflection_gamma). The engines differ only in how they evaluate this integral:

• closed-form real images (image, image_2layer, image_nlayer);
• closed-form complex images (cim);
• direct numerical quadrature (mom_sommerfeld);
• volume PDE (fem).

ADR-0002 (docs/adr/0002-engine-family.md) records the selection heuristic between the engines.

Auto-dispatch

Engine.solve automatically forwards backend="image":

• to image_2layer if world.soil is a TwoLayerSoil;
• to image_nlayer if world.soil is a MultiLayerSoil.

Notebooks therefore keep working unchanged when the soil model is swapped.

Example

import groundfield as gf

# Build a small world (single ring electrode in two-layer soil).
soil  = gf.TwoLayerSoil(rho_1=100.0, rho_2=500.0, h_1=2.0)
world = gf.create_world(soil=soil)
gf.create_electrode(
    world, "ring", name="g1",
    center=(0.0, 0.0, 0.8), radius=5.0, wire_radius=0.005,
)
gf.create_source(world, attached_to="g1", magnitude=1.0)

# Create an engine and solve. Auto-dispatch hands `image` over to
# `image_2layer` because the soil is two-layer.
engine = gf.create_engine(
    backend="image",
    frequencies=[50.0, 150.0, 250.0],
    segment_length=0.05,
)
result = world.solve(engine)

# Cluster impedance of the ring electrode.
print(result.cluster_impedance("g1"))

The same World can be solved with any of the eight backends — compare_engines(world, engines={"image": ..., "mom": ...}) reports their cluster-impedance agreement (cross-validation rules below).

Frequency-list order

Engine.frequencies is order-preserving. The list is iterated verbatim and the same order propagates to FieldResult.frequencies and to every per-frequency column in the post-processing helpers (postprocess.sweep, the CSV writers, fit_to_sympy, etc.). A non-monotonic list — e.g. [5000, 50] — is accepted but raises a UserWarning so the convention is visible. Use the explicit opt-in if the order is intentional:

engine = (
    gf.create_engine(backend="image")
    .with_frequencies(5000.0, 50.0, preserve_order=True)
)

Engine.with_frequencies returns a fresh Engine instance; the receiver is not mutated.

Cross-engine validation

groundfield.compare_engines(world, engines={...}) runs the same world through several engines and reports the cluster-impedance agreement. The same rule is enforced by tests/test_cross_engines_extended.py:

• For homogeneous worlds every engine must agree to within 5 % (10 % for fem).
• For 2-layer worlds the closed-form / image / MoM engines must agree with each other to within 5 %.
• A monotone \(\rho_2\) sweep at fixed \(\rho_1\) must produce a monotonically increasing cluster impedance for every engine.

solver

Numerical field solver.

This subpackage forms the computational core of groundfield. The key quantity is the complex potential phi(r, f) in the soil and on the conductor surfaces, evaluated per frequency in the phasor domain. The default solution method for homogeneous soil is the closed-form image-charge sum; for layered soil the Tagg/Sunde image series. A Method-of-Moments backend with the layered Green's function and a finite-element backend (scikit-fem) are reserved.

Contents

Engine Top-level configuration of the numerical kernel: backend choice, frequency list, mesh resolution, tolerances. solve(world) runs the simulation. Backend Literal type listing the available backends ("image", "image_2layer", "mom", "fem"). FieldResult, PointSource Result objects. solve_image, solve_image_2layer Backend entry points (usually called via Engine.solve).

Guiding principle

The PDE / field model is a reference, not the end product. Every solution must expose the quantities required by groundinsight (input impedance, transfer impedances, rho-f curve) for the reduction step.

Engine

Bases: BaseModel

Configuration of the numerical kernel.

Attributes:

Name	Type	Description
backend	Backend	Numerical method. One of image, image_2layer, image_nlayer, cim, mom, mom_sommerfeld, bem, fem. See the module docstring above for the family overview, and ADR-0002 for the selection heuristic.
frequencies	list[float]	Frequency list in Hz. Default [50.0].
segment_length	float	Maximum segment length used to discretise the geometry, in m.
tolerance	float	Relative convergence threshold for iterative solvers.
max_iterations	int	Maximum iterations for iterative solvers.
earth_inductive_model	EarthInductiveModel	How the earth contributes to the inductive coupling between distributed-conductor segments (only effective when at least one conductor sets inductance_model = "neumann"): "perfect_mirror" (ADR-0004, the pre-0.13 default) — the earth is a perfect magnetic mirror. Cheap, frequency-independent inductance assembly, but its additive image is far from the f < 1 kHz physics of buried conductors (audit finding F5); kept as a fast diagnostic baseline. "carson_series" (ADR-0005) — the earth has finite conductivity. Adds Carson 1926's per-meter earth-return correction \(\Delta Z_\text{Carson}(\omega)\) on top of the perfect-mirror Neumann matrix, scaled by the geometric length of each segment-pair. Asymptotically correct for long parallel wires over homogeneous earth; an approximation for short wires or layered soils. "sommerfeld" (ADR-0006, Pollaczek kernel since the 2026-07-08 audit, WP-C1) — geometric integration of the σ-dependent vector-potential Green's function over the actual segment-pair geometry, dispatched by conductor side: buried pairs carry the in-soil attenuation \(e^{-\gamma R}\) (and with it the positive earth-return resistance \(\approx \omega\mu_0/8\) per metre in the long-wire limit), overhead pairs use the air-side reflection (Carson's geometry). Supports layered earth. Limits: \(\sigma \to 0\) → free space (no image); \(\sigma \to \infty\) → screened (buried) / anti-parallel PEC image (overhead). Validated against Pollaczek's analytic \(K_0\) mutual and Carson's line impedance (tests/test_earth_return_benchmark.py). The default since 0.13.0 — the rigorous choice whenever a quantitative Z(f) matters.

solve

solve(world: 'World') -> FieldResult

Run the simulation with the configured backend.

Dispatches the assembled world to the backend selected by :attr:backend. When backend == "image" the backend is auto-forwarded to "image_2layer" for a :class:~groundfield.soil.models.TwoLayerSoil and to "image_nlayer" for a :class:~groundfield.soil.models.MultiLayerSoil, so notebooks written for the homogeneous case keep working when the soil model is replaced.

Parameters:

Name	Type	Description	Default
world	World	The assembled world (soil, electrodes, conductors and at least one source) the backend should solve.	required

Returns:

Type	Description
FieldResult	Container with per-frequency potentials, currents, cluster impedances and post-processing metadata.

Raises:

Type	Description
ValueError	If world.soil is unset, or world contains no electrodes.

Source code in src/groundfield/solver/engine.py

def solve(self, world: "World") -> FieldResult:
    """Run the simulation with the configured backend.

    Dispatches the assembled ``world`` to the backend selected by
    :attr:`backend`. When ``backend == "image"`` the backend is
    auto-forwarded to ``"image_2layer"`` for a
    :class:`~groundfield.soil.models.TwoLayerSoil` and to
    ``"image_nlayer"`` for a
    :class:`~groundfield.soil.models.MultiLayerSoil`, so notebooks
    written for the homogeneous case keep working when the soil
    model is replaced.

    Parameters
    ----------
    world : World
        The assembled world (soil, electrodes, conductors and at
        least one source) the backend should solve.

    Returns
    -------
    FieldResult
        Container with per-frequency potentials, currents, cluster
        impedances and post-processing metadata.

    Raises
    ------
    ValueError
        If ``world.soil`` is unset, or ``world`` contains no
        electrodes.
    """
    if world.soil is None:
        raise ValueError(
            "World has no soil model. Set one before calling solve()."
        )
    if not world.electrodes:
        raise ValueError("World contains no electrodes.")

    _log.info(
        "Engine.solve: backend=%s, n_electrodes=%d, n_freq=%d",
        self.backend,
        len(world.electrodes),
        len(self.frequencies),
    )

    # Pre-flight discretisation check (advisory, never raises).
    # Surfaces thin-wire violations, degenerate electrode sizes
    # and segment-budget overruns as warnings *before* the solve
    # spends minutes on a misconfigured world.
    from groundfield.diagnostics import check_segment_resolution

    for finding in check_segment_resolution(world, self):
        warnings.warn(finding, UserWarning, stacklevel=2)

    # Auto-forwarding: ``backend="image"`` transparently picks
    # ``image_2layer`` for a 2-layer soil and ``image_nlayer`` for
    # a multilayer soil. Notebooks therefore do not need to change
    # the backend string when the soil model changes.
    from groundfield.soil.models import MultiLayerSoil, TwoLayerSoil

    effective_backend = self.backend
    if effective_backend == "image":
        if isinstance(world.soil, TwoLayerSoil):
            effective_backend = "image_2layer"
            _log.info(
                "Engine.solve: TwoLayerSoil detected, "
                "switching automatically to 'image_2layer'."
            )
        elif isinstance(world.soil, MultiLayerSoil):
            effective_backend = "image_nlayer"
            _log.info(
                "Engine.solve: MultiLayerSoil detected, "
                "switching automatically to 'image_nlayer'."
            )

    if effective_backend == "image":
        from groundfield.solver.image import solve_image

        return solve_image(world, self)

    if effective_backend == "image_2layer":
        from groundfield.solver.image_2layer import solve_image_2layer

        return solve_image_2layer(
            world, self,
            max_terms=self.image_max_terms, tol=self.image_series_tol,
        )

    if effective_backend == "image_nlayer":
        from groundfield.solver.image_nlayer import solve_image_nlayer

        return solve_image_nlayer(world, self)

    if effective_backend == "cim":
        from groundfield.solver.cim import solve_cim

        return solve_cim(world, self)

    if effective_backend == "mom":
        from groundfield.solver.mom import solve_mom

        return solve_mom(world, self)

    if effective_backend == "mom_sommerfeld":
        from groundfield.solver.mom_sommerfeld import solve_mom_sommerfeld

        return solve_mom_sommerfeld(world, self)

    if effective_backend == "bem":
        from groundfield.solver.bem import solve_bem

        return solve_bem(world, self)

    if effective_backend == "fem":
        from groundfield.solver.fem import solve_fem

        return solve_fem(world, self)

    raise ValueError(f"Unknown backend '{self.backend}'.")

with_frequencies

with_frequencies(
    *frequencies: float, preserve_order: bool = False
) -> "Engine"

Return a copy of this engine with a new frequencies list.

Parameters:

Name	Type	Description	Default
*frequencies	float	One or more frequencies in Hz. Variadic to make the common case readable: engine.with_frequencies(50, 5000).	()
preserve_order	bool	If True, the frequencies are stored verbatim and no non-monotonic warning is raised even for a decreasing or duplicate-bearing list. This is the explicit opt-in for sweeps that require a specific iteration order. If False (default) the standard :func:_validate_frequencies checks run.	False

Returns:

Type	Description
Engine	A new :class:Engine instance — the receiver is not mutated.

Examples:

>>> eng = Engine(backend="image").with_frequencies(50, 5000,
...                                                preserve_order=True)
>>> eng.frequencies
[50.0, 5000.0]

Source code in src/groundfield/solver/engine.py

def with_frequencies(
    self,
    *frequencies: float,
    preserve_order: bool = False,
) -> "Engine":
    """Return a copy of this engine with a new ``frequencies`` list.

    Parameters
    ----------
    *frequencies
        One or more frequencies in Hz. Variadic to make the
        common case readable: ``engine.with_frequencies(50, 5000)``.
    preserve_order
        If ``True``, the frequencies are stored verbatim and no
        non-monotonic warning is raised even for a decreasing or
        duplicate-bearing list. This is the explicit opt-in for
        sweeps that require a specific iteration order. If ``False``
        (default) the standard :func:`_validate_frequencies` checks
        run.

    Returns
    -------
    Engine
        A new :class:`Engine` instance — the receiver is **not**
        mutated.

    Examples
    --------
    >>> eng = Engine(backend="image").with_frequencies(50, 5000,
    ...                                                preserve_order=True)
    >>> eng.frequencies
    [50.0, 5000.0]
    """
    freqs = [float(f) for f in frequencies]
    if preserve_order:
        # Bypass the validator's monotone-check by using
        # ``model_construct`` for the field. We still want type /
        # positivity checks, so do them explicitly here.
        if not freqs:
            raise ValueError(
                "with_frequencies(): at least one frequency required."
            )
        for f in freqs:
            if not (f >= 0.0):
                raise ValueError(
                    "with_frequencies(): frequencies must be "
                    f"non-negative, got {f}."
                )
        data = self.model_dump()
        data["frequencies"] = freqs
        # Re-validate the rest, but inject the preserved list
        # directly so the field-validator's monotone-warning is
        # silenced.
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", EngineFrequencyOrderWarning)
            return self.__class__.model_validate(data)
    return self.model_copy(update={"frequencies": freqs})

FieldResult

Bases: BaseModel

Result of a field computation.

cluster_impedance

cluster_impedance(electrode_name: str) -> list[complex]

Grounding impedance of the galvanic cluster containing an electrode.

Definition: \(Z_{\text{cluster}}(f) = \varphi_{\text{cluster}}/\sum_{e \in c} I_e\). For a stand-alone electrode this is identical to :meth:grounding_impedance. For connected electrodes it corresponds to the parallel combination of the individual grounding admittances.

Parameters:

Name	Type	Description	Default
electrode_name	str	Name of any electrode in the target cluster. The cluster members are looked up in :attr:clusters; if the electrode does not appear there it is treated as a stand-alone cluster [electrode_name].	required

Returns:

Type	Description
list[complex]	Complex cluster impedance per frequency in ohms, one entry per entry in :attr:frequencies. Frequencies for which the summed cluster current is exactly zero return nan.

Raises:

Type	Description
KeyError	If the resolved cluster is empty.

Source code in src/groundfield/solver/result.py

def cluster_impedance(self, electrode_name: str) -> list[complex]:
    """Grounding impedance of the galvanic cluster containing an electrode.

    Definition: $Z_{\\text{cluster}}(f) =
    \\varphi_{\\text{cluster}}/\\sum_{e \\in c} I_e$. For a
    stand-alone electrode this is identical to
    :meth:`grounding_impedance`. For connected electrodes it
    corresponds to the parallel combination of the individual
    grounding admittances.

    Parameters
    ----------
    electrode_name : str
        Name of any electrode in the target cluster. The cluster
        members are looked up in :attr:`clusters`; if the
        electrode does not appear there it is treated as a
        stand-alone cluster ``[electrode_name]``.

    Returns
    -------
    list[complex]
        Complex cluster impedance per frequency in ohms, one entry
        per entry in :attr:`frequencies`. Frequencies for which
        the summed cluster current is exactly zero return ``nan``.

    Raises
    ------
    KeyError
        If the resolved cluster is empty.
    """
    members = self.clusters.get(electrode_name, [electrode_name])
    if not members:
        raise KeyError(f"No cluster for '{electrode_name}'.")
    u = self.electrode_potentials[members[0]]
    i_sum = [
        sum(self.electrode_currents[m][k] for m in members)
        for k in range(len(u))
    ]
    return [
        (uk / ik) if ik != 0 else complex("nan")
        for uk, ik in zip(u, i_sum)
    ]

grounding_impedance

grounding_impedance(electrode_name: str) -> list[complex]

Input impedance \(Z(f) = U/I\) of an electrode.

Parameters:

Name	Type	Description	Default
electrode_name	str	Name of the electrode whose input impedance is requested. Must be present in both :attr:electrode_potentials and :attr:electrode_currents.	required

Returns:

Type	Description
list[complex]	Complex impedance per frequency in ohms, one entry per entry in :attr:frequencies. Frequencies for which the electrode current is exactly zero return nan.

Raises:

Type	Description
KeyError	If electrode_name is unknown to the result.

Notes

For galvanically connected electrodes (cluster with more than one member) this quantity is the cluster potential divided by the electrode's share of the total current. The physically meaningful quantity is the cluster impedance (see :meth:cluster_impedance).

Source code in src/groundfield/solver/result.py

def grounding_impedance(self, electrode_name: str) -> list[complex]:
    """Input impedance $Z(f) = U/I$ of an electrode.

    Parameters
    ----------
    electrode_name : str
        Name of the electrode whose input impedance is requested.
        Must be present in both :attr:`electrode_potentials` and
        :attr:`electrode_currents`.

    Returns
    -------
    list[complex]
        Complex impedance per frequency in ohms, one entry per
        entry in :attr:`frequencies`. Frequencies for which the
        electrode current is exactly zero return ``nan``.

    Raises
    ------
    KeyError
        If ``electrode_name`` is unknown to the result.

    Notes
    -----
    For galvanically connected electrodes (cluster with more than
    one member) this quantity is the cluster potential divided by
    the electrode's share of the total current. The physically
    meaningful quantity is the **cluster impedance** (see
    :meth:`cluster_impedance`).
    """
    if electrode_name not in self.electrode_potentials:
        raise KeyError(f"No potential data for '{electrode_name}'.")
    if electrode_name not in self.electrode_currents:
        raise KeyError(f"No current data for '{electrode_name}'.")
    u = self.electrode_potentials[electrode_name]
    i = self.electrode_currents[electrode_name]
    return [
        (uk / ik) if ik != 0 else complex("nan") for uk, ik in zip(u, i)
    ]

potential

potential(
    points: np.ndarray,
    frequency_index: int = 0,
    min_distance: float = 0.001,
) -> np.ndarray

Evaluate the potential at field points (image-charge sum).

Picks the appropriate Green's-function kernel automatically:

• homogeneous (:class:HomogeneousSoil or no self.soil set but soil_resistivity available): classic image-charge sum \(1/r + 1/r_{\text{img}}\).
• 2-layer (:class:TwoLayerSoil, or :class:MultiLayerSoil with exactly two layers): Tagg/Sunde series with adaptive truncation (tolerance \(10^{-6}\), at most 100 terms).
• 1-layer (:class:MultiLayerSoil with a single layer): degenerate case, dispatched to the homogeneous kernel.

Parameters:

Name	Type	Description	Default
points	ndarray	Field points, array of shape (M, 3) in metres.	required
frequency_index	int	Index into :attr:frequencies. Default 0.	0
min_distance	float	Numerical cutoff for 1/r singularities, in metres.	0.001

Returns:

Name	Type	Description
phi	(ndarray, shape(M))	Complex potential in V.

Raises:

Type

Description

NotImplementedError

If self.soil is a :class:MultiLayerSoil with three or more layers. The post-solve potential path does not yet carry an n-layer Green's-function kernel. The solve itself is still correct via the cim / mom_sommerfeld / bem backends — only the explicit :meth:potential evaluation is unavailable for n ≥ 3. Cluster impedances, electrode currents and potentials stored in :attr:electrode_potentials / :attr:cluster_impedance remain accessible.

Notes

Prior to 0.2.0 a :class:MultiLayerSoil with n ≥ 3 silently fell through to the homogeneous kernel, returning incorrect potentials, profiles, touch / step voltages and VTK exports. The error is now raised explicitly so the regime is visible rather than silent.

Source code in src/groundfield/solver/result.py

def potential(
    self,
    points: np.ndarray,
    frequency_index: int = 0,
    min_distance: float = 1e-3,
) -> np.ndarray:
    """Evaluate the potential at field points (image-charge sum).

    Picks the appropriate Green's-function kernel automatically:

    - **homogeneous** (:class:`HomogeneousSoil` or no ``self.soil``
      set but ``soil_resistivity`` available): classic image-charge
      sum $1/r + 1/r_{\\text{img}}$.
    - **2-layer** (:class:`TwoLayerSoil`, or :class:`MultiLayerSoil`
      with exactly two layers): Tagg/Sunde series with adaptive
      truncation (tolerance $10^{-6}$, at most 100 terms).
    - **1-layer** (:class:`MultiLayerSoil` with a single layer):
      degenerate case, dispatched to the homogeneous kernel.

    Parameters
    ----------
    points
        Field points, array of shape ``(M, 3)`` in metres.
    frequency_index
        Index into :attr:`frequencies`. Default 0.
    min_distance
        Numerical cutoff for 1/r singularities, in metres.

    Returns
    -------
    phi : np.ndarray, shape (M,)
        Complex potential in V.

    Raises
    ------
    NotImplementedError
        If ``self.soil`` is a :class:`MultiLayerSoil` with three or
        more layers. The post-solve potential path does not yet
        carry an n-layer Green's-function kernel. The solve itself
        is still correct via the ``cim`` / ``mom_sommerfeld`` /
        ``bem`` backends — only the explicit
        :meth:`potential` evaluation is unavailable for n ≥ 3.
        Cluster impedances, electrode currents and potentials
        stored in :attr:`electrode_potentials` /
        :attr:`cluster_impedance` remain accessible.

    Notes
    -----
    Prior to 0.2.0 a :class:`MultiLayerSoil` with n ≥ 3 silently
    fell through to the homogeneous kernel, returning incorrect
    potentials, profiles, touch / step voltages and VTK exports.
    The error is now raised explicitly so the regime is visible
    rather than silent.
    """
    if not self.point_sources:
        raise RuntimeError(
            "FieldResult has no point_sources — the backend did not "
            "populate the post-processing data (stub?)."
        )

    pts = np.asarray(points, dtype=float)
    if pts.ndim == 1:
        pts = pts[None, :]
    if pts.shape[1] != 3:
        raise ValueError(
            f"points must have shape (M, 3), got {pts.shape}."
        )

    sources = np.array([ps.position for ps in self.point_sources])
    currents = np.array(
        [ps.current[frequency_index] for ps in self.point_sources],
        dtype=complex,
    )

    # Dispatch to the kernel that matches the effective number of
    # soil layers. A ``MultiLayerSoil`` with one or two layers is
    # mathematically equivalent to the homogeneous / two-layer
    # case and is cast on the fly. For n ≥ 3 there is no
    # closed-form image-charge series (Γ_1(λ) is no longer
    # constant in λ — see ``solver/image_nlayer.py``), so the
    # call is rejected explicitly rather than silently falling
    # back to the homogeneous kernel.
    if isinstance(self.soil, TwoLayerSoil):
        return self._potential_two_layer(
            pts, sources, currents, self.soil, min_distance
        )
    if isinstance(self.soil, MultiLayerSoil):
        n = len(self.soil.layers)
        if n == 1:
            # Degenerate 1-layer MultiLayerSoil: use the single
            # layer's resistivity as the homogeneous case.
            return self._potential_homogeneous(
                pts,
                sources,
                currents,
                min_distance,
                rho_override=float(self.soil.layers[0].resistivity),
            )
        if n == 2:
            # Two-layer MultiLayerSoil: cast to TwoLayerSoil and
            # reuse the Tagg/Sunde kernel.
            two_layer = TwoLayerSoil(
                rho_1=float(self.soil.layers[0].resistivity),
                rho_2=float(self.soil.layers[1].resistivity),
                h_1=float(self.soil.layers[0].thickness),
            )
            return self._potential_two_layer(
                pts, sources, currents, two_layer, min_distance
            )
        raise NotImplementedError(
            f"FieldResult.potential: MultiLayerSoil with n={n} ≥ 3 "
            "layers is not supported on the post-solve potential "
            "path. The solve itself is correct (use the 'cim', "
            "'mom_sommerfeld', or 'bem' backend), but a closed-form "
            "n-layer Green's-function kernel is not yet wired into "
            "FieldResult.potential. Access cluster impedances and "
            "electrode potentials directly via "
            "result.electrode_potentials / result.cluster_impedance "
            "instead, or reduce the soil model to "
            "TwoLayerSoil / HomogeneousSoil for explicit field "
            "evaluation."
        )
    return self._potential_homogeneous(pts, sources, currents, min_distance)

summary

summary() -> str

Compact one-line description of the result.

Returns:

Type	Description
str	Human-readable summary listing the backend name, the number of frequencies, the number of electrodes for which potentials were stored and the number of discretised point sources. Intended for logging and notebook output; the format is informational and not stable across versions.

Source code in src/groundfield/solver/result.py

def summary(self) -> str:
    """Compact one-line description of the result.

    Returns
    -------
    str
        Human-readable summary listing the backend name, the
        number of frequencies, the number of electrodes for which
        potentials were stored and the number of discretised
        point sources. Intended for logging and notebook output;
        the format is informational and not stable across
        versions.
    """
    return (
        f"FieldResult(backend='{self.backend}', "
        f"n_freq={len(self.frequencies)}, "
        f"n_electrodes={len(self.electrode_potentials)}, "
        f"n_segments={len(self.point_sources)})"
    )

PointSource

Bases: BaseModel

A discretised point current source (segment midpoint).

Filled by the backend; not instantiated directly by the user.

solve_image

solve_image(
    world: "World", engine: "Engine"
) -> FieldResult

Image-charge solver for homogeneous soil.

Parameters:

Name	Type	Description	Default
world	'World'	World whose soil must be a :class:HomogeneousSoil.	required
engine	'Engine'	Engine configuration; relevant fields are segment_length and frequencies.	required

Source code in src/groundfield/solver/image.py

def solve_image(world: "World", engine: "Engine") -> FieldResult:
    """Image-charge solver for homogeneous soil.

    Parameters
    ----------
    world
        World whose ``soil`` must be a :class:`HomogeneousSoil`.
    engine
        Engine configuration; relevant fields are ``segment_length``
        and ``frequencies``.
    """
    if not isinstance(world.soil, HomogeneousSoil):
        raise TypeError(
            "Backend 'image' requires HomogeneousSoil. "
            f"Got: {type(world.soil).__name__}. "
            "For layered models pick backend='image_2layer' "
            "(Tagg/Sunde) or backend='mom' (planned)."
        )
    if not world.electrodes:
        raise ValueError("World contains no electrodes.")

    rho = world.soil.resistivity
    ds = engine.segment_length

    # 1) Discretisation of the electrodes
    all_segments: list[_Segment] = []
    elec_to_segidx: dict[str, list[int]] = {}
    for e in world.electrodes:
        segs = _discretize_electrode(e, ds)
        elec_to_segidx[e.name] = list(range(len(all_segments),
                                            len(all_segments) + len(segs)))
        all_segments.extend(segs)

    # 2) Per-electrode input currents from the configured sources
    elec_input_current: dict[str, complex] = {
        e.name: 0j for e in world.electrodes
    }
    n_ignored_sources = 0
    for src in world.sources:
        if src.kind != "current":
            # Voltage sources are not supported by this backend.
            n_ignored_sources += 1
            continue
        i_complex = src.magnitude * np.exp(1j * np.deg2rad(src.phase_deg))
        if src.attached_to in elec_input_current:
            elec_input_current[src.attached_to] += i_complex
    if n_ignored_sources:
        warnings.warn(
            f"solve_image: {n_ignored_sources} non-current source(s) "
            "ignored — only CurrentSource is supported. A world driven "
            "solely by voltage sources solves to all-zero.",
            UserWarning,
            stacklevel=2,
        )

    # 3) Cluster building: electrodes joined by an *ideal* conductor
    #    share a common potential. Conductors with a finite series
    #    resistance enter the linear system as branches of the
    #    nodal-analysis solver instead.
    cluster_id = _build_clusters(world.electrodes, world.conductors)
    finite_branches = _build_finite_branches(world.conductors, cluster_id)

    # 3b) Distributed-conductor topology (ADR-0003).
    #     Conductors with a finite ``discretize_segment_length`` are
    #     split into sub-pieces; ``coupling_to_soil="galvanic"`` adds
    #     midpoint pseudo-electrode segments to the Z-matrix, while
    #     the longitudinal-segment chain is appended to the branch
    #     list. Pseudo-nodes get one-element entries in
    #     ``elec_to_segidx`` and are flagged as their own clusters.
    cond_segs, distributed_branches_objs, interior_nodes = _build_distributed_topology(
        world.conductors, cluster_id
    )
    pseudo_owners: list[str] = []
    for s in cond_segs:
        pn = s.electrode_name
        elec_to_segidx[pn] = [len(all_segments)]
        all_segments.append(s)
        cluster_id[pn] = pn
        pseudo_owners.append(pn)
    # Anonymous interior nodes from isolated distributed conductors
    # (no leakage segment, only KCL participation) are also flagged
    # as standalone clusters.
    for n in interior_nodes:
        if n not in cluster_id:
            cluster_id[n] = n
            pseudo_owners.append(n)
            elec_to_segidx[n] = []  # no segment, no leakage
    n_lumped_branches = len(finite_branches)
    distributed_branch_tuples = [
        (db.node_a, db.node_b, db.R) for db in distributed_branches_objs
    ]
    finite_branches = list(finite_branches) + distributed_branch_tuples

    # ADR-0004 + ADR-0005: assemble the partial-inductance matrix
    # for the active distributed-conductor branches (lumped branches
    # stay purely resistive and contribute zero entries). When
    # ``engine.earth_inductive_model == "carson_series"`` we also
    # receive a closure that builds dZ_carson(omega) per frequency.
    earth_inductive_model = getattr(
        engine, "earth_inductive_model", "perfect_mirror"
    )
    sigma_earth_for_carson: float | None = None
    layered_earth_for_sommerfeld: object = None
    if earth_inductive_model == "carson_series":
        from groundfield.coupling import resolve_earth_conductivity

        sigma_earth_for_carson = resolve_earth_conductivity(world.soil)
    elif earth_inductive_model == "sommerfeld":
        from groundfield.coupling import resolve_earth_layers

        layered_earth_for_sommerfeld = resolve_earth_layers(world.soil)
    inductance_matrix_full, has_inductance, carson_builder = _assemble_inductance_matrix(
        distributed_branches_objs,
        n_lumped_branches=n_lumped_branches,
        n_total_branches=len(finite_branches),
        earth_model=earth_inductive_model,
        sigma_earth=sigma_earth_for_carson,
        layered_earth=layered_earth_for_sommerfeld,
    )

    n_segments = len(all_segments)
    seg_points = np.array([s.midpoint for s in all_segments])  # (N, 3)
    seg_lengths = np.array([s.length for s in all_segments])    # (N,)

    # 4) Current sharing within clusters via the multi-port matrix
    wire_radii = np.array([s.wire_radius for s in all_segments])
    # ADR-0012 V2: per-segment concrete-shell coefficient (zero for
    # every segment that does not belong to a concrete-encased
    # foundation electrode — historic case).
    seg_shell_coeffs = np.array(
        [s.concrete_shell_coefficient_ohm_m for s in all_segments],
        dtype=float,
    )

    def _homogeneous_self(seg_pts, seg_lens, wr, currents):
        """Closure: homogeneous self-action with fixed rho."""
        return _self_corrected_kernel(seg_pts, seg_lens, wr, currents, rho)

    n_freq = len(engine.frequencies)
    omegas = [2.0 * np.pi * float(f) for f in engine.frequencies]
    real_electrode_names = {e.name for e in world.electrodes}

    # ADR-0010 Tier 1 (WP-F): the multi-port grounding matrix Z is
    # frequency-independent — share it across the per-frequency calls.
    _mp_cache: dict = {}

    def _solve_at(omega: float) -> tuple[dict[str, complex], np.ndarray]:
        """Solve once at a given angular frequency.

        Returns
        -------
        elec_total : dict
            Per-owner total leakage current.
        seg_currents : np.ndarray
            Per-segment current distribution (uniform per unit length).
        """
        carson_dz = (
            carson_builder(omega) if (has_inductance and carson_builder is not None)
            else None
        )
        elec_total = _solve_cluster_currents(
            electrodes=world.electrodes,
            elec_input_current=elec_input_current,
            cluster_id=cluster_id,
            seg_points=seg_points,
            seg_lengths=seg_lengths,
            wire_radii=wire_radii,
            elec_to_segidx=elec_to_segidx,
            self_kernel=_homogeneous_self,
            finite_branches=finite_branches,
            pseudo_owners=pseudo_owners,
            omega=omega if has_inductance else 0.0,
            inductance_matrix=inductance_matrix_full if has_inductance else None,
            carson_correction=carson_dz,
            shell_coefficients=seg_shell_coeffs,
            multiport_cache=_mp_cache,
        )
        sc = np.zeros(n_segments, dtype=complex)
        for ename, idxs in elec_to_segidx.items():
            if not idxs:
                continue
            I_total = elec_total.get(ename, 0j)
            if I_total == 0j:
                continue
            L_total = seg_lengths[idxs].sum()
            sc[idxs] = I_total * seg_lengths[idxs] / L_total
        return elec_total, sc

    def _phi_batch(sc_list: list[np.ndarray]) -> list[np.ndarray]:
        """Segment-midpoint potentials for a list of current vectors.

        All real/imaginary parts are stacked into one
        ``(n_segments, 2·len(sc_list))`` excitation matrix so the
        kernel's O(N²) geometry tensors are built exactly once for
        the whole frequency set (ADR-0010 Tier 1).
        """
        k = len(sc_list)
        stacked = np.zeros((n_segments, 2 * k))
        for m, sc in enumerate(sc_list):
            stacked[:, m] = sc.real
            stacked[:, k + m] = sc.imag
        if not stacked.any():
            return [np.zeros(n_segments, dtype=complex)] * k
        phi = _self_corrected_kernel(
            seg_points, seg_lengths, wire_radii, stacked, rho,
        )
        # ADR-0012 V2: same shell augmentation as in
        # _solve_cluster_currents so electrode_potentials /
        # cluster_impedance reflect the concrete-shell drop.
        if np.any(seg_shell_coeffs > 0.0):
            with np.errstate(divide="ignore", invalid="ignore"):
                shell_diag = np.where(
                    seg_lengths > 0.0,
                    seg_shell_coeffs / seg_lengths,
                    0.0,
                )
            phi = phi + shell_diag[:, None] * stacked
        return [phi[:, m] + 1j * phi[:, k + m] for m in range(k)]

    # Frequency loop. With no inductive coupling the system is
    # frequency-independent, so we solve once and replicate.
    elec_per_freq: list[dict[str, complex]] = []
    sc_per_freq: list[np.ndarray] = []
    phi_per_freq: list[np.ndarray] = []
    if has_inductance:
        for omega in omegas:
            et, sc = _solve_at(omega)
            elec_per_freq.append(et)
            sc_per_freq.append(sc)
        phi_per_freq = _phi_batch(sc_per_freq)
    else:
        et, sc = _solve_at(0.0)
        ph = _phi_batch([sc])[0]
        elec_per_freq = [et] * n_freq
        sc_per_freq = [sc] * n_freq
        phi_per_freq = [ph] * n_freq

    # Build the FieldResult mappings.
    electrode_potentials: dict[str, list[complex]] = {}
    electrode_currents: dict[str, list[complex]] = {}
    conductor_currents: dict[str, list[complex]] = {}
    conductor_potentials: dict[str, list[complex]] = {}
    for ename, idxs in elec_to_segidx.items():
        if not idxs:
            continue
        u_list = [
            complex(np.mean(phi_per_freq[k][idxs])) for k in range(n_freq)
        ]
        i_list = [elec_per_freq[k][ename] for k in range(n_freq)]
        if ename in real_electrode_names:
            electrode_potentials[ename] = u_list
            electrode_currents[ename] = i_list
        else:
            conductor_potentials[ename] = u_list
            conductor_currents[ename] = i_list

    # 7) Point-source list for post-processing (plots, profiles)
    point_sources = [
        PointSource(
            position=tuple(seg_points[i].tolist()),
            current=[complex(sc_per_freq[k][i]) for k in range(n_freq)],
            electrode_name=all_segments[i].electrode_name,
            length=float(seg_lengths[i]),
        )
        for i in range(n_segments)
    ]

    # Cluster map: electrode_name -> sorted list of cluster members
    # (only real electrodes are surfaced).
    cluster_members: dict[str, list[str]] = {}
    for ename in real_electrode_names:
        cluster_members[ename] = sorted(
            n for n in cluster_id
            if cluster_id[n] == cluster_id[ename] and n in real_electrode_names
        )

    metadata: dict = {
        "world_name": world.name,
        "n_segments": n_segments,
        "segment_length": ds,
        "stub": False,
        "earth_inductive_model": earth_inductive_model,
    }
    # ADR-0005 §"Eindringtiefen-Diagnostik": expose the
    # electromagnetic skin depth in soil at every solved frequency,
    # so notebooks and benchmarks can answer "is my geometry small
    # or large compared to delta(omega)?" without re-deriving the
    # formula. Only active for engines that ran a frequency loop.
    if has_inductance and sigma_earth_for_carson is not None:
        from groundfield.coupling.carson import skin_depth

        metadata["penetration_depth"] = {
            float(f): skin_depth(2.0 * np.pi * f, sigma_earth_for_carson)
            for f in engine.frequencies
        }
    elif has_inductance and isinstance(world.soil, HomogeneousSoil):
        # No Carson active, but homogeneous soil — still useful as a
        # *reference* skin depth, even though the perfect-mirror
        # solver does not actually use it.
        from groundfield.coupling.carson import skin_depth

        sigma_ref = 1.0 / float(world.soil.resistivity)
        metadata["penetration_depth"] = {
            float(f): skin_depth(2.0 * np.pi * f, sigma_ref)
            for f in engine.frequencies
        }
    if conductor_currents:
        metadata["conductor_node_currents"] = conductor_currents
        metadata["conductor_node_potentials"] = conductor_potentials

    return FieldResult(
        backend="image",
        frequencies=list(engine.frequencies),
        electrode_potentials=electrode_potentials,
        electrode_currents=electrode_currents,
        point_sources=point_sources,
        soil_resistivity=float(rho),
        soil=world.soil,
        clusters=cluster_members,
        metadata=metadata,
    )

solve_image_2layer

solve_image_2layer(
    world: "World",
    engine: "Engine",
    *,
    max_terms: int = 100,
    tol: float = 1e-06
) -> FieldResult

Image-charge solver for 2-layer soil (Tagg / Sunde).

Parameters:

Name	Type	Description	Default
world	'World'	World whose soil is a :class:TwoLayerSoil. All electrodes must lie inside the upper layer.	required
engine	'Engine'	Engine configuration. engine.segment_length controls the discretisation as in the homogeneous backend.	required
max_terms	int	Maximum number of series terms.	100
tol	float	Series truncation: stop as soon as \(|K|^n < \text{tol}\).	1e-06

Returns:

Type	Description
FieldResult	Result object. :attr:FieldResult.metadata exposes diagnostic fields (backend, K, n_terms_used, converged).

Source code in src/groundfield/solver/image_2layer.py

def solve_image_2layer(
    world: "World",
    engine: "Engine",
    *,
    max_terms: int = 100,
    tol: float = 1e-6,
) -> FieldResult:
    """Image-charge solver for 2-layer soil (Tagg / Sunde).

    Parameters
    ----------
    world
        World whose ``soil`` is a :class:`TwoLayerSoil`. All electrodes
        must lie inside the upper layer.
    engine
        Engine configuration. ``engine.segment_length`` controls the
        discretisation as in the homogeneous backend.
    max_terms
        Maximum number of series terms.
    tol
        Series truncation: stop as soon as $|K|^n < \\text{tol}$.

    Returns
    -------
    FieldResult
        Result object. :attr:`FieldResult.metadata` exposes diagnostic
        fields (``backend``, ``K``, ``n_terms_used``, ``converged``).
    """
    if not isinstance(world.soil, TwoLayerSoil):
        raise TypeError(
            "Backend 'image_2layer' requires TwoLayerSoil. "
            f"Got: {type(world.soil).__name__}."
        )
    if not world.electrodes:
        raise ValueError("World contains no electrodes.")

    soil: TwoLayerSoil = world.soil
    h_1 = soil.h_1
    K = soil.reflection_coefficient
    ds = engine.segment_length

    _log.info(
        "image_2layer: rho_1=%.1f, rho_2=%.1f, h_1=%.2f, K=%.4f",
        soil.rho_1, soil.rho_2, h_1, K,
    )

    # 1) Discretisation — identical to the homogeneous backend, but
    #    with segments split at the layer interface (ADR-0007 step 1,
    #    audit 2026-07-08 WP-D2) so no segment straddles h_1.
    all_segments: list[_Segment] = []
    elec_to_segidx: dict[str, list[int]] = {}
    for e in world.electrodes:
        segs = _discretize_electrode(e, ds, layer_interfaces=(h_1,))
        elec_to_segidx[e.name] = list(range(len(all_segments),
                                            len(all_segments) + len(segs)))
        all_segments.extend(segs)

    # 2) Per-electrode input currents from sources
    elec_input_current: dict[str, complex] = {
        e.name: 0j for e in world.electrodes
    }
    n_ignored_sources = 0
    for src in world.sources:
        if src.kind != "current":
            n_ignored_sources += 1
            continue
        i_complex = src.magnitude * np.exp(1j * np.deg2rad(src.phase_deg))
        if src.attached_to in elec_input_current:
            elec_input_current[src.attached_to] += i_complex
    if n_ignored_sources:
        warnings.warn(
            f"image_2layer: {n_ignored_sources} non-current source(s) "
            "ignored — only CurrentSource is supported.",
            UserWarning,
            stacklevel=2,
        )

    # 3) Cluster building (ideal conductors only) and finite-impedance
    #    branch list (passed into the nodal-analysis solver).
    cluster_id = _build_clusters(world.electrodes, world.conductors)
    finite_branches = _build_finite_branches(world.conductors, cluster_id)

    # 3b) Distributed-conductor topology (ADR-0003) + ADR-0004
    #     inductive coupling assembly.
    cond_segs, distributed_branches_objs, interior_nodes = _build_distributed_topology(
        world.conductors, cluster_id
    )
    pseudo_owners: list[str] = []
    for s in cond_segs:
        pn = s.electrode_name
        elec_to_segidx[pn] = [len(all_segments)]
        all_segments.append(s)
        cluster_id[pn] = pn
        pseudo_owners.append(pn)
    for n_ in interior_nodes:
        if n_ not in cluster_id:
            cluster_id[n_] = n_
            pseudo_owners.append(n_)
            elec_to_segidx[n_] = []
    n_lumped_branches = len(finite_branches)
    distributed_branch_tuples = [
        (db.node_a, db.node_b, db.R) for db in distributed_branches_objs
    ]
    finite_branches = list(finite_branches) + distributed_branch_tuples
    earth_inductive_model = getattr(
        engine, "earth_inductive_model", "perfect_mirror"
    )
    sigma_earth_for_carson: float | None = None
    layered_earth_for_sommerfeld: object = None
    if earth_inductive_model == "carson_series":
        from groundfield.coupling import resolve_earth_conductivity

        sigma_earth_for_carson = resolve_earth_conductivity(world.soil)
    elif earth_inductive_model == "sommerfeld":
        from groundfield.coupling import resolve_earth_layers

        layered_earth_for_sommerfeld = resolve_earth_layers(world.soil)
    inductance_matrix_full, has_inductance, carson_builder = _assemble_inductance_matrix(
        distributed_branches_objs,
        n_lumped_branches=n_lumped_branches,
        n_total_branches=len(finite_branches),
        earth_model=earth_inductive_model,
        sigma_earth=sigma_earth_for_carson,
        layered_earth=layered_earth_for_sommerfeld,
    )

    seg_points = np.array([s.midpoint for s in all_segments])
    seg_lengths = np.array([s.length for s in all_segments])
    wire_radii = np.array([s.wire_radius for s in all_segments])
    # ADR-0012 V2: per-segment radial shell coefficient. Zero for
    # any segment that does not belong to a concrete-encased
    # foundation electrode — historic case preserved bit-exact.
    seg_shell_coeffs = np.array(
        [s.concrete_shell_coefficient_ohm_m for s in all_segments],
        dtype=float,
    )

    # ADR-0007: detect cross-layer geometry. When all segments live
    # inside the upper layer, the historic Tagg/Sunde image series
    # is used (bit-exact regression). When any segment is at or
    # below the interface, the kernel automatically dispatches to
    # the rigorous Sommerfeld path.
    z_max = seg_points[:, 2].max()
    cross_layer = bool(z_max >= h_1)
    if cross_layer:
        _log.info(
            "image_2layer: cross-layer geometry detected "
            "(z_max = %.3f m, h_1 = %.3f m). Using ADR-0007 "
            "Sommerfeld fallback for all source/observer pairs.",
            z_max, h_1,
        )
    elif h_1 - z_max < 2.0 * ds:
        # Audit 2026-07-08 WP-D2: within ~one segment of the
        # interface the point-source treatment of the K-weighted
        # nearest image (distance 2(h_1 − z)) degrades. Surface the
        # margin so the user can refine ds or accept the bias.
        warnings.warn(
            f"image_2layer: deepest segment ends {h_1 - z_max:.2f} m "
            f"above the layer interface (h_1 = {h_1:.2f} m), which is "
            f"less than 2×segment_length ({2.0 * ds:.2f} m). The "
            "image-series self-terms near the interface carry an "
            "elevated discretisation error there — consider a finer "
            "segment_length.",
            UserWarning,
            stacklevel=2,
        )
    # Tag every segment with its layer index (0 = upper, 1 = lower).
    for s in all_segments:
        s.layer_index = 0 if s.midpoint[2] < h_1 else 1

    # 4) Self-kernel closure and frequency loop
    self_kernel = _two_layer_self_kernel_factory(
        soil, max_terms, tol, allow_cross_layer=cross_layer,
    )
    n_segments = len(all_segments)
    n_freq = len(engine.frequencies)
    omegas = [2.0 * np.pi * float(f) for f in engine.frequencies]
    real_electrode_names = {e.name for e in world.electrodes}

    # ADR-0010 Tier 1 (WP-F): the multi-port grounding matrix Z is
    # frequency-independent — share it across the per-frequency calls.
    _mp_cache: dict = {}

    def _solve_at(omega: float) -> tuple[dict[str, complex], np.ndarray]:
        carson_dz = (
            carson_builder(omega) if (has_inductance and carson_builder is not None)
            else None
        )
        elec_total = _solve_cluster_currents(
            electrodes=world.electrodes,
            elec_input_current=elec_input_current,
            cluster_id=cluster_id,
            seg_points=seg_points,
            seg_lengths=seg_lengths,
            wire_radii=wire_radii,
            elec_to_segidx=elec_to_segidx,
            self_kernel=self_kernel,
            finite_branches=finite_branches,
            pseudo_owners=pseudo_owners,
            omega=omega if has_inductance else 0.0,
            inductance_matrix=inductance_matrix_full if has_inductance else None,
            carson_correction=carson_dz,
            shell_coefficients=seg_shell_coeffs,
            multiport_cache=_mp_cache,
        )
        sc = np.zeros(n_segments, dtype=complex)
        for ename, idxs in elec_to_segidx.items():
            if not idxs:
                continue
            I_total = elec_total.get(ename, 0j)
            if I_total == 0j:
                continue
            L_total = seg_lengths[idxs].sum()
            sc[idxs] = I_total * seg_lengths[idxs] / L_total
        return elec_total, sc

    def _phi_batch(sc_list: list[np.ndarray]) -> list[np.ndarray]:
        """Segment-midpoint potentials for a list of current vectors.

        Stacks all real/imaginary parts into one excitation matrix so
        the kernel's O(N²·M_series) tensors are built exactly once
        for the whole frequency set (ADR-0010 Tier 1).
        """
        k = len(sc_list)
        stacked = np.zeros((n_segments, 2 * k))
        for m, sc in enumerate(sc_list):
            stacked[:, m] = sc.real
            stacked[:, k + m] = sc.imag
        if not stacked.any():
            return [np.zeros(n_segments, dtype=complex)] * k
        phi = self_kernel(seg_points, seg_lengths, wire_radii, stacked)
        # ADR-0012 V2: post-solve potential evaluation must carry the
        # same shell augmentation as the inverse-problem assembly
        # inside _solve_cluster_currents. Without this step the
        # electrode_potentials and cluster_impedance ignore the shell
        # drop and report the bare-soil values, hiding the V2 effect.
        if np.any(seg_shell_coeffs > 0.0):
            with np.errstate(divide="ignore", invalid="ignore"):
                shell_diag = np.where(
                    seg_lengths > 0.0,
                    seg_shell_coeffs / seg_lengths,
                    0.0,
                )
            phi = phi + shell_diag[:, None] * stacked
        return [phi[:, m] + 1j * phi[:, k + m] for m in range(k)]

    elec_per_freq: list[dict[str, complex]] = []
    sc_per_freq: list[np.ndarray] = []
    phi_per_freq: list[np.ndarray] = []
    if has_inductance:
        for omega in omegas:
            et, sc = _solve_at(omega)
            elec_per_freq.append(et)
            sc_per_freq.append(sc)
        phi_per_freq = _phi_batch(sc_per_freq)
    else:
        et, sc = _solve_at(0.0)
        ph = _phi_batch([sc])[0]
        elec_per_freq = [et] * n_freq
        sc_per_freq = [sc] * n_freq
        phi_per_freq = [ph] * n_freq

    n_terms_used_self = 0
    if any(np.any(sc) for sc in sc_per_freq):
        _, _, _, n_terms_used_self = _two_layer_image_offsets(
            K, h_1, max_terms, tol
        )

    electrode_potentials: dict[str, list[complex]] = {}
    electrode_currents: dict[str, list[complex]] = {}
    conductor_currents: dict[str, list[complex]] = {}
    conductor_potentials: dict[str, list[complex]] = {}
    for ename, idxs in elec_to_segidx.items():
        if not idxs:
            continue
        u_list = [
            complex(np.mean(phi_per_freq[k][idxs])) for k in range(n_freq)
        ]
        i_list = [elec_per_freq[k][ename] for k in range(n_freq)]
        if ename in real_electrode_names:
            electrode_potentials[ename] = u_list
            electrode_currents[ename] = i_list
        else:
            conductor_potentials[ename] = u_list
            conductor_currents[ename] = i_list

    # 7) Point-source list for post-processing
    point_sources = [
        PointSource(
            position=tuple(seg_points[i].tolist()),
            current=[complex(sc_per_freq[k][i]) for k in range(n_freq)],
            electrode_name=all_segments[i].electrode_name,
            length=float(seg_lengths[i]),
        )
        for i in range(n_segments)
    ]

    cluster_members: dict[str, list[str]] = {}
    for ename in real_electrode_names:
        cluster_members[ename] = sorted(
            n for n in cluster_id
            if cluster_id[n] == cluster_id[ename] and n in real_electrode_names
        )

    converged = (
        _series_tail_bound(abs(K), n_terms_used_self) < tol
        if n_terms_used_self else True
    )
    if not converged:
        tail = _series_tail_bound(abs(K), n_terms_used_self)
        _log.warning(
            "image_2layer: max_terms=%d reached, series tail bound "
            "|K|^(n+1)/(1-|K|) = %.2e > tol=%.2e. Result may be "
            "inaccurate.",
            max_terms, tail, tol,
        )
        warnings.warn(
            f"image_2layer: the Tagg/Sunde image series was truncated "
            f"at max_terms={max_terms} with a remaining tail bound of "
            f"{tail:.2e} > tol={tol:.2e} (|K| = {abs(K):.3f}). Raise "
            "Engine.image_max_terms for high layer contrasts.",
            SeriesTruncationWarning,
            stacklevel=2,
        )

    metadata = {
        "world_name": world.name,
        "n_segments": n_segments,
        "segment_length": ds,
        "K": float(K),
        "rho_1": float(soil.rho_1),
        "rho_2": float(soil.rho_2),
        "h_1": float(h_1),
        "n_terms_used": n_terms_used_self,
        "converged": converged,
        "stub": False,
        "earth_inductive_model": earth_inductive_model,
    }
    if has_inductance and sigma_earth_for_carson is not None:
        from groundfield.coupling.carson import skin_depth

        metadata["penetration_depth"] = {
            float(f): skin_depth(2.0 * np.pi * f, sigma_earth_for_carson)
            for f in engine.frequencies
        }
    elif has_inductance:
        from groundfield.coupling.carson import skin_depth

        sigma_ref = 1.0 / float(soil.rho_1)
        metadata["penetration_depth"] = {
            float(f): skin_depth(2.0 * np.pi * f, sigma_ref)
            for f in engine.frequencies
        }
    if conductor_currents:
        metadata["conductor_node_currents"] = conductor_currents
        metadata["conductor_node_potentials"] = conductor_potentials

    return FieldResult(
        backend="image_2layer",
        frequencies=list(engine.frequencies),
        electrode_potentials=electrode_potentials,
        electrode_currents=electrode_currents,
        point_sources=point_sources,
        soil_resistivity=float(soil.rho_1),
        soil=soil,
        clusters=cluster_members,
        metadata=metadata,
    )

solve_mom

solve_mom(
    world: "World",
    engine: "Engine",
    *,
    two_layer_max_terms: int | None = None,
    two_layer_tol: float | None = None
) -> FieldResult

Galerkin Method-of-Moments backend.

Supports :class:HomogeneousSoil and :class:TwoLayerSoil. The Green's-function kernel is taken from the matching image backend; only the resolution scheme (per-segment Galerkin instead of per-electrode average potential) differs.

Parameters:

Name	Type	Description	Default
world	'World'	World to evaluate.	required
engine	'Engine'	Engine configuration; engine.segment_length controls the discretisation.	required
two_layer_max_terms	int | None	Truncation parameters for the Tagg/Sunde series, only used when world.soil is a :class:TwoLayerSoil. None (default) takes engine.image_max_terms / engine.image_series_tol, so the engine knob reaches this backend consistently with image_2layer.	None
two_layer_tol	int | None	Truncation parameters for the Tagg/Sunde series, only used when world.soil is a :class:TwoLayerSoil. None (default) takes engine.image_max_terms / engine.image_series_tol, so the engine knob reaches this backend consistently with image_2layer.	None

Returns:

Type	Description
FieldResult	Result object compatible with :class:FieldResult from the image backends. metadata['solver'] is set to 'galerkin'.

Source code in src/groundfield/solver/mom.py

def solve_mom(
    world: "World",
    engine: "Engine",
    *,
    two_layer_max_terms: int | None = None,
    two_layer_tol: float | None = None,
) -> FieldResult:
    """Galerkin Method-of-Moments backend.

    Supports :class:`HomogeneousSoil` and :class:`TwoLayerSoil`. The
    Green's-function kernel is taken from the matching image backend;
    only the resolution scheme (per-segment Galerkin instead of
    per-electrode average potential) differs.

    Parameters
    ----------
    world
        World to evaluate.
    engine
        Engine configuration; ``engine.segment_length`` controls the
        discretisation.
    two_layer_max_terms, two_layer_tol
        Truncation parameters for the Tagg/Sunde series, only used
        when ``world.soil`` is a :class:`TwoLayerSoil`. ``None``
        (default) takes ``engine.image_max_terms`` /
        ``engine.image_series_tol``, so the engine knob reaches this
        backend consistently with ``image_2layer``.

    Returns
    -------
    FieldResult
        Result object compatible with :class:`FieldResult` from the
        image backends. ``metadata['solver']`` is set to ``'galerkin'``.
    """
    if not world.electrodes:
        raise ValueError("World contains no electrodes.")
    if not isinstance(world.soil, (HomogeneousSoil, TwoLayerSoil)):
        raise TypeError(
            "Backend 'mom' currently supports HomogeneousSoil and "
            f"TwoLayerSoil only. Got: {type(world.soil).__name__}."
        )
    _reject_concrete_shells(world, "mom")
    _warn_ignored_sources(world, "mom")
    if two_layer_max_terms is None:
        two_layer_max_terms = engine.image_max_terms
    if two_layer_tol is None:
        two_layer_tol = engine.image_series_tol

    ds = engine.segment_length
    _log.info(
        "mom: soil=%s, segment_length=%.3f",
        type(world.soil).__name__, ds,
    )

    # 1) Discretisation — identical to the image backends; for
    #    two-layer soil segments are split at the interface
    #    (ADR-0007 step 1, audit 2026-07-08 WP-D2).
    interfaces = (
        (world.soil.h_1,) if isinstance(world.soil, TwoLayerSoil) else None
    )
    all_segments: list[_Segment] = []
    elec_to_segidx: dict[str, list[int]] = {}
    for e in world.electrodes:
        segs = _discretize_electrode(e, ds, layer_interfaces=interfaces)
        elec_to_segidx[e.name] = list(range(len(all_segments),
                                            len(all_segments) + len(segs)))
        all_segments.extend(segs)

    # 2) Per-electrode input currents from the configured sources.
    elec_input_current: dict[str, complex] = {
        e.name: 0j for e in world.electrodes
    }
    for src in world.sources:
        if src.kind != "current":
            continue
        i_complex = src.magnitude * np.exp(1j * np.deg2rad(src.phase_deg))
        if src.attached_to in elec_input_current:
            elec_input_current[src.attached_to] += i_complex

    # 3) Cluster building (ideal conductors only), finite-impedance
    #    branch list, plus the distributed-conductor topology
    #    (ADR-0003) — all fed into the augmented Galerkin system.
    cluster_id = _build_clusters(world.electrodes, world.conductors)
    finite_branches = _build_finite_branches(world.conductors, cluster_id)
    cond_segs, distributed_branches_objs, interior_nodes = _build_distributed_topology(
        world.conductors, cluster_id
    )
    for s in cond_segs:
        pn = s.electrode_name
        elec_to_segidx[pn] = [len(all_segments)]
        all_segments.append(s)
        cluster_id[pn] = pn
    for n_ in interior_nodes:
        if n_ not in cluster_id:
            cluster_id[n_] = n_
            elec_to_segidx[n_] = []
    n_lumped_branches = len(finite_branches)
    distributed_branch_tuples = [
        (db.node_a, db.node_b, db.R) for db in distributed_branches_objs
    ]
    finite_branches = list(finite_branches) + distributed_branch_tuples
    earth_inductive_model = getattr(
        engine, "earth_inductive_model", "perfect_mirror"
    )
    sigma_earth_for_carson: float | None = None
    layered_earth_for_sommerfeld: object = None
    if earth_inductive_model == "carson_series":
        from groundfield.coupling import resolve_earth_conductivity

        sigma_earth_for_carson = resolve_earth_conductivity(world.soil)
    elif earth_inductive_model == "sommerfeld":
        from groundfield.coupling import resolve_earth_layers

        layered_earth_for_sommerfeld = resolve_earth_layers(world.soil)
    inductance_matrix_full, has_inductance, carson_builder = _assemble_inductance_matrix(
        distributed_branches_objs,
        n_lumped_branches=n_lumped_branches,
        n_total_branches=len(finite_branches),
        earth_model=earth_inductive_model,
        sigma_earth=sigma_earth_for_carson,
        layered_earth=layered_earth_for_sommerfeld,
    )

    n_segments = len(all_segments)
    seg_points = np.array([s.midpoint for s in all_segments])
    seg_lengths = np.array([s.length for s in all_segments])
    wire_radii = np.array([s.wire_radius for s in all_segments])

    # 4) Cross-layer support (audit 2026-07-08, WP-D3): the shared
    #    self-kernel factory dispatches to the rigorous ADR-0007
    #    Sommerfeld path when any segment sits at or below h_1, so
    #    the historic hard rejection is gone. mom is the engine that
    #    *solves* the per-segment current distribution — exactly the
    #    tool to quantify the image family's uniform-leakage bias
    #    across the interface.
    if isinstance(world.soil, TwoLayerSoil):
        z_max = seg_points[:, 2].max()
        if z_max >= world.soil.h_1:
            _log.info(
                "mom: cross-layer geometry (z_max = %.3f m >= h_1 = "
                "%.3f m) — using the ADR-0007 Sommerfeld kernel.",
                z_max, world.soil.h_1,
            )

    # 5) Assemble the reaction matrix Z (kernel depends on soil model).
    n_terms_used = 0
    if isinstance(world.soil, HomogeneousSoil):
        Z = _build_Z_homogeneous(
            seg_points, seg_lengths, wire_radii, world.soil.resistivity
        )
    else:
        # TwoLayerSoil
        Z = _build_Z_two_layer(
            seg_points, seg_lengths, wire_radii,
            world.soil, two_layer_max_terms, two_layer_tol,
        )
        # Approximate the term count from the same offset generator.
        from groundfield.solver.image_2layer import _two_layer_image_offsets

        _, _, _, n_terms_used = _two_layer_image_offsets(
            world.soil.reflection_coefficient, world.soil.h_1,
            two_layer_max_terms, two_layer_tol,
        )

    # 6) Frequency loop. With no inductive coupling, solve once and
    #    replicate across frequencies; otherwise loop.
    n_freq = len(engine.frequencies)
    omegas = [2.0 * np.pi * float(f) for f in engine.frequencies]
    real_electrode_names = {e.name for e in world.electrodes}

    if isinstance(world.soil, HomogeneousSoil):
        rho_self = world.soil.resistivity
        phi_kernel = lambda sc_real: _self_corrected_kernel(
            seg_points, seg_lengths, wire_radii, sc_real, rho_self,
        )
    else:
        self_kernel = _two_layer_self_kernel_factory(
            world.soil, two_layer_max_terms, two_layer_tol,
            allow_cross_layer=True,
        )
        phi_kernel = lambda sc_real: self_kernel(
            seg_points, seg_lengths, wire_radii, sc_real,
        )

    def _solve_at(omega: float) -> tuple[np.ndarray, np.ndarray]:
        carson_dz = (
            carson_builder(omega) if (has_inductance and carson_builder is not None)
            else None
        )
        sc, _ = _galerkin_solve(
            Z=Z,
            elec_input_current=elec_input_current,
            cluster_id=cluster_id,
            elec_to_segidx=elec_to_segidx,
            n_segments=n_segments,
            finite_branches=finite_branches,
            omega=omega if has_inductance else 0.0,
            inductance_matrix=inductance_matrix_full if has_inductance else None,
            carson_correction=carson_dz,
        )
        ph = np.zeros(n_segments, dtype=complex)
        if sc.any():
            ph = phi_kernel(sc.real) + 1j * phi_kernel(sc.imag)
        return sc, ph

    sc_per_freq: list[np.ndarray] = []
    phi_per_freq: list[np.ndarray] = []
    if has_inductance:
        for omega in omegas:
            sc, ph = _solve_at(omega)
            sc_per_freq.append(sc)
            phi_per_freq.append(ph)
    else:
        sc, ph = _solve_at(0.0)
        sc_per_freq = [sc] * n_freq
        phi_per_freq = [ph] * n_freq

    electrode_potentials: dict[str, list[complex]] = {}
    electrode_currents: dict[str, list[complex]] = {}
    conductor_currents: dict[str, list[complex]] = {}
    conductor_potentials: dict[str, list[complex]] = {}
    for ename, idxs in elec_to_segidx.items():
        if not idxs:
            continue
        u_list = [
            complex(np.mean(phi_per_freq[k][idxs])) for k in range(n_freq)
        ]
        i_list = [
            complex(sc_per_freq[k][idxs].sum()) for k in range(n_freq)
        ]
        if ename in real_electrode_names:
            electrode_potentials[ename] = u_list
            electrode_currents[ename] = i_list
        else:
            conductor_potentials[ename] = u_list
            conductor_currents[ename] = i_list

    # 8) Point-source list and cluster map for post-processing.
    point_sources = [
        PointSource(
            position=tuple(seg_points[i].tolist()),
            current=[complex(sc_per_freq[k][i]) for k in range(n_freq)],
            electrode_name=all_segments[i].electrode_name,
            length=float(seg_lengths[i]),
        )
        for i in range(n_segments)
    ]
    cluster_members: dict[str, list[str]] = {}
    for ename in real_electrode_names:
        cluster_members[ename] = sorted(
            n for n in cluster_id
            if cluster_id[n] == cluster_id[ename] and n in real_electrode_names
        )

    # 9) Soil-specific metadata.
    if isinstance(world.soil, HomogeneousSoil):
        soil_resistivity = float(world.soil.resistivity)
        soil_meta: dict[str, float | int | bool] = {}
    else:
        K = world.soil.reflection_coefficient
        soil_resistivity = float(world.soil.rho_1)
        converged = (
            _series_tail_bound(abs(K), n_terms_used) < two_layer_tol
            if n_terms_used else True
        )
        if not converged:
            warnings.warn(
                f"mom: the Tagg/Sunde image series was truncated at "
                f"max_terms={two_layer_max_terms} with tail bound > "
                f"tol={two_layer_tol:.2e} (|K| = {abs(K):.3f}). Raise "
                "Engine.image_max_terms for high layer contrasts.",
                SeriesTruncationWarning,
                stacklevel=2,
            )
        soil_meta = {
            "K": float(K),
            "rho_1": float(world.soil.rho_1),
            "rho_2": float(world.soil.rho_2),
            "h_1": float(world.soil.h_1),
            "n_terms_used": n_terms_used,
            "converged": bool(converged),
        }

    return FieldResult(
        backend="mom",
        frequencies=list(engine.frequencies),
        electrode_potentials=electrode_potentials,
        electrode_currents=electrode_currents,
        point_sources=point_sources,
        soil_resistivity=soil_resistivity,
        soil=world.soil,
        clusters=cluster_members,
        metadata={
            "world_name": world.name,
            "n_segments": n_segments,
            "segment_length": ds,
            "solver": "galerkin",
            "stub": False,
            "earth_inductive_model": earth_inductive_model,
            **soil_meta,
            **_carson_skin_metadata(
                has_inductance, sigma_earth_for_carson, world.soil,
                engine.frequencies,
            ),
            **(
                {
                    "conductor_node_currents": conductor_currents,
                    "conductor_node_potentials": conductor_potentials,
                }
                if conductor_currents else {}
            ),
        },
    )

solve_mutual_field

solve_mutual_field(
    world: "World",
    engine: "Engine",
    anchors: list[str],
    field_points: np.ndarray,
    *,
    frequency_index: int = 0,
    matrix_min_distance: float = 0.001,
    matrix_max_terms: int = 100,
    matrix_tol: float = 1e-06
) -> np.ndarray

Galvanic Green-function matrix: potential at arbitrary field points per node.

Returns R of shape (M, nG) with R[i, j] = potential at field_points[i] (x, y, depth) under unit current (1 A) injected into anchors[j] while all other clusters stay floating. This is the off-diagonal field evaluation of :func:solve_mutual_matrix, generalised to an arbitrary number of field points (decoupled from nG) and without the self/diagonal term — the building block for the frequency-dependent surface-potential distribution by superposition::

phi(field_points, f) = R @ I_leak(f)

with I_leak(f) = Y_G @ u(f) from the reduced network (Y_L(f) + Y_G) u = i. Uses the same one-shot assembly as :func:solve_mutual_matrix (reaction matrix and per-excitation leakage columns built once); only the field-point evaluation is generalised. The expensive 3D part is therefore frequency-independent and computed once, so a frequency sweep over I_leak(f) is cheap.

world must already carry all nG grounding clusters (every anchor in anchors must be a real electrode in world) and the soil model. Any current source is ignored — the excitations are applied internally, one cluster at a time, exactly as the historic per-excitation loop did.

Parameters:

Name	Type	Description	Default
world	'World'	Assembled world holding all nG grounding clusters and the soil model. Sources, if any, are ignored.	required
engine	'Engine'	Engine configuration. segment_length controls the discretisation; frequencies selects the solved frequency via frequency_index; image_max_terms / image_series_tol feed the 2-layer self-kernel (the reaction-matrix assembly), exactly as Engine.solve would.	required
anchors	list[str]	Length-nG list of electrode names, one per node. Anchor j receives the unit current of excitation j and fills column j of R.	required
field_points	ndarray	(M, 3) array of (x, y, depth) evaluation points, decoupled from nG: M may differ from nG and the points are arbitrary (surface or buried). Row i becomes row i of R.	required
frequency_index	int	Index into engine.frequencies. Default 0. Only relevant when distributed conductors make the assembled system frequency-dependent.	0
matrix_min_distance	float	Truncation / clamp parameters of the field evaluation. The defaults 1e-3 / 100 / 1e-6 reproduce the hard-coded defaults of :meth:FieldResult.potential / _potential_two_layer exactly. (The reaction-matrix assembly uses engine.image_max_terms / image_series_tol independently, just like the live solver.)	0.001
matrix_max_terms	float	Truncation / clamp parameters of the field evaluation. The defaults 1e-3 / 100 / 1e-6 reproduce the hard-coded defaults of :meth:FieldResult.potential / _potential_two_layer exactly. (The reaction-matrix assembly uses engine.image_max_terms / image_series_tol independently, just like the live solver.)	0.001
matrix_tol	float	Truncation / clamp parameters of the field evaluation. The defaults 1e-3 / 100 / 1e-6 reproduce the hard-coded defaults of :meth:FieldResult.potential / _potential_two_layer exactly. (The reaction-matrix assembly uses engine.image_max_terms / image_series_tol independently, just like the live solver.)	0.001

Returns:

Name	Type	Description
R	(ndarray, shape(M, nG), complex)	Galvanic Green-function matrix; R[i, j] is the potential at field_points[i] under unit excitation of cluster j.

Raises:

Type	Description
ValueError	If world has no soil model, holds no electrodes, or field_points does not have shape (M, 3).
KeyError	If an entry of anchors is not an electrode in world.
NotImplementedError	If the soil is a MultiLayerSoil with three or more layers (only homogeneous and 2-layer soils are supported).
TypeError	If the soil model type is not supported.

See Also

solve_mutual_matrix : Full nG × nG mutual grounding-impedance matrix.

Source code in src/groundfield/solver/mutual.py

def solve_mutual_field(
    world: "World",
    engine: "Engine",
    anchors: list[str],
    field_points: np.ndarray,
    *,
    frequency_index: int = 0,
    matrix_min_distance: float = 1e-3,
    matrix_max_terms: int = 100,
    matrix_tol: float = 1e-6,
) -> np.ndarray:
    """Galvanic Green-function matrix: potential at arbitrary field points per node.

    Returns ``R`` of shape ``(M, nG)`` with ``R[i, j]`` = potential at
    ``field_points[i]`` (``x, y, depth``) under unit current (1 A) injected
    into ``anchors[j]`` while all other clusters stay floating. This is the
    off-diagonal field evaluation of :func:`solve_mutual_matrix`, generalised
    to an **arbitrary** number of field points (decoupled from ``nG``) and
    without the self/diagonal term — the building block for the
    frequency-dependent surface-potential distribution by superposition::

        phi(field_points, f) = R @ I_leak(f)

    with ``I_leak(f) = Y_G @ u(f)`` from the reduced network
    ``(Y_L(f) + Y_G) u = i``. Uses the same one-shot assembly as
    :func:`solve_mutual_matrix` (reaction matrix and per-excitation leakage
    columns built once); only the field-point evaluation is generalised. The
    expensive 3D part is therefore frequency-independent and computed once, so
    a frequency sweep over ``I_leak(f)`` is cheap.

    ``world`` must already carry **all** ``nG`` grounding clusters (every
    anchor in ``anchors`` must be a real electrode in ``world``) and the soil
    model. Any current source is ignored — the excitations are applied
    internally, one cluster at a time, exactly as the historic per-excitation
    loop did.

    Parameters
    ----------
    world
        Assembled world holding all ``nG`` grounding clusters and the soil
        model. Sources, if any, are ignored.
    engine
        Engine configuration. ``segment_length`` controls the
        discretisation; ``frequencies`` selects the solved frequency via
        ``frequency_index``; ``image_max_terms`` / ``image_series_tol`` feed
        the 2-layer **self**-kernel (the reaction-matrix assembly), exactly
        as ``Engine.solve`` would.
    anchors
        Length-``nG`` list of electrode names, one per node. Anchor ``j``
        receives the unit current of excitation ``j`` and fills column ``j``
        of ``R``.
    field_points
        ``(M, 3)`` array of ``(x, y, depth)`` evaluation points, decoupled
        from ``nG``: ``M`` may differ from ``nG`` and the points are
        arbitrary (surface or buried). Row ``i`` becomes row ``i`` of ``R``.
    frequency_index
        Index into ``engine.frequencies``. Default 0. Only relevant when
        distributed conductors make the assembled system
        frequency-dependent.
    matrix_min_distance, matrix_max_terms, matrix_tol
        Truncation / clamp parameters of the field evaluation. The defaults
        ``1e-3 / 100 / 1e-6`` reproduce the hard-coded defaults of
        :meth:`FieldResult.potential` / ``_potential_two_layer`` exactly.
        (The reaction-matrix assembly uses ``engine.image_max_terms`` /
        ``image_series_tol`` independently, just like the live solver.)

    Returns
    -------
    R : np.ndarray, shape (M, nG), complex
        Galvanic Green-function matrix; ``R[i, j]`` is the potential at
        ``field_points[i]`` under unit excitation of cluster ``j``.

    Raises
    ------
    ValueError
        If ``world`` has no soil model, holds no electrodes, or
        ``field_points`` does not have shape ``(M, 3)``.
    KeyError
        If an entry of ``anchors`` is not an electrode in ``world``.
    NotImplementedError
        If the soil is a ``MultiLayerSoil`` with three or more layers
        (only homogeneous and 2-layer soils are supported).
    TypeError
        If the soil model type is not supported.

    See Also
    --------
    solve_mutual_matrix : Full ``nG × nG`` mutual grounding-impedance matrix.
    """
    if world.soil is None:
        raise ValueError("World has no soil model.")
    if not world.electrodes:
        raise ValueError("World contains no electrodes.")

    soil = world.soil
    nG = len(anchors)
    field_points = np.asarray(field_points, dtype=float)
    if field_points.ndim != 2 or field_points.shape[1] != 3:
        raise ValueError(f"field_points must have shape (M, 3), got {field_points.shape}.")

    # 1) Discretisation + topology — identical to solve_mutual_matrix.
    ds = engine.segment_length
    all_segments: list[_Segment] = []
    elec_to_segidx: dict[str, list[int]] = {}
    interfaces = (
        (world.soil.h_1,) if isinstance(world.soil, TwoLayerSoil) else None
    )
    for e in world.electrodes:
        segs = _discretize_electrode(e, ds, layer_interfaces=interfaces)
        elec_to_segidx[e.name] = list(range(len(all_segments), len(all_segments) + len(segs)))
        all_segments.extend(segs)

    cluster_id = _build_clusters(world.electrodes, world.conductors)
    finite_branches = _build_finite_branches(world.conductors, cluster_id)
    cond_segs, distributed_branches_objs, interior_nodes = (
        _build_distributed_topology(world.conductors, cluster_id)
    )
    pseudo_owners: list[str] = []
    for s in cond_segs:
        pn = s.electrode_name
        elec_to_segidx[pn] = [len(all_segments)]
        all_segments.append(s)
        cluster_id[pn] = pn
        pseudo_owners.append(pn)
    for n_ in interior_nodes:
        if n_ not in cluster_id:
            cluster_id[n_] = n_
            pseudo_owners.append(n_)
            elec_to_segidx[n_] = []
    n_lumped_branches = len(finite_branches)
    distributed_branch_tuples = [(db.node_a, db.node_b, db.R) for db in distributed_branches_objs]
    finite_branches = list(finite_branches) + distributed_branch_tuples

    earth_inductive_model = getattr(engine, "earth_inductive_model", "perfect_mirror")
    sigma_earth_for_carson: float | None = None
    layered_earth_for_sommerfeld: object = None
    if earth_inductive_model == "carson_series":
        from groundfield.coupling import resolve_earth_conductivity
        sigma_earth_for_carson = resolve_earth_conductivity(soil)
    elif earth_inductive_model == "sommerfeld":
        from groundfield.coupling import resolve_earth_layers
        layered_earth_for_sommerfeld = resolve_earth_layers(soil)
    inductance_matrix_full, has_inductance, carson_builder = _assemble_inductance_matrix(
        distributed_branches_objs,
        n_lumped_branches=n_lumped_branches,
        n_total_branches=len(finite_branches),
        earth_model=earth_inductive_model,
        sigma_earth=sigma_earth_for_carson,
        layered_earth=layered_earth_for_sommerfeld,
    )

    n_segments = len(all_segments)
    seg_points = np.array([s.midpoint for s in all_segments])
    seg_lengths = np.array([s.length for s in all_segments])
    wire_radii = np.array([s.wire_radius for s in all_segments])
    seg_shell_coeffs = np.array(
        [s.concrete_shell_coefficient_ohm_m for s in all_segments], dtype=float
    )

    # 2) Soil self-kernel — same case split as solve_mutual_matrix.
    rho_homog: float | None = None
    soil_two_layer: TwoLayerSoil | None = None
    if isinstance(soil, TwoLayerSoil):
        soil_two_layer = soil
    elif isinstance(soil, MultiLayerSoil):
        n_layers = len(soil.layers)
        if n_layers == 1:
            rho_homog = float(soil.layers[0].resistivity)
        elif n_layers == 2:
            soil_two_layer = TwoLayerSoil(
                rho_1=float(soil.layers[0].resistivity),
                rho_2=float(soil.layers[1].resistivity),
                h_1=float(soil.layers[0].thickness),
            )
        else:
            raise NotImplementedError(
                "solve_mutual_field supports homogeneous and 2-layer soils; "
                f"got MultiLayerSoil with n={n_layers} layers."
            )
    elif isinstance(soil, HomogeneousSoil):
        rho_homog = float(soil.resistivity)
    else:
        raise TypeError(f"Unsupported soil type {type(soil).__name__}.")

    if soil_two_layer is not None:
        self_max_terms = int(getattr(engine, "image_max_terms", 100))
        self_tol = float(getattr(engine, "image_series_tol", 1e-6))
        z_max = float(seg_points[:, 2].max())
        cross_layer = bool(z_max >= soil_two_layer.h_1)
        self_kernel = _two_layer_self_kernel_factory(
            soil_two_layer, self_max_terms, self_tol, allow_cross_layer=cross_layer
        )
    else:
        rho_for_kernel = float(rho_homog)

        def self_kernel(seg_pts, seg_lens, wr, currents):  # noqa: ANN001
            return _self_corrected_kernel(seg_pts, seg_lens, wr, currents, rho_for_kernel)

    eye = np.eye(n_segments)
    Z_seg_full = self_kernel(seg_points, seg_lengths, wire_radii, eye)
    omega = (
        2.0 * np.pi * float(engine.frequencies[frequency_index]) if has_inductance else 0.0
    )
    carson_dz = (
        carson_builder(omega) if (has_inductance and carson_builder is not None) else None
    )

    def _cached_bulk_self_kernel(seg_pts, seg_lens, wr, currents):  # noqa: ANN001
        return Z_seg_full @ currents

    # 3) Per-excitation leakage columns T (no diagonal/self term needed).
    T = np.zeros((n_segments, nG), dtype=complex)
    for j, anchor_j in enumerate(anchors):
        if anchor_j not in elec_to_segidx:
            raise KeyError(f"Anchor '{anchor_j}' is not an electrode in world.")
        elec_input_current = {e.name: 0j for e in world.electrodes}
        elec_input_current[anchor_j] = 1.0 + 0j
        elec_total = _solve_cluster_currents(
            electrodes=world.electrodes,
            elec_input_current=elec_input_current,
            cluster_id=cluster_id,
            seg_points=seg_points,
            seg_lengths=seg_lengths,
            wire_radii=wire_radii,
            elec_to_segidx=elec_to_segidx,
            self_kernel=_cached_bulk_self_kernel,
            finite_branches=finite_branches,
            pseudo_owners=pseudo_owners,
            omega=omega if has_inductance else 0.0,
            inductance_matrix=inductance_matrix_full if has_inductance else None,
            carson_correction=carson_dz,
            shell_coefficients=seg_shell_coeffs,
        )
        sc = np.zeros(n_segments, dtype=complex)
        for ename, idxs in elec_to_segidx.items():
            if not idxs:
                continue
            I_total = elec_total.get(ename, 0j)
            if I_total == 0j:
                continue
            L_total = seg_lengths[idxs].sum()
            sc[idxs] = I_total * seg_lengths[idxs] / L_total
        T[:, j] = sc

    # 4) Potential at the M field points per excitation (nG) — one kernel pass.
    if soil_two_layer is not None:
        R = _probe_potential_two_layer(
            field_points, seg_points, T, soil_two_layer,
            matrix_min_distance, max_terms=matrix_max_terms, tol=matrix_tol,
        )
    else:
        K_probe = _probe_kernel_homogeneous(
            field_points, seg_points, float(rho_homog), matrix_min_distance
        )
        R = K_probe @ T.real + 1j * (K_probe @ T.imag)
    return R

solve_mutual_matrix

solve_mutual_matrix(
    world: "World",
    engine: "Engine",
    anchors: list[str],
    probe_points: np.ndarray,
    *,
    frequency_index: int = 0,
    symmetrize: bool = True,
    matrix_min_distance: float = 0.001,
    matrix_max_terms: int = 100,
    matrix_tol: float = 1e-06
) -> np.ndarray

Full mutual grounding-impedance matrix Z_G in one assembly.

Reproduces the historic per-excitation construction (one Engine.solve per gnode; diagonal from :meth:FieldResult.cluster_impedance, off-diagonal from :meth:FieldResult.potential at the probe points) bit-for-bit, but assembles the segment reaction matrix and the probe-point field kernel only once.

world must already carry all nG grounding clusters (every anchor in anchors must be a real electrode in world). It must not carry any current source — the excitations are applied internally, one cluster at a time, exactly as the historic loop did. Finite-impedance / distributed conductors between distinct clusters are out of scope for the catalogue use case (each gnode is an isolated galvanic cluster) and are handled via the standard :func:_solve_cluster_currents machinery only in so far as they remain confined to a single excited cluster.

Parameters:

Name	Type	Description	Default
world	'World'	Assembled world holding all nG grounding clusters and the soil model. Sources, if any, are ignored.	required
engine	'Engine'	Engine configuration. segment_length controls the discretisation; frequencies selects the solved frequency via frequency_index; image_max_terms / image_series_tol feed the 2-layer self-kernel (the reaction-matrix assembly), exactly as Engine.solve would.	required
anchors	list[str]	Length-nG list of electrode names, one per gnode. Anchor j receives the unit current of excitation j and its cluster supplies the diagonal entry Z[j, j].	required
probe_points	ndarray	(nG, 3) array of off-diagonal probe coordinates: row i is the (x, y, depth) point at which the potential of every other excitation is sampled to fill column i of the off-diagonal. Mirrors res.potential([[x_i, y_i, depth]]).	required
frequency_index	int	Index into engine.frequencies. Default 0.	0
symmetrize	bool	If True (default) return 0.5 · (Z + Zᵀ) — the historic post-processing.	True
matrix_min_distance	float	Truncation / clamp parameters of the off-diagonal field evaluation. The defaults 1e-3 / 100 / 1e-6 reproduce the hard-coded defaults of :meth:FieldResult.potential / _potential_two_layer exactly. (The reaction-matrix assembly uses engine.image_max_terms / image_series_tol independently, just like the live solver.)	0.001
matrix_max_terms	float	Truncation / clamp parameters of the off-diagonal field evaluation. The defaults 1e-3 / 100 / 1e-6 reproduce the hard-coded defaults of :meth:FieldResult.potential / _potential_two_layer exactly. (The reaction-matrix assembly uses engine.image_max_terms / image_series_tol independently, just like the live solver.)	0.001
matrix_tol	float	Truncation / clamp parameters of the off-diagonal field evaluation. The defaults 1e-3 / 100 / 1e-6 reproduce the hard-coded defaults of :meth:FieldResult.potential / _potential_two_layer exactly. (The reaction-matrix assembly uses engine.image_max_terms / image_series_tol independently, just like the live solver.)	0.001

Returns:

Name	Type	Description
Z	(ndarray, shape(nG, nG), complex)	The mutual grounding-impedance matrix.

Source code in src/groundfield/solver/mutual.py

def solve_mutual_matrix(
    world: "World",
    engine: "Engine",
    anchors: list[str],
    probe_points: np.ndarray,
    *,
    frequency_index: int = 0,
    symmetrize: bool = True,
    matrix_min_distance: float = 1e-3,
    matrix_max_terms: int = 100,
    matrix_tol: float = 1e-6,
) -> np.ndarray:
    """Full mutual grounding-impedance matrix ``Z_G`` in one assembly.

    Reproduces the historic per-excitation construction (one
    ``Engine.solve`` per gnode; diagonal from
    :meth:`FieldResult.cluster_impedance`, off-diagonal from
    :meth:`FieldResult.potential` at the probe points) **bit-for-bit**,
    but assembles the segment reaction matrix and the probe-point field
    kernel only once.

    ``world`` must already carry **all** ``nG`` grounding clusters
    (every anchor in ``anchors`` must be a real electrode in
    ``world``). It must **not** carry any current source — the
    excitations are applied internally, one cluster at a time, exactly
    as the historic loop did. Finite-impedance / distributed conductors
    *between distinct clusters* are out of scope for the catalogue use
    case (each gnode is an isolated galvanic cluster) and are handled
    via the standard :func:`_solve_cluster_currents` machinery only in
    so far as they remain confined to a single excited cluster.

    Parameters
    ----------
    world
        Assembled world holding all ``nG`` grounding clusters and the
        soil model. Sources, if any, are ignored.
    engine
        Engine configuration. ``segment_length`` controls the
        discretisation; ``frequencies`` selects the solved frequency
        via ``frequency_index``; ``image_max_terms`` /
        ``image_series_tol`` feed the 2-layer **self**-kernel (the
        reaction-matrix assembly), exactly as ``Engine.solve`` would.
    anchors
        Length-``nG`` list of electrode names, one per gnode. Anchor
        ``j`` receives the unit current of excitation ``j`` and its
        cluster supplies the diagonal entry ``Z[j, j]``.
    probe_points
        ``(nG, 3)`` array of off-diagonal probe coordinates: row ``i``
        is the ``(x, y, depth)`` point at which the potential of every
        other excitation is sampled to fill column ``i`` of the
        off-diagonal. Mirrors ``res.potential([[x_i, y_i, depth]])``.
    frequency_index
        Index into ``engine.frequencies``. Default 0.
    symmetrize
        If ``True`` (default) return ``0.5 · (Z + Zᵀ)`` — the historic
        post-processing.
    matrix_min_distance, matrix_max_terms, matrix_tol
        Truncation / clamp parameters of the **off-diagonal** field
        evaluation. The defaults ``1e-3 / 100 / 1e-6`` reproduce the
        hard-coded defaults of :meth:`FieldResult.potential` /
        ``_potential_two_layer`` exactly. (The reaction-matrix
        assembly uses ``engine.image_max_terms`` / ``image_series_tol``
        independently, just like the live solver.)

    Returns
    -------
    Z : np.ndarray, shape (nG, nG), complex
        The mutual grounding-impedance matrix.
    """
    from groundfield.solver.engine import Engine  # noqa: F401  (typing only)

    if world.soil is None:
        raise ValueError("World has no soil model.")
    if not world.electrodes:
        raise ValueError("World contains no electrodes.")

    soil = world.soil
    nG = len(anchors)
    probe_points = np.asarray(probe_points, dtype=float)
    if probe_points.shape != (nG, 3):
        raise ValueError(
            f"probe_points must have shape ({nG}, 3), got "
            f"{probe_points.shape}."
        )

    # ------------------------------------------------------------------
    # 1) Discretisation + topology — identical to the live backends.
    # ------------------------------------------------------------------
    ds = engine.segment_length
    all_segments: list[_Segment] = []
    elec_to_segidx: dict[str, list[int]] = {}
    interfaces = (
        (world.soil.h_1,) if isinstance(world.soil, TwoLayerSoil) else None
    )
    for e in world.electrodes:
        segs = _discretize_electrode(e, ds, layer_interfaces=interfaces)
        elec_to_segidx[e.name] = list(
            range(len(all_segments), len(all_segments) + len(segs))
        )
        all_segments.extend(segs)

    cluster_id = _build_clusters(world.electrodes, world.conductors)
    finite_branches = _build_finite_branches(world.conductors, cluster_id)

    cond_segs, distributed_branches_objs, interior_nodes = (
        _build_distributed_topology(world.conductors, cluster_id)
    )
    pseudo_owners: list[str] = []
    for s in cond_segs:
        pn = s.electrode_name
        elec_to_segidx[pn] = [len(all_segments)]
        all_segments.append(s)
        cluster_id[pn] = pn
        pseudo_owners.append(pn)
    for n_ in interior_nodes:
        if n_ not in cluster_id:
            cluster_id[n_] = n_
            pseudo_owners.append(n_)
            elec_to_segidx[n_] = []
    n_lumped_branches = len(finite_branches)
    distributed_branch_tuples = [
        (db.node_a, db.node_b, db.R) for db in distributed_branches_objs
    ]
    finite_branches = list(finite_branches) + distributed_branch_tuples

    # Inductive coupling assembly (kept for parity; the catalogue probe
    # case has no distributed conductors, so ``has_inductance`` is False
    # and the system stays real/frequency-independent — matching the
    # historic ``image`` / ``image_2layer`` no-inductance fast path).
    earth_inductive_model = getattr(
        engine, "earth_inductive_model", "perfect_mirror"
    )
    sigma_earth_for_carson: float | None = None
    layered_earth_for_sommerfeld: object = None
    if earth_inductive_model == "carson_series":
        from groundfield.coupling import resolve_earth_conductivity

        sigma_earth_for_carson = resolve_earth_conductivity(soil)
    elif earth_inductive_model == "sommerfeld":
        from groundfield.coupling import resolve_earth_layers

        layered_earth_for_sommerfeld = resolve_earth_layers(soil)
    inductance_matrix_full, has_inductance, carson_builder = (
        _assemble_inductance_matrix(
            distributed_branches_objs,
            n_lumped_branches=n_lumped_branches,
            n_total_branches=len(finite_branches),
            earth_model=earth_inductive_model,
            sigma_earth=sigma_earth_for_carson,
            layered_earth=layered_earth_for_sommerfeld,
        )
    )

    n_segments = len(all_segments)
    seg_points = np.array([s.midpoint for s in all_segments])   # (N, 3)
    seg_lengths = np.array([s.length for s in all_segments])     # (N,)
    wire_radii = np.array([s.wire_radius for s in all_segments])  # (N,)
    seg_shell_coeffs = np.array(
        [s.concrete_shell_coefficient_ohm_m for s in all_segments],
        dtype=float,
    )

    # ------------------------------------------------------------------
    # 2) Soil-specific self-kernel closure — same dispatch as
    #    ``Engine.solve`` (image -> image_2layer for TwoLayerSoil).
    # ------------------------------------------------------------------
    rho_homog: float | None = None
    soil_two_layer: TwoLayerSoil | None = None
    if isinstance(soil, TwoLayerSoil):
        soil_two_layer = soil
    elif isinstance(soil, MultiLayerSoil):
        n_layers = len(soil.layers)
        if n_layers == 1:
            rho_homog = float(soil.layers[0].resistivity)
        elif n_layers == 2:
            soil_two_layer = TwoLayerSoil(
                rho_1=float(soil.layers[0].resistivity),
                rho_2=float(soil.layers[1].resistivity),
                h_1=float(soil.layers[0].thickness),
            )
        else:
            raise NotImplementedError(
                "solve_mutual_matrix supports homogeneous and 2-layer "
                f"soils; got MultiLayerSoil with n={n_layers} layers."
            )
    elif isinstance(soil, HomogeneousSoil):
        rho_homog = float(soil.resistivity)
    else:
        raise TypeError(f"Unsupported soil type {type(soil).__name__}.")

    if soil_two_layer is not None:
        # 2-layer reaction-matrix self-kernel: feed the engine's
        # image-series truncation, exactly like ``solve_image_2layer``.
        self_max_terms = int(getattr(engine, "image_max_terms", 100))
        self_tol = float(getattr(engine, "image_series_tol", 1e-6))
        z_max = float(seg_points[:, 2].max())
        cross_layer = bool(z_max >= soil_two_layer.h_1)
        self_kernel = _two_layer_self_kernel_factory(
            soil_two_layer, self_max_terms, self_tol,
            allow_cross_layer=cross_layer,
        )
    else:
        rho_for_kernel = float(rho_homog)

        def self_kernel(seg_pts, seg_lens, wr, currents):  # noqa: ANN001
            return _self_corrected_kernel(
                seg_pts, seg_lens, wr, currents, rho_for_kernel
            )

    # ------------------------------------------------------------------
    # 3) Reaction matrix ONCE.  ``Z_seg_full = self_kernel(eye)`` is the
    #    N×N segment reaction matrix used by ``mom._build_Z_*``; column
    #    k is ``self_kernel(e_k)``. Reusing it for every excitation makes
    #    the diagonal ``ph`` evaluation and the cluster-current assembly
    #    a sequence of matvecs instead of nG full kernel calls.
    # ------------------------------------------------------------------
    eye = np.eye(n_segments)
    Z_seg_full = self_kernel(seg_points, seg_lengths, wire_radii, eye)  # (N, N)

    # ADR-0012 V2 shell augmentation (added to the *diagonal* of the
    # reaction matrix exactly as in ``_solve_cluster_currents`` and the
    # post-solve ``ph`` step).
    shell_diag: np.ndarray | None = None
    if np.any(seg_shell_coeffs > 0.0):
        with np.errstate(divide="ignore", invalid="ignore"):
            shell_diag = np.where(
                seg_lengths > 0.0,
                seg_shell_coeffs / seg_lengths,
                0.0,
            )

    # A linear ``self_kernel``-equivalent backed by the precomputed
    # matrix: ``Z_seg_full @ currents`` reproduces ``self_kernel(...)``
    # for the unit-current assembly inside ``_solve_cluster_currents``.
    # The shell augmentation is applied on top (the assembly inside
    # ``_solve_cluster_currents`` does the same: bulk kernel then
    # ``phi += shell_diag · currents``), so we deliberately pass a
    # *bulk-only* kernel here and let ``shell_coefficients`` re-add the
    # shell — keeping the exact same operation order as the live solver.
    def _cached_bulk_self_kernel(seg_pts, seg_lens, wr, currents):  # noqa: ANN001
        return Z_seg_full @ currents

    omega = (
        2.0 * np.pi * float(engine.frequencies[frequency_index])
        if has_inductance
        else 0.0
    )
    carson_dz = (
        carson_builder(omega)
        if (has_inductance and carson_builder is not None)
        else None
    )

    # ------------------------------------------------------------------
    # 4) Per-excitation leakage distribution.  Each gnode is an isolated
    #    galvanic cluster, so injecting 1 A into cluster ``j`` activates
    #    only that cluster; the resulting per-segment leakage column
    #    ``sc_j`` is non-zero on cluster-``j`` segments only — exactly
    #    what the historic single-source solve produced.
    # ------------------------------------------------------------------
    real_electrode_names = {e.name for e in world.electrodes}
    # cluster members per real electrode (== FieldResult.clusters), used
    # for the diagonal ``cluster_impedance`` average and current sum.
    cluster_members: dict[str, list[str]] = {}
    for ename in real_electrode_names:
        cluster_members[ename] = sorted(
            n
            for n in cluster_id
            if cluster_id[n] == cluster_id[ename] and n in real_electrode_names
        )

    T = np.zeros((n_segments, nG), dtype=complex)  # column j = sc_j
    Z = np.zeros((nG, nG), dtype=complex)

    for j, anchor_j in enumerate(anchors):
        if anchor_j not in elec_to_segidx:
            raise KeyError(f"Anchor '{anchor_j}' is not an electrode in world.")
        elec_input_current = {e.name: 0j for e in world.electrodes}
        elec_input_current[anchor_j] = 1.0 + 0j

        elec_total = _solve_cluster_currents(
            electrodes=world.electrodes,
            elec_input_current=elec_input_current,
            cluster_id=cluster_id,
            seg_points=seg_points,
            seg_lengths=seg_lengths,
            wire_radii=wire_radii,
            elec_to_segidx=elec_to_segidx,
            self_kernel=_cached_bulk_self_kernel,
            finite_branches=finite_branches,
            pseudo_owners=pseudo_owners,
            omega=omega if has_inductance else 0.0,
            inductance_matrix=inductance_matrix_full if has_inductance else None,
            carson_correction=carson_dz,
            shell_coefficients=seg_shell_coeffs,
        )

        # Reconstruct the per-segment leakage column (image.py:1455-1463).
        sc = np.zeros(n_segments, dtype=complex)
        for ename, idxs in elec_to_segidx.items():
            if not idxs:
                continue
            I_total = elec_total.get(ename, 0j)
            if I_total == 0j:
                continue
            L_total = seg_lengths[idxs].sum()
            sc[idxs] = I_total * seg_lengths[idxs] / L_total
        T[:, j] = sc

        # -- Diagonal Z[j, j] == cluster_impedance(anchor_j) -----------
        # Build ``ph`` (potential at every segment midpoint) exactly as
        # the post-solve step of the live solver: bulk self-kernel on
        # the real and imaginary leakage, then the shell augmentation.
        members = cluster_members.get(anchor_j, [anchor_j])
        i_sum = sum(elec_total.get(m, 0j) for m in members)
        member0_idxs = elec_to_segidx[members[0]]
        if sc.any():
            phi_re = Z_seg_full @ sc.real
            phi_im = Z_seg_full @ sc.imag
            if shell_diag is not None:
                phi_re = phi_re + shell_diag * sc.real
                phi_im = phi_im + shell_diag * sc.imag
            ph_member0 = (phi_re + 1j * phi_im)[member0_idxs]
            u_cluster = complex(np.mean(ph_member0))
        else:
            u_cluster = 0j
        Z[j, j] = (u_cluster / i_sum) if i_sum != 0 else complex("nan")

    # ------------------------------------------------------------------
    # 5) Off-diagonal block.  ``Z[i, j] = potential at probe_i due to the
    #    leakage column sc_j`` — the probe kernel is built ONCE and
    #    applied to all excitation columns at once.
    # ------------------------------------------------------------------
    if soil_two_layer is not None:
        Phi = _probe_potential_two_layer(
            probe_points,
            seg_points,
            T,
            soil_two_layer,
            matrix_min_distance,
            max_terms=matrix_max_terms,
            tol=matrix_tol,
        )  # (nG, nG): row i = probe i, col j = excitation j
    else:
        K_probe = _probe_kernel_homogeneous(
            probe_points, seg_points, float(rho_homog), matrix_min_distance
        )  # (nG, N)
        Phi = K_probe @ T.real + 1j * (K_probe @ T.imag)

    for j in range(nG):
        for i in range(nG):
            if i == j:
                continue
            Z[i, j] = Phi[i, j]

    if symmetrize:
        return 0.5 * (Z + Z.T)
    return Z