Skip to content

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.

The series controls Engine.image_max_terms (default 300) and Engine.image_series_tol (default \(10^{-6}\)) are forwarded on both dispatch paths. Up to 0.14.x image_nlayer was called without them and fell back to its own signature defaults (200, \(10^{-6}\)), so raising image_max_terms — the documented remedy for a high layer contrast, \(|K| \to 1\), where the geometric tail bound \(|K|^{n+1} / (1 - |K|)\) decays slowly — had no effect when the same soil was spelled as a MultiLayerSoil instead of a TwoLayerSoil. Since 0.15.0 both spellings return the identical impedance.

mom_sommerfeld accuracy controls

Since 0.15.0 Engine.solve also forwards four fields to mom_sommerfeld. Up to 0.14.x that backend was called as solve_mom_sommerfeld(world, self) with no keywords, so lambda_max_factor, epsabs and epsrel were unreachable from the public API — a user who hit a non-converged corner in the package's reference engine had nothing to tighten.

Engine field backend keyword default meaning
sommerfeld_lambda_max_factor lambda_max_factor 200.0 truncation of the reflected-remainder quadrature in units of \(1/d_\text{rem}\); saturates at 60
sommerfeld_epsabs epsabs \(10^{-9}\) absolute threshold of the grid-refinement convergence test
sommerfeld_epsrel epsrel \(10^{-7}\) relative threshold, measured against \(\lvert G\rvert\)
sommerfeld_max_osc_panels max_osc_panels 200 000 budget for the \(\pi/s\) panel family resolving \(J_0(\lambda s)\)

The defaults equal the previous signature defaults, so existing callers get bit-identical results. All four are echoed into FieldResult.metadata. They are not exposed through gf.create_engine(...) — construct the model directly:

import groundfield as gf

engine = gf.Engine(
    backend="mom_sommerfeld",
    segment_length=0.5,
    sommerfeld_epsrel=1e-9,
    sommerfeld_max_osc_panels=2_000_000,
)

lambda_max_factor has two different units

On the top-layer path it is in units of \(1/d_\text{rem}\) (the remainder decay length) and saturates at 60, so raising it is harmless. On the delegated \(n = 2\) cross-layer path (z > h_1 or z_s > h_1, forwarded to coupling.layered_green.two_layer_real_space_kernel) it is in units of \(1/\ell_\text{char}\) (a geometric length), is not saturated, and raising it exhausts that kernel's Hankel node budget and triggers SommerfeldResolutionWarning — there the remedy is to lower it. See mom_sommerfeld for the full note.

SommerfeldConvergenceWarning

groundfield.solver.mom_sommerfeld.SommerfeldConvergenceWarning is a UserWarning subclass raised when the reflected-remainder quadrature exhausts its refinement ladder without two successive grids agreeing within max(epsabs, epsrel · |G|). It reports the achieved versus the requested tolerance, the offending \((s, \lvert z - z_s\rvert, z + z_s)\), whether the oscillation panel budget bound and how much of the required \(\lambda\)-range was covered, that the returned value's sign may be wrong, and which parameter to raise. An \(O(N^2)\) matrix assembly aggregates all failures into one warning naming the worst pair and the number of affected entries.

Up to 0.15.0 the same condition was only logger.debug, and the least accurate rung of the refinement ladder was the value returned; that combination could return a negative Green's function (−9.47e-05 where the far-field asymptote requires +5.00e-04) with zero Python-level warnings. Two mechanisms reach the warning:

  • oscillation budget — sommerfeld_max_osc_panels binds for \(s/d_\text{rem} \gtrsim 1.9\cdot10^{4}\) (a centimetre-thin intermediate layer plus a kilometre-scale separation). Remedy: raise it; cost grows linearly.
  • layer contrast — \(\lvert K_1\rvert \to 1\). Remedy: none of the knobs; use image_2layer, whose geometric series is exact for \(n = 2\).

Silence with warnings.simplefilter("ignore", SommerfeldConvergenceWarning), or promote it with simplefilter("error", ...) in a validation harness.

Sibling categories in the same family: solver.image_2layer.SeriesTruncationWarning (image series hit image_max_terms) and coupling.layered_green.SommerfeldResolutionWarning (Hankel grid cannot resolve \(J_0\)).

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.

Diagnostics raised by the galvanic image family

ShallowSegmentWarning

Exported from groundfield.solver (and from groundfield.solver.image, where it is defined):

import warnings
from groundfield.solver import ShallowSegmentWarning

warnings.simplefilter("ignore", ShallowSegmentWarning)

A UserWarning subclass emitted when a leakage segment is longer than four times its burial depth, \(L_i > 4 z_i\). The reaction-matrix diagonal then uses the point image at distance \(2 z_i\) where the segment's image line would be needed, which overestimates that entry by more than 13 % (unboundedly as \(z \to 0\)), so the grounding impedance is biased high and keeps moving under mesh refinement. The message quantifies the bias of the worst offender with the horizontal-segment closed form.

What triggers it, honestly:

  • Not only distributed conductors. Any geometry whose segments are long compared with the burial depth qualifies — a ring of radius 25 m at \(z = 0.8\) m discretised at segment_length = 5 m gives arc segments of 4.909 m and warns "by 26 %"; so do coarse meshes, counterpoises and surface-near tapes. Measured across the test suite it fires at 20 test nodes in 9 modules (test_carson_coupling, test_conductors_distributed, test_conductors_finite, test_convergence, test_diagnostics, test_earth_return_benchmark, test_inductive_coupling, test_sommerfeld_inductance, test_pass9_mom_sommerfeld).
  • Only from image and image_2layer. mom, bem, cim, mom_sommerfeld and mutual share the same biased diagonal through solver.image._self_corrected_kernel, but they do not emit the warning, so a compare_engines run warns asymmetrically about a bias all of them carry. The reasoning is recorded in _check_segment_depths ("Emission scope"): the shared kernel is re-entered per frequency and per excitation, so emitting there would repeat the message with a stacklevel pointing into solver internals.

The hard depth guards are not asymmetric: a segment midpoint at \(z \le 0\) or with \(2|z| \le a\) raises ValueError from the shared kernel, so all backends of the family reject it identically (0.15.0; before that only image / image_2layer rejected \(z \le 0\), and mom / bem / cim / mom_sommerfeld returned a number for an airborne electrode). A segment that straddles \(z = 0\) — midpoint buried, upper end in the air — needs the segment direction and is therefore rejected in image / image_2layer only.

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:

The reported deviation is the worst pairwise relative difference of the cluster impedances, (max Z - min Z) / min |Z| — see Validation for the definition and for the change of metric in 0.15.0.

  • 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). ShallowSegmentWarning Diagnostic warning category of the galvanic image family: a leakage segment is longer than four times its burial depth, so the point-image stand-in for the segment's image line biases the reaction-matrix diagonal high. Re-exported here so callers can write warnings.simplefilter("ignore", ShallowSegmentWarning) without importing from the backend module (groundfield.solver.image).

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.
image_max_terms, image_series_tol

Controls of the Tagg / Sunde image-charge series, forwarded to every layered backend.
sommerfeld_lambda_max_factor, sommerfeld_epsabs, sommerfeld_epsrel, sommerfeld_max_osc_panels

Accuracy controls of the mom_sommerfeld reflected-remainder quadrature, forwarded to :func:groundfield.solver.mom_sommerfeld.solve_mom_sommerfeld (new in 0.15.0 — up to 0.14.x that backend was called with no keywords at all, so none of these was reachable through the public API). The defaults reproduce the previous behaviour exactly. sommerfeld_max_osc_panels is the escape hatch for the corner in which the backend raises :class:~groundfield.solver.mom_sommerfeld.SommerfeldConvergenceWarning. Not exposed through :func:groundfield.create_engine; construct the engine directly, e.g. gf.Engine(backend="mom_sommerfeld", sommerfeld_max_osc_panels=2_000_000).

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 
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
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

        # ``image_nlayer`` re-dispatches to ``solve_image_2layer``
        # for n = 2, so the series controls must be forwarded here
        # as well — otherwise the backend silently falls back to
        # its own signature defaults (200, 1e-6) and the engine
        # configuration promised by ``Engine.image_max_terms`` is
        # discarded on the multilayer path.
        return solve_image_nlayer(
            world, self,
            max_terms=self.image_max_terms, tol=self.image_series_tol,
        )

    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

        # Up to 0.14.x this call passed no keywords, so the backend
        # fell back to its own signature defaults and the accuracy
        # knobs of the *reference* engine were unreachable from the
        # public API — a user who hit the non-converged corner had
        # nothing to tighten. Forwarded following the
        # ``image_max_terms`` / ``image_series_tol`` precedent above;
        # the field defaults equal the old signature defaults, so
        # existing callers are unaffected.
        return solve_mom_sommerfeld(
            world, self,
            lambda_max_factor=self.sommerfeld_lambda_max_factor,
            epsabs=self.sommerfeld_epsabs,
            epsrel=self.sommerfeld_epsrel,
            max_osc_panels=self.sommerfeld_max_osc_panels,
        )

    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 
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
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]
    # The scalar contract holds for *both* paths. Checking it here
    # as well as in the validator keeps the error attributable to
    # ``with_frequencies()`` rather than to a nested model rebuild.
    _check_frequency_values(freqs)
    data = self.model_dump()
    data["frequencies"] = freqs
    if preserve_order:
        # Re-validate the rest, but suppress the field validator's
        # monotone warning — a deliberately ordered sweep is the
        # documented opt-in.
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", EngineFrequencyOrderWarning)
            return self.__class__.model_validate(data)
    # Default path: run the full field validation, including the
    # monotone-order warning. ``model_copy(update=...)`` would skip
    # every validator (pydantic v2 semantics), which is the
    # pre-0.14.1 defect. Re-emit any warning at the caller's frame
    # so the notebook line that built the engine is the one shown.
    with warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter("always")
        engine = self.__class__.model_validate(data)
    for entry in caught:
        warnings.warn(entry.message, entry.category, stacklevel=2)
    return engine

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, or if none of its members carries potential and current data.
Notes

The cluster potential is the length-weighted (Galerkin) average over all members,

.. math:: \varphi_c \;=\; \frac{\sum_{e \in c} L_e\, \varphi_e} {\sum_{e \in c} L_e}, \qquad L_e = \sum_{k \in e} L_k ,

with the per-electrode lengths taken from :attr:point_sources. Since every member potential :math:\varphi_e is itself the length-weighted average over that electrode's segments, :math:\varphi_c is exactly the length-weighted average over the whole cluster surface — the quantity the augmented system constrains to one value (see solver/image.py, block 1 of the multi-port system).

Up to 0.14.1 this method read members[0], i.e. the alphabetically first member name. Because the reported member potentials were unweighted segment means they did not agree with each other, so renaming an electrode changed the reported cluster impedance of an unchanged geometry (19.4393 Ω → 19.2314 Ω in the audit case, 2026-07-28, F11). With the weighted reduction the members agree to round-off for the image family and the result is invariant under renaming for every backend.

Numerator and denominator run over different member sets, deliberately:

  • :math:\varphi_c averages over members that carry both potential and current data — a member with no potential cannot contribute a term to the average.
  • :math:\sum_e I_e sums over every member that carries current data, whether or not its potential was stored, because its injected current is part of the cluster total regardless.

Every backend in this package populates both mappings for every member, so the two sets coincide in practice. They differ only for hand-built :class:FieldResult instances, where restricting the current sum to the potential-carrying members would silently under-count the cluster current (audit 2026-07-29).

Source code in src/groundfield/solver/result.py 
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
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, or if none of its members
        carries potential *and* current data.

    Notes
    -----
    The cluster potential is the **length-weighted (Galerkin)
    average** over all members,

    .. math::
        \\varphi_c \\;=\\;
        \\frac{\\sum_{e \\in c} L_e\\, \\varphi_e}
              {\\sum_{e \\in c} L_e},
        \\qquad
        L_e = \\sum_{k \\in e} L_k ,

    with the per-electrode lengths taken from
    :attr:`point_sources`. Since every member potential
    :math:`\\varphi_e` is itself the length-weighted average over
    that electrode's segments, :math:`\\varphi_c` is exactly the
    length-weighted average over the whole cluster surface — the
    quantity the augmented system constrains to one value (see
    ``solver/image.py``, block 1 of the multi-port system).

    Up to 0.14.1 this method read ``members[0]``, i.e. the
    *alphabetically first* member name. Because the reported
    member potentials were unweighted segment means they did not
    agree with each other, so renaming an electrode changed the
    reported cluster impedance of an unchanged geometry (19.4393 Ω
    → 19.2314 Ω in the audit case, 2026-07-28, F11). With the
    weighted reduction the members agree to round-off for the
    image family and the result is invariant under renaming for
    every backend.

    Numerator and denominator run over **different** member sets,
    deliberately:

    * :math:`\\varphi_c` averages over members that carry *both*
      potential and current data — a member with no potential
      cannot contribute a term to the average.
    * :math:`\\sum_e I_e` sums over every member that carries
      current data, whether or not its potential was stored,
      because its injected current is part of the cluster total
      regardless.

    Every backend in this package populates both mappings for
    every member, so the two sets coincide in practice. They
    differ only for hand-built :class:`FieldResult` instances,
    where restricting the current sum to the potential-carrying
    members would silently under-count the cluster current
    (audit 2026-07-29).
    """
    members = self.clusters.get(electrode_name, [electrode_name])
    if not members:
        raise KeyError(f"No cluster for '{electrode_name}'.")
    usable = [
        m for m in members
        if m in self.electrode_potentials and m in self.electrode_currents
    ]
    if not usable:
        raise KeyError(
            f"No potential/current data for any member of the cluster "
            f"of '{electrode_name}' (members: {members})."
        )

    # Per-electrode total leakage length (weights of the Galerkin
    # average). Falls back to equal weights when the backend did
    # not populate ``point_sources`` (stub results).
    lengths: dict[str, float] = {}
    for ps in self.point_sources:
        lengths[ps.electrode_name] = (
            lengths.get(ps.electrode_name, 0.0) + float(ps.length)
        )
    w = np.array([lengths.get(m, 0.0) for m in usable], dtype=float)
    if w.sum() <= 0.0:
        w = np.ones(len(usable), dtype=float)

    n_freq = len(self.electrode_potentials[usable[0]])
    u = [
        complex(
            sum(
                float(wm) * self.electrode_potentials[m][k]
                for wm, m in zip(w, usable)
            )
            / float(w.sum())
        )
        for k in range(n_freq)
    ]
    # The *denominator* runs over every member that has current
    # data, not over ``usable``: a member without potential data
    # cannot contribute to the weighted average, but its injected
    # current is still part of the cluster's total. Restricting the
    # sum to ``usable`` (0.15.0 development state) silently dropped
    # it and returned a too-large impedance instead of the
    # ``KeyError`` v0.14.1 raised — e.g. potentials {"a"} with
    # currents {"a": 1 A, "b": 1 A} and cluster ["a", "b"] gave
    # 10 V / 1 A = 10 Ω where 10 V / 2 A = 5 Ω is meant
    # (audit 2026-07-29). Not reachable through any shipped
    # backend, but ``FieldResult`` is public and user-constructible.
    current_members = [
        m for m in members
        if m in self.electrode_currents
        and len(self.electrode_currents[m]) >= n_freq
    ]
    i_sum = [
        sum(self.electrode_currents[m][k] for m in current_members)
        for k in range(n_freq)
    ]
    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 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
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 
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
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 
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
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.

ShallowSegmentWarning

Bases: UserWarning

A leakage segment is shallow compared with its own length.

The diagonal self-image entry of the reaction matrix uses the point image at distance \(2z\) instead of the image line of the segment, i.e. the \(L \ll 4 z\) limit

.. math:: \frac{1}{2 z} \;=\; \lim_{L \to 0} \frac{2}{L}\,\operatorname{arsinh}!\frac{L}{4 z}.

For \(L \gtrsim 4 z\) the point form overestimates that entry by more than 13 % (and unboundedly as \(z \to 0\)), so the reported grounding impedance is biased high and keeps moving under mesh refinement. Typical triggers: a counterpoise or PEN conductor discretised with segments much longer than its burial depth, surface-near tapes, and rings or meshes whose arc/wire segments are longer than four times the burial depth (a ring of radius 25 m at \(z = 0.8\) m with segment_length = 5 m gives 4.909 m arc segments and warns "by 26 %") — the warning is not restricted to distributed conductors.

Scope of the quantified bias

The percentage quoted in the message is the horizontal-segment figure \(\bigl(L/(4z)\bigr)/\operatorname{arsinh} \bigl(L/(4z)\bigr) - 1\). That is the applicable figure for every segment that can still reach the warning: a segment inclined by more than \(30°\) to the horizontal and satisfying \(L > 4 z\) would have its upper end above \(z = 0\), which :func:_check_segment_depths rejects with a ValueError first (proof: \(L > 4 z\) together with \(L\,|e_z| \le 2 z\) forces \(|e_z| < 1/2\)). The message states the inclination it assumed so the reader can check that reasoning.

Emission scope

The warning is raised by :func:solve_image and :func:~groundfield.solver.image_2layer.solve_image_2layer only, even though mom, bem, cim, mom_sommerfeld and mutual share the biased diagonal through :func:_self_corrected_kernel. See the "Emission scope" note in :func:_check_segment_depths for why the hard guard is shared but the advisory is not.

Silence with warnings.simplefilter("ignore", ShallowSegmentWarning) once the bias has been accepted, or refine Engine.segment_length / Conductor.discretize_segment_length until the impedance stabilises.

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 
2182
2183
2184
2185
2186
2187
2188
2189
2190
2191
2192
2193
2194
2195
2196
2197
2198
2199
2200
2201
2202
2203
2204
2205
2206
2207
2208
2209
2210
2211
2212
2213
2214
2215
2216
2217
2218
2219
2220
2221
2222
2223
2224
2225
2226
2227
2228
2229
2230
2231
2232
2233
2234
2235
2236
2237
2238
2239
2240
2241
2242
2243
2244
2245
2246
2247
2248
2249
2250
2251
2252
2253
2254
2255
2256
2257
2258
2259
2260
2261
2262
2263
2264
2265
2266
2267
2268
2269
2270
2271
2272
2273
2274
2275
2276
2277
2278
2279
2280
2281
2282
2283
2284
2285
2286
2287
2288
2289
2290
2291
2292
2293
2294
2295
2296
2297
2298
2299
2300
2301
2302
2303
2304
2305
2306
2307
2308
2309
2310
2311
2312
2313
2314
2315
2316
2317
2318
2319
2320
2321
2322
2323
2324
2325
2326
2327
2328
2329
2330
2331
2332
2333
2334
2335
2336
2337
2338
2339
2340
2341
2342
2343
2344
2345
2346
2347
2348
2349
2350
2351
2352
2353
2354
2355
2356
2357
2358
2359
2360
2361
2362
2363
2364
2365
2366
2367
2368
2369
2370
2371
2372
2373
2374
2375
2376
2377
2378
2379
2380
2381
2382
2383
2384
2385
2386
2387
2388
2389
2390
2391
2392
2393
2394
2395
2396
2397
2398
2399
2400
2401
2402
2403
2404
2405
2406
2407
2408
2409
2410
2411
2412
2413
2414
2415
2416
2417
2418
2419
2420
2421
2422
2423
2424
2425
2426
2427
2428
2429
2430
2431
2432
2433
2434
2435
2436
2437
2438
2439
2440
2441
2442
2443
2444
2445
2446
2447
2448
2449
2450
2451
2452
2453
2454
2455
2456
2457
2458
2459
2460
2461
2462
2463
2464
2465
2466
2467
2468
2469
2470
2471
2472
2473
2474
2475
2476
2477
2478
2479
2480
2481
2482
2483
2484
2485
2486
2487
2488
2489
2490
2491
2492
2493
2494
2495
2496
2497
2498
2499
2500
2501
2502
2503
2504
2505
2506
2507
2508
2509
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,)

    # 3c) Burial-depth validation (audit 2026-07-28, F26/F29): the
    #     galvanic image construction needs a finite image separation
    #     2 z >> a. Reject surface-laid / airborne conductors and warn
    #     for the shallow regime where the point-image stand-in biases
    #     the diagonal.
    _check_segment_depths(all_segments, where="solve_image")

    # 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 = _weighted_node_potential(
            phi_per_freq, idxs, seg_lengths, 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 
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
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])

    # Burial-depth validation (audit 2026-07-28, F26/F29). The air
    # mirror of the Tagg/Sunde series is the same perfect mirror as in
    # the homogeneous backend, so the same finite-image-separation
    # requirement 2 z >> a applies to every leakage segment.
    _check_segment_depths(all_segments, where="solve_image_2layer")
    # 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
        # Length-weighted (Galerkin) node potential — the same
        # reduction the multi-port matrix enforces inside
        # _solve_cluster_currents (audit 2026-07-28, F10/F11). The
        # rod-split-at-the-interface case of ADR-0007 is exactly where
        # the historic unweighted mean diverged from the solved value.
        u_list = _weighted_node_potential(
            phi_per_freq, idxs, seg_lengths, 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 
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
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 
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
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 
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
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