image — homogeneous image-charge sum

Physical context

A grounding electrode embedded in a homogeneous half-space of resistivity \(\rho\) injects a current \(I\) into the soil. The quasi-static potential field is governed by Laplace's equation (\(-\nabla \cdot (\sigma \nabla \varphi) = q\) with \(\sigma = 1/\rho\) and \(q\) the source-current density), with an insulating boundary at the soil surface (\(\partial\varphi/\partial z = 0\) at \(z = 0\)) and \(\varphi \to 0\) at infinity.

For frequencies \(f < 1\,\text{kHz}\) the displacement-current term is negligible — the relaxation time of moist soil (\(\tau = \varepsilon/\sigma\)) is on the order of 100 ns, well below the millisecond regime. The static potential field is therefore representative of the entire quasi-static frequency window.

Governing equation: image-charge solution

A point current source \(I\) at depth \(z_s > 0\) in a homogeneous half-space satisfies the Neumann boundary at \(z = 0\) exactly through the image-charge construction: place a virtual source of identical strength at the mirror position \(z = -z_s\). The superposition of source and image gives

\[ \varphi(x, y, z) \;=\; \frac{\rho\, I}{4\pi} \left( \frac{1}{r} + \frac{1}{r'} \right), \qquad r = \sqrt{(x{-}x_s)^2 + (y{-}y_s)^2 + (z{-}z_s)^2}, \]

with \(r' = \sqrt{(x{-}x_s)^2 + (y{-}y_s)^2 + (z{+}z_s)^2}\) the distance to the air-mirrored image.

This is the smallest, cleanest closed form in the engine family and provides the baseline against which every layered engine collapses when its layer contrast vanishes.

Numerical strategy

Wire-segment discretisation

A finite electrode (rod, ring, mesh) is discretised into \(N\) collinear segments of length \(L_i \le \Delta s\) (the engine.segment_length parameter). Each segment carries one point current source at its midpoint. The total current of an electrode is distributed uniformly per unit length across its segments — i.e. the segment current \(I_i = I_{\text{electrode}} \cdot L_i / \sum_j L_j\).

The uniform-current ansatz is an approximation: the true current distribution along a wire is non-uniform, with end-point concentrations on the order of \(\sim 5\,\%\). This residual is handled either by accepting a \(\sim 5\,\%\) Dwight-bias on the input impedance (cheap), or by switching to the mom backend which solves for the actual distribution at \(O(N^3)\) cost.

Self-action correction

For the average-potential evaluation at segment midpoints, the diagonal of the kernel matrix carries a \(1/r\) singularity that the point-source representation cannot handle. We replace the direct-source self-distance by the analytical line self-potential

\[ \varphi_{\text{self,line}} \;=\; \frac{\rho\, I_i}{2\pi\, L_i}\, \ln\!\frac{L_i}{a_i}, \]

with \(a_i\) the wire radius. This is the classic Howe / Sunde average-potential formula for a thin wire of finite length.

The image contribution at the same segment is evaluated as a point image at distance \(2 z_i\),

\[ K_{ii}^{\text{img}} \;=\; \frac{1}{2 z_i} \;=\; \lim_{L_i \to 0}\; \frac{2}{L_i}\, \operatorname{arsinh}\!\frac{L_i}{4 z_i}, \]

i.e. the \(L_i \ll 4 z_i\) limit of the exact potential of the segment's image line at the segment midpoint. Two consequences are enforced explicitly since 0.15.0 (they were silent before):

• Hard error for vanishing image separation. As \(z_i \to 0\) the image line merges with the conductor and the self-image term diverges — like \(1/z\) in the point form, like \(\ln\bigl(L/(2z)\bigr)\) in the exact form. A conductor lying in the surface plane has no image separation at all, and its correct model is a coincident source and image (twice the free-space line self-potential), which is a different kernel from the one assembled here. Any segment whose midpoint has \(z \le 0\) or \(2|z| \le a\) therefore raises ValueError — the same criterion the inductive path uses in build_inductance_matrix. Up to 0.14.1 the term was instead clamped at _MIN_DISTANCE = 1 mm, which returned a mesh-dependent grounding resistance up to ~15x too high (a 10 m tape at \(z = 0\): 810 Ω at \(\Delta s = 1\) m, 413 Ω at 0.5 m, 98 Ω at 0.1 m — never converging, against a physical ~23 Ω), and an electrode at negative \(z\) (in air) was silently mirrored into the soil.

Both midpoint conditions are checked inside the shared reaction-matrix kernel _self_corrected_kernel, so mom, bem, cim, mom_sommerfeld and mutual reject them identically. In 0.15.0 development the \(z \le 0\) half still sat in solve_image / solve_image_2layer alone, so compare_engines on an airborne ring (\(z = -100\) m) raised for the image family and returned 5.5162 Ω for the other four.

Who is named in the message. solve_image / solve_image_2layer pass owner labels down, so the message reads electrode 'g1' — or conductor 'pen' for a galvanically coupled distributed conductor, whose leakage segments carry the internal pseudo-node name __cond_pen__seg_0 that the user never wrote. The sibling backends call the kernel positionally and get the segment index plus its (x, y, z) coordinates instead; the subject of the sentence degrades to "the model" there. - Hard error for a segment that straddles the surface. A segment whose midpoint is buried but whose upper end \(z_{\text{top}} = z_i - \tfrac{1}{2} L_i |e_{z,i}|\) is negative has part of its conductor in the air. Checking midpoints only let RodElectrode(position=(0, 0, -0.5), length=1.2) solve to 106.2444 Ω with nothing but a ShallowSegmentWarning. This check needs the segment tangent (_Segment.direction), which the shared kernel does not receive, so it runs in solve_image / solve_image_2layer only. A segment whose upper end is exactly at \(z = 0\) stays legal: that is the ordinary driven rod (position=(x, y, 0)), whose source-plus-image is the length-\(2L\) line Dwight's closed form is derived from. - ShallowSegmentWarning for \(L_i > 4 z_i\). The regime is legal but the point image then overestimates \(K_{ii}^{\text{img}}\) by more than 13 % (unboundedly as \(z \to 0\)), so the impedance is biased high and still moves under refinement. The warning quantifies the bias of the worst offender with the horizontal-segment formula \(\bigl(L/4z\bigr)/\operatorname{arsinh}\bigl(L/4z\bigr) - 1\); measured at \(z = 0.01\) m, \(\Delta s = 0.5\) m: 35.5 Ω against Dwight's 21.0 Ω (+69 %).

That formula is the applicable one for every segment that can still reach the warning, because the straddling guard above removes the steep ones: \(L > 4 z\) together with \(L |e_z| \le 2 z\) forces \(|e_z| < 1/2\), i.e. an inclination below \(30°\) to the horizontal. The message states the inclination it assumed. Up to 0.15.0 development the same horizontal figure was quoted for vertical worst offenders, where the vertical closed form \((1/L)\ln\bigl((4z+L)/(4z-L)\bigr)\) has no real value at all (\(4z < L\) makes the argument negative) and the point form under-estimates instead of over-estimating — the sign of the quoted bias was wrong for exactly the class that triggered it.

The warning is not restricted to distributed conductors. Any geometry with segments long compared with the burial depth triggers it: a ring of radius 25 m at \(z = 0.8\) m discretised at \(\Delta s = 5\) m produces 4.909 m arc segments and warns ("by 26 %"), and so do coarse meshes, counterpoises, PEN conductors and surface-near tapes. Measured across the test suite it fires at 20 test nodes (≈35 emissions) in 9 modules.

It is also asymmetric across the engine family. mom, bem, cim, mom_sommerfeld and mutual inherit the identical biased diagonal from the shared kernel but do not emit the warning, so a compare_engines run warns for image and stays silent for the others while all of them carry the bias. The check stays in solve_image* because the shared kernel is re-entered per frequency and per excitation inside the current-solving backends, where the advisory would repeat several times per solve with a stacklevel pointing into solver internals rather than at user code.

Residual bias of the point-image diagonal (not fixed in 0.15.0)

The orientation-aware image-line term is still not used on the diagonal. The bias is real and large: a 10 m tape at \(z = 0.01\) m gives 54.361 / 35.491 / 26.603 / 22.039 / 21.043 Ω for \(\Delta s = 1.0 / 0.5 / 0.25 / 0.1 / 0.05\) m, i.e. +69 % at \(\Delta s = 0.5\) m against the converged ≈21 Ω.

Three points a reader should have straight:

1. Replacing \(1/(2z)\) by the horizontal image-line form unconditionally would make the vertical case worse, not better. At \(z = L = 0.5\) m the true vertical value is 1.021651, the point form 1.000000, the horizontal arsinh form 0.989866 — the point form is the closer of the two for a rod.
2. The expression \((2/L)\operatorname{arsinh}\bigl(L/(4z)\bigr)\) quoted above is the collocation value — the image line's potential at the segment midpoint. The reduction the rest of the pipeline uses since 0.15.0 (see "Node potential" below) is the Galerkin, surface-averaged one, $\bigl(2/L^2\bigr)\bigl[L \operatorname{arsinh}(L/2z)
3. \sqrt{L^2 + 4z^2} + 2z\bigr]$, which is 0.980575 at \(L = z = 0.5\) m against the collocation 0.989866. A future fix should use the Galerkin form for consistency with _weighted_node_potential.
4. The mechanical blocker previously recorded here — "the kernel is shared and receives no segment directions, so this needs a signature change for the whole family" — overstates the case. _Segment, _discretize_electrode and _self_corrected_kernel all live in solver/image.py, and an optional seg_directions=None keyword (point form when None) would leave every sibling backend bit-exact. What actually defers the fix is that switching the image family onto a different diagonal changes every stored reference impedance in the benchmark catalogue at once; that belongs in its own release step, not in a review-pass bug fix.

Until then the guards plus the warning delimit the validity envelope instead of hiding the bias.

Cluster constraints

Multiple electrodes connected by a Conductor form a galvanic cluster with a shared (unknown) cluster potential \(\varphi_c\) and a known total injected current \(I_{c,\text{in}} = \sum_{e\in c} I_{\text{src},e}\). The current sharing within the cluster is solved through the multi-port grounding matrix \(Z_{ij}\) (average potential at electrode \(i\) for unit current at electrode \(j\), reduced with the length-weighted average of the next subsection):

\[ \begin{bmatrix} Z & -C \\ C^{\top} & 0 \end{bmatrix} \begin{bmatrix} I \\ \varphi_c \end{bmatrix} = \begin{bmatrix} 0 \\ I_{\text{in}} \end{bmatrix}, \]

with \(C\) the cluster-membership indicator. The first \(N\) rows enforce \(\varphi_i = \varphi_c\) for every electrode in cluster \(c\); the last \(K\) rows enforce \(\sum_{i \in c} I_i = I_{c,\text{in}}\).

Node potential: the length-weighted (Galerkin) average

The "potential of electrode \(i\)" that enters the row reduction of \(Z_{ij}\) — and that is reported back as FieldResult.electrode_potentials — is the length-weighted average of the segment-midpoint potentials,

\[ \varphi_e \;=\; \frac{\sum_{k \in e} L_k\, \varphi_k}{\sum_{k \in e} L_k}, \]

which is the pairing dual to the uniform-per-unit-length current ansatz \(I_k = I_e L_k / \sum_j L_j\): it averages over the electrode surface, not over its segment list, and it restores exact discrete reciprocity \(Z_{ij} = Z_{ji}\). Consequently grounding_impedance(e) equals the diagonal \(Z_{ii}\) of the matrix the solver actually inverted, to round-off, and cluster_impedance(e) uses the cluster-wide weighted potential \(\varphi_c = \sum_e L_e \varphi_e / \sum_e L_e\) — a physical quantity, invariant under renaming of the members.

For an electrode with equal segment lengths — every plain rod, ring or strip — the weighted and unweighted averages coincide, so those results are unchanged. They differ wherever segment lengths are mixed: mesh electrodes with \(\Delta x/n_x \neq \Delta y/n_y\), polylines with legs of incommensurable length, and rods split at a soil-layer interface (ADR-0007). Up to 0.14.1 the reduction inside the solve was weighted (0.11.0, WP-B3) but the reported potential was a plain mean, so the impedance the user read back was not the impedance the solver enforced: measured \(-10.9\,\%\) on a 4.509 m rod split at \(h_1 = 5\) m (45.686 Ω enforced, 40.701 Ω reported), and ideally bonded electrodes reported different potentials for one and the same constrained node.

Postprocessing

After the cluster currents are known, every segment current is fixed by the uniform-per-unit-length rule. Field-point evaluations (profiles, contours, transferred potentials) reuse the same kernel \(1/r + 1/r'\) at the actual field point. The FieldResult.potential helper is a thin wrapper around this evaluation.

Validity envelope

Property	Range / value
Soil model	HomogeneousSoil only
Frequency	quasi-static, \(f < 1\,\text{kHz}\)
Wire radius	\(a \ll L_i\) (thin-wire)
Segment length	\(L_i \lesssim a_{\text{eq}} / 5\) for stable averaging
Burial depth (midpoints)	\(z > 0\) and \(2 z > a\) for every segment midpoint — hard ValueError otherwise, raised in the shared reaction-matrix kernel so mom / bem / cim / mom_sommerfeld / mutual reject the same geometries
Burial depth (extent)	\(z_i - \tfrac{1}{2} L_i \lvert e_{z,i}\rvert \ge 0\): no segment may cross the surface. Hard ValueError, needs the segment tangent and therefore fires in image / image_2layer only. An upper end exactly at \(z = 0\) is legal (ordinary driven rod)
Self-image bias	\(L_i \le 4 z_i\) for an unbiased self-image term; ShallowSegmentWarning otherwise, quantified with the horizontal-segment formula. Emitted by image / image_2layer only although the bias is shared by the whole family
Air boundary	insulating (Neumann at \(z = 0\))
Far-field	\(\varphi \to 0\) as $

Convergence and cost

• Discretisation error. The uniform-per-unit-length ansatz carries a \(\sim 4{-}5\,\%\) residual compared to the Sunde rod formula at the canonical 1.5 m / 5 mm rod, and shrinks to \(< 1\,\%\) at sub-centimetre segment lengths and short rods.
• Computational cost. \(O(N^2)\) matrix build for the cluster reaction matrix; \(O(K^3)\) for the constraint solve, where \(K\) is the cluster count (typically 1–3). For typical geometries with \(N \le 10^3\), the homogeneous engine completes in milliseconds.
• Numerical singularity. Distances below _MIN_DISTANCE = 1 mm are clamped at the floor to keep the kernel finite during plot evaluations near the wire axis, and off-diagonal reaction-matrix pairs that hit the clamp raise a warning (two distinct conductors overlap). The diagonal image term is no longer clamped: a vanishing image separation is a modelling error, not a numerical one, and is rejected outright (see "Self-action correction").

Cross-validation notes

Counterpart	Expected agreement	What is checked
Dwight 1936 closed forms	\(\le 10\,\%\)	rod / ring / mesh DC resistance
mom (Galerkin)	\(\le 2\,\%\)	same kernel, different test function
image_2layer at \(K = 0\)	bit-exact	layered family collapses to homogeneous
cim at \(n = 1\)	bit-exact	matrix-pencil fit returns \(P = 0\), kernel matches
mom_sommerfeld at \(n = 1\)	bit-exact	quadrature short-circuits to closed form
fem (axisymmetric volume PDE)	\(\le 10\,\%\)	reduction to equivalent hemisphere

These bounds are codified as parametric pytest fixtures; see tests/test_cross_engines.py and tests/test_cross_engines_extended.py.

References

• Sunde, E. D. (1968). Earth Conduction Effects in Transmission Systems, Dover. Chapter 2 — image-charge construction and average-potential method.
• Dwight, H. B. (1936). Calculation of resistances to ground. AIEE Transactions 55. Reference DC resistances for canonical geometries.
• Tagg, G. F. (1964). Earth Resistances, Pitman. The practitioner's reference for image methods.

Example

import groundfield as gf

soil = gf.HomogeneousSoil(resistivity=100.0)
world = gf.create_world(soil=soil)
gf.create_electrode(world, "rod", name="g1",
                    position=(0.0, 0.0, 0.0), length=1.5)
gf.create_source(world, attached_to="g1", magnitude=1.0)

engine = gf.create_engine(backend="image",
                          segment_length=0.05,
                          frequencies=[50.0])
result = world.solve(engine)
print(result.cluster_impedance("g1")[0])

API reference

image

Image-charge backend for homogeneous soil.

Computes the potential field of an arbitrary grounding system in a homogeneous half-space (resistivity \(\rho\), soil surface at \(z = 0\), \(z\) axis pointing into the soil) using the classical image-charge method.

Notes

A point current source \(I\) at \(r_s = (x_s, y_s, z_s)\) with \(z_s > 0\) (inside the soil) produces, in a homogeneous half-space with an insulating soil surface, the potential $$ \varphi(r) \;=\; \frac{\rho\, I}{4\pi}\, \Big(\frac{1}{|r - r_s|} + \frac{1}{|r - r_s'|}\Big), $$ with the image \(r_s' = (x_s, y_s, -z_s)\) mirrored at the soil surface. An extended electrode is discretised into \(N\) segments; each segment carries one point current source at its midpoint. The total current \(I_e\) of an electrode is distributed uniformly per unit length across its segments — a surprisingly good approximation for wire electrodes at low frequencies (cf. Sunde 1968, Tagg 1964).

The input impedance of an electrode is computed as the average of the potential on its own segment midpoints (average-potential method). For a single driven rod the backend reproduces the Sunde formula within a few per cent.

Further properties of this backend:

• Frequency-independent: in the quasi-static range \(f < 1\,\mathrm{kHz}\) the backend returns the same real solution per frequency. Complex extensions (Carson, frequency-dependent soil) come in later backends.
• Multiple electrodes: each electrode has its own total current (sum of the current sources attached to it). An electrode without a source carries zero current and acts purely as a passive observer.

References

.. [1] E. D. Sunde, Earth Conduction Effects in Transmission Systems, Dover, 1968. .. [2] G. F. Tagg, Earth Resistances, Pitman, 1964.

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

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,
    )

Related material

• ADR-0001 documents why this homogeneous engine sits at the root of the engine family.