Postprocess

The postprocess subpackage turns a raw :class:FieldResult into quantities the user actually needs: potential / EPR plots, a fitted rho-f standard model, and a vector fit that produces a closed-form SymPy expression compatible with :class:groundinsight.BusType.impedance_formula.

Plotting helpers

The standard plot suite (plot_potential_contour, plot_potential_profile, plot_potential_radial) is re-exported at the package level for convenience and is documented in Quickstart. Below we focus on the reduced model path that closes the field → network bridge.

Vector fitting and SymPy export

Physical background

For a passive, linear, time-invariant grounding cluster, the driving-point impedance \(Z(s)\) at the feed-in electrode is a rational function of the Laplace variable \(s = j\omega\). Under the \(f \le 1\,\mathrm{kHz}\) quasi-static assumption, the function is well approximated by a low-order partial-fraction expansion

\[ Z(s) \;\approx\; R_\infty \;+\; s\,L_\infty \;+\; \sum_{k=1}^{N_p}\,\frac{r_k}{s - p_k} \]

with poles \(p_k\) on the negative real axis (damped RC modes) and residues \(r_k\). Complex-conjugate pole pairs are admissible and produce damped resonant LC-like behaviour. The fit is constructed by the Vector Fitting algorithm of Gustavsen & Semlyen 1999, which iterates pole locations to a stable solution and is the de-facto standard for transmission-line and grounding modelling.

The output is a SymPy expression in a single free symbol s with complex-conjugate pole pairs combined into real second-order terms. With the symbolic substitution \(s \to j\,2\pi f\) the expression matches the groundinsight.BusType.impedance_formula parser, so a field-grade Z(s) can be turned into a network-grade BusType without any manual reformatting.

Example

import numpy as np
import groundfield as gf
from groundfield.postprocess.vector_fitting import (
    vector_fit, fit_to_sympy, rho_f_from_field_result,
)

# 1. Run any groundfield engine and obtain Z(f) at the feed cluster.
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)

frequencies = np.array([10.0, 50.0, 150.0, 250.0, 500.0, 1000.0])
engine = gf.create_engine(backend="image", frequencies=frequencies.tolist())
result = world.solve(engine)

# 2. Fit a rational Z(s) directly from the FieldResult.
fit = rho_f_from_field_result(result, electrode_name="g1", n_poles=3)

# 3. Render as a SymPy expression and extract the formula string
#    that groundinsight.BusType.impedance_formula consumes.
expr = fit_to_sympy(fit, decimals=6)
print(expr)            # printable rational expression
formula_str = str(expr)

A pre-tabulated \(Z(f)\) array can be fitted directly with :func:vector_fit(frequencies, Z_values, n_poles=3). The returned :class:VectorFitResult carries poles, residues, and an evaluate(frequencies) method for re-evaluation.

For the export of the fit as a BusType (JSON file or live groundinsight instance) see the IO reference and ADR-0008.

API reference — vector fitting

vector_fitting

Vector fitting and SymPy export for the rho-f reduced model.

This module turns a frequency response computed by groundfield into a closed-form rational impedance expression that :mod:groundinsight can consume as BusType.impedance_formula. It is the bridge between the field-grade reference computation (groundfield) and the reduced equivalent-network model (groundinsight).

Mathematical background

For a passive, linear, time-invariant grounding cluster, the driving-point impedance \(Z(s)\) at the feed-in electrode is a rational function of the Laplace variable \(s = j\omega\). Under the \(f \le 1\,\mathrm{kHz}\) assumption the function is well approximated by a low-order partial-fraction expansion

\[ Z(s) \;\approx\; R_\infty \;+\; s\,L_\infty \;+\; \sum_{k=1}^{N_p}\,\frac{r_k}{s - p_k}, \]

with poles \(p_k\) on the negative real axis (damped RC modes) and residues \(r_k\). Complex-conjugate pole pairs are admissible too — they produce damped resonant LC-like behaviour. The fit is constructed by the Vector Fitting algorithm of Gustavsen & Semlyen 1999, which iterates pole locations to a stable solution and is the de-facto standard for transmission-line and grounding modelling.

Implementation

This module implements a clean, dependency-free version of vector fitting in NumPy:

1. Initial pole guess: linearly spaced on the negative real axis between the smallest and largest sampled frequencies.
2. Pole relocation via Sanathanan-Koerner-style least squares with a shared denominator (the canonical Gustavsen iteration).
3. Residue solve given the converged poles.
4. Stability enforcement: any pole ending up in the right half plane is reflected.

The exported SymPy expression follows the groundinsight.BusType.impedance_formula convention: a single free symbol s and arbitrary numeric constants.

References

• Gustavsen, B. & Semlyen, A. (1999). Rational approximation of frequency domain responses by Vector Fitting. IEEE Trans. Power Delivery 14(3), 1052–1061.
• Gustavsen, B. (2006). Improving the pole relocating properties of vector fitting. IEEE Trans. Power Delivery 21(3), 1587–1592.

VectorFitConvergenceWarning

Bases: UserWarning

The pole-relocation loop did not converge.

Emitted when the relative pole displacement of the last relocation step is still above 1 % after n_iter iterations (audit 2026-07-08, WP-B4 — the historic implementation ran a fixed number of blind iterations with no convergence check). Inspect VectorFitResult.pole_shift / n_iter_used and consider more iterations, a different n_poles, or weighting="inverse_magnitude" for wide-dynamic-range data.

VectorFitResult dataclass

VectorFitResult(
    poles: np.ndarray,
    residues: np.ndarray,
    R_inf: float,
    L_inf: float,
    rms_error: float,
    fit_frequencies: np.ndarray,
    fit_values: np.ndarray,
    n_iter_used: int = 0,
    pole_shift: float = 0.0,
)

Result of a vector fit.

Attributes:

Name	Type	Description
poles	ndarray	Complex poles \(p_k\) (1-D array, length \(N_p\)). Stable poles have poles.real <= 0.
residues	ndarray	Complex residues \(r_k\) (same shape as poles).
R_inf	float	Real constant offset \(R_\infty\) (DC residual).
L_inf	float	Real proportional-to-\(s\) term \(L_\infty\) (high-frequency slope). 0 if the fit was performed without it.
rms_error	float	Root-mean-square residual fitting error in \(\Omega\), taken over the input frequency points.
fit_frequencies	ndarray	Frequencies (Hz) used for the fit, kept for diagnostics.
fit_values	ndarray	Original \(Z(s)\) values used for the fit.

evaluate

evaluate(frequencies: Sequence[float]) -> np.ndarray

Evaluate the fitted \(Z(s)\) at arbitrary frequencies.

Source code in src/groundfield/postprocess/vector_fitting.py

def evaluate(self, frequencies: Sequence[float]) -> np.ndarray:
    """Evaluate the fitted $Z(s)$ at arbitrary frequencies."""
    s = 2j * math.pi * np.asarray(frequencies, dtype=float)
    Z = np.full_like(s, self.R_inf, dtype=complex)
    Z = Z + self.L_inf * s
    for p, r in zip(self.poles, self.residues):
        Z = Z + r / (s - p)
    return Z

VectorFitUnderdeterminedWarning

Bases: UserWarning

Vector fit is under- or exactly-determined.

A vector fit with n_poles poles has 2 * n_poles + extras real unknowns (residues plus optional R_inf / L_inf); the real-imag stacking of an input of length N gives 2 N equations. n_poles >= N therefore leads to an under- or exactly-determined least-squares problem. The fit can succeed numerically but is effectively an identity interpolation that carries no information beyond the input samples.

A dedicated warning category lets users opt into / out of the diagnostic with warnings.simplefilter("once", VectorFitUnderdeterminedWarning).

fit_to_sympy

fit_to_sympy(fit: VectorFitResult, *, decimals: int = 6)

Convert a :class:VectorFitResult to a SymPy expression.

The resulting expression is

\[ Z(s) \;=\; R_\infty \;+\; L_\infty \cdot s \;+\; \sum_k\,\frac{r_k}{s - p_k}. \]

Complex-conjugate pole/residue pairs are combined into a single real second-order term to keep the formula compact and interpretable for groundinsight:

\[ \frac{r}{s - p} + \frac{r^*}{s - p^*} \;=\; \frac{2\,\Re(r)\,(s - \Re(p)) - 2\,\Im(r)\,\Im(p)} {(s - \Re(p))^2 + \Im(p)^2}. \]

The coefficients are rounded to decimals digits to keep the printed formula readable.

Returns a :class:sympy.Expr. Free symbol: s.

Source code in src/groundfield/postprocess/vector_fitting.py

def fit_to_sympy(fit: VectorFitResult, *, decimals: int = 6):
    """Convert a :class:`VectorFitResult` to a SymPy expression.

    The resulting expression is

    $$
    Z(s) \\;=\\; R_\\infty \\;+\\; L_\\infty \\cdot s \\;+\\;
    \\sum_k\\,\\frac{r_k}{s - p_k}.
    $$

    Complex-conjugate pole/residue pairs are combined into a single
    real second-order term to keep the formula compact and
    interpretable for ``groundinsight``:

    $$
    \\frac{r}{s - p} + \\frac{r^*}{s - p^*}
    \\;=\\; \\frac{2\\,\\Re(r)\\,(s - \\Re(p)) -
                  2\\,\\Im(r)\\,\\Im(p)}
                 {(s - \\Re(p))^2 + \\Im(p)^2}.
    $$

    The coefficients are rounded to ``decimals`` digits to keep
    the printed formula readable.

    Returns a :class:`sympy.Expr`. Free symbol: ``s``.
    """
    import sympy as sp

    s = sp.Symbol("s", complex=True)
    expr = sp.Float(fit.R_inf, decimals)
    if fit.L_inf != 0.0:
        expr = expr + sp.Float(fit.L_inf, decimals) * s

    # Bucket poles into real and complex-conjugate pairs. Vector
    # fitting's eigenvalue solve does not return exactly mirrored
    # imaginary parts; we use a loose relative tolerance and
    # **symmetrise** the matched pair (taking the average of |Im|)
    # so the resulting expression is rigorously real-valued for
    # real s.
    used = np.zeros(len(fit.poles), dtype=bool)
    rel_tol = 1e-3
    for k, (p, r) in enumerate(zip(fit.poles, fit.residues)):
        if used[k]:
            continue
        if abs(p.imag) < rel_tol * max(abs(p.real), 1.0):
            # Effectively real pole
            expr = expr + sp.Float(r.real, decimals) / (
                s - sp.Float(p.real, decimals)
            )
            used[k] = True
            continue
        # Look for the conjugate (loose tolerance, opposite-sign imag)
        conj_idx = -1
        best_score = float("inf")
        for k2 in range(k + 1, len(fit.poles)):
            if used[k2]:
                continue
            d_real = abs(fit.poles[k2].real - p.real) / max(abs(p.real), 1.0)
            d_imag = abs(fit.poles[k2].imag + p.imag) / max(abs(p.imag), 1.0)
            score = d_real + d_imag
            if d_real < rel_tol and d_imag < rel_tol and score < best_score:
                conj_idx = k2
                best_score = score
        if conj_idx < 0:
            # No conjugate found — fall back to a single complex term
            expr = expr + sp.nsimplify(r, rational=False) / (
                s - sp.nsimplify(p, rational=False)
            )
            used[k] = True
            continue
        # Symmetrise: take averages so the pair is exactly conjugate.
        used[k] = True
        used[conj_idx] = True
        p2 = fit.poles[conj_idx]
        r2 = fit.residues[conj_idx]
        Re_p_avg = 0.5 * (p.real + p2.real)
        Im_p_avg = 0.5 * (abs(p.imag) + abs(p2.imag))  # use positive imag for the canonical pair
        Re_r_avg = 0.5 * (r.real + r2.real)
        # The residue conjugate carries opposite imag; canonical pair
        # has +Im_r corresponding to +Im_p. If p.imag > 0, r.imag is
        # the "canonical" imag; otherwise flip.
        Im_r_avg = 0.5 * (abs(r.imag) + abs(r2.imag))
        if p.imag < 0:
            Im_r_avg *= -1.0  # match sign convention
        # Real second-order term:
        # 2·Re(r)·(s - Re(p)) - 2·Im(r)·Im(p)  /  (s - Re(p))² + Im(p)²
        Re_p_sym = sp.Float(Re_p_avg, decimals)
        Im_p_sym = sp.Float(Im_p_avg, decimals)
        Re_r_sym = sp.Float(Re_r_avg, decimals)
        Im_r_sym = sp.Float(Im_r_avg, decimals)
        num = 2 * Re_r_sym * (s - Re_p_sym) - 2 * Im_r_sym * Im_p_sym
        den = (s - Re_p_sym) ** 2 + Im_p_sym ** 2
        expr = expr + num / den

    expr = sp.simplify(expr)

    # Defensive guard: every rho-f formula must depend on ``s``. A
    # constant Z(rho, f) silently breaks the groundinsight BusType
    # extraction, so we surface that as a ``UserWarning`` here. This
    # cannot trigger under the current ``vector_fit`` because
    # ``n_poles >= 1`` is enforced — the guard catches programmatically
    # constructed ``VectorFitResult`` objects with zero poles and zero
    # ``L_inf``.
    if s not in expr.free_symbols:
        import warnings

        warnings.warn(
            "fit_to_sympy: the resulting expression has no dependency "
            f"on the frequency-domain variable s ({expr!r}). Downstream "
            "groundinsight BusType.impedance_formula consumers expect "
            "a frequency-dependent rho-f model. Re-run vector_fit with "
            "n_poles >= 1 and at least one non-zero residue, or pass "
            "include_L_inf=True if you genuinely want a constant + "
            "inductive shoulder.",
            UserWarning,
            stacklevel=2,
        )

    return expr

rho_f_from_field_result

rho_f_from_field_result(
    result,
    electrode_name: str,
    *,
    n_poles: int = 4,
    n_iter: int = 8,
    include_R_inf: bool = True,
    include_L_inf: bool = False,
    complex_poles: bool = True
) -> VectorFitResult

Fit \(Z(s)\) for one electrode of a :class:FieldResult.

Computes \(Z_k = U_k / I_k\) at every frequency in result.frequencies for result.electrode_potentials[electrode_name] and result.electrode_currents[electrode_name], then runs :func:vector_fit on the resulting series.

Parameters:

Name	Type	Description	Default
result		:class:groundfield.FieldResult produced by an Engine.solve call with at least 4–8 frequency samples.	required
electrode_name	str	Name of the electrode to extract.	required
n_poles	int	Forwarded to :func:vector_fit.	4
n_iter	int	Forwarded to :func:vector_fit.	4
include_R_inf	int	Forwarded to :func:vector_fit.	4
include_L_inf	int	Forwarded to :func:vector_fit.	4
complex_poles	int	Forwarded to :func:vector_fit.	4

Returns:

Type	Description
VectorFitResult

Source code in src/groundfield/postprocess/vector_fitting.py

def rho_f_from_field_result(
    result,
    electrode_name: str,
    *,
    n_poles: int = 4,
    n_iter: int = 8,
    include_R_inf: bool = True,
    include_L_inf: bool = False,
    complex_poles: bool = True,
) -> VectorFitResult:
    """Fit $Z(s)$ for one electrode of a :class:`FieldResult`.

    Computes
    $Z_k = U_k / I_k$
    at every frequency in ``result.frequencies`` for
    ``result.electrode_potentials[electrode_name]`` and
    ``result.electrode_currents[electrode_name]``, then runs
    :func:`vector_fit` on the resulting series.

    Parameters
    ----------
    result
        :class:`groundfield.FieldResult` produced by an
        ``Engine.solve`` call with at least 4–8 frequency samples.
    electrode_name
        Name of the electrode to extract.
    n_poles, n_iter, include_R_inf, include_L_inf, complex_poles
        Forwarded to :func:`vector_fit`.

    Returns
    -------
    VectorFitResult
    """
    if electrode_name not in result.electrode_potentials:
        raise KeyError(
            f"electrode '{electrode_name}' not in result; "
            f"available: {list(result.electrode_potentials)}"
        )
    freqs = np.asarray(result.frequencies, dtype=float)
    U = np.asarray(result.electrode_potentials[electrode_name], dtype=complex)
    I = np.asarray(result.electrode_currents[electrode_name], dtype=complex)
    if I.shape != freqs.shape:
        raise ValueError(
            "Inconsistent shapes between frequencies and electrode currents"
        )
    # Audit 2026-07-08, WP-B4 (finding N23): frequencies where the
    # electrode carries no current (passive electrode / dead column)
    # are MASKED instead of being injected as Z = 0 samples — a
    # single dead frequency used to pull the whole fit through the
    # origin silently.
    alive = np.abs(I) > 0.0
    if not alive.any():
        raise ValueError(
            f"rho_f_from_field_result: electrode '{electrode_name}' "
            "carries no current at any frequency — Z = U/I is "
            "undefined. Fit the driven electrode instead."
        )
    if not alive.all():
        import warnings as _warnings

        _warnings.warn(
            f"rho_f_from_field_result: electrode '{electrode_name}' "
            f"carries no current at {int((~alive).sum())} of "
            f"{alive.size} frequencies; those samples are excluded "
            "from the fit (historic behaviour injected Z = 0 there).",
            UserWarning,
            stacklevel=2,
        )
    Z = U[alive] / I[alive]
    return vector_fit(
        freqs[alive], Z,
        n_poles=n_poles, n_iter=n_iter,
        include_R_inf=include_R_inf, include_L_inf=include_L_inf,
        complex_poles=complex_poles,
    )

vector_fit

vector_fit(
    frequencies: Sequence[float],
    Z_values: Sequence[complex],
    *,
    n_poles: int = 4,
    n_iter: int = 8,
    include_R_inf: bool = True,
    include_L_inf: bool = False,
    complex_poles: bool = True,
    weighting: str = "uniform"
) -> VectorFitResult

Fit a rational \(Z(s)\) approximation to a sampled frequency response.

Implements the Vector Fitting algorithm of Gustavsen & Semlyen 1999 in its single-output, equal-weight form. The denominator polynomial is shared between numerator and denominator estimation (Sanathanan-Koerner with iterative pole relocation); after n_iter outer iterations the poles are taken as fixed and the residues are solved from a final linear system.

Parameters:

Name	Type	Description	Default
frequencies	Sequence[float]	Frequency samples in Hz, length \(N\).	required
Z_values	Sequence[complex]	Complex impedance samples \(Z(j2\pi f_k)\) in \(\Omega\), length \(N\).	required
n_poles	int	Target number of poles. Complex-conjugate pairs count as 2 poles each.	4
n_iter	int	Number of pole-relocation iterations.	8
include_R_inf	bool	Whether to fit a constant offset \(R_\infty\) (DC residual).	True
include_L_inf	bool	Whether to fit a \(s\cdot L_\infty\) proportional term.	False
complex_poles	bool	True (default) initialises with complex-conjugate pole pairs, suitable for resonant behaviour. Set False for purely-resistive / monotonic responses.	True
weighting	str	"uniform" (default, historic behaviour) or "inverse_magnitude": weights every sample by 1/|Z_n|, the standard remedy when |Z| spans several decades and the unweighted LS would fit only the largest magnitudes (audit 2026-07-08, WP-B4).	'uniform'

Notes

Since the 2026-07-08 audit (WP-B4) the least-squares systems use the conjugate-pair real basis (see :func:_pair_basis): complex pole pairs carry genuinely complex residues, matching the cited Gustavsen & Semlyen algorithm. The pole-relocation loop additionally monitors the relative pole displacement, breaks early below 1e-8 and warns (:class:VectorFitConvergenceWarning) when it is still above 1 % after n_iter iterations; the diagnostics are exposed as VectorFitResult.n_iter_used / pole_shift.

Returns:

Type	Description
VectorFitResult

Source code in src/groundfield/postprocess/vector_fitting.py

def vector_fit(
    frequencies: Sequence[float],
    Z_values: Sequence[complex],
    *,
    n_poles: int = 4,
    n_iter: int = 8,
    include_R_inf: bool = True,
    include_L_inf: bool = False,
    complex_poles: bool = True,
    weighting: str = "uniform",
) -> VectorFitResult:
    """Fit a rational $Z(s)$ approximation to a sampled frequency response.

    Implements the Vector Fitting algorithm of Gustavsen & Semlyen
    1999 in its single-output, equal-weight form. The denominator
    polynomial is shared between numerator and denominator
    estimation (Sanathanan-Koerner with iterative pole relocation);
    after ``n_iter`` outer iterations the poles are taken as fixed
    and the residues are solved from a final linear system.

    Parameters
    ----------
    frequencies
        Frequency samples in Hz, length $N$.
    Z_values
        Complex impedance samples $Z(j2\\pi f_k)$ in $\\Omega$,
        length $N$.
    n_poles
        Target number of poles. Complex-conjugate pairs count as
        2 poles each.
    n_iter
        Number of pole-relocation iterations.
    include_R_inf
        Whether to fit a constant offset $R_\\infty$ (DC residual).
    include_L_inf
        Whether to fit a $s\\cdot L_\\infty$ proportional term.
    complex_poles
        ``True`` (default) initialises with complex-conjugate pole
        pairs, suitable for resonant behaviour. Set ``False`` for
        purely-resistive / monotonic responses.
    weighting
        ``"uniform"`` (default, historic behaviour) or
        ``"inverse_magnitude"``: weights every sample by
        ``1/|Z_n|``, the standard remedy when ``|Z|`` spans several
        decades and the unweighted LS would fit only the largest
        magnitudes (audit 2026-07-08, WP-B4).

    Notes
    -----
    Since the 2026-07-08 audit (WP-B4) the least-squares systems use
    the **conjugate-pair real basis** (see :func:`_pair_basis`):
    complex pole pairs carry genuinely complex residues, matching
    the cited Gustavsen & Semlyen algorithm. The pole-relocation
    loop additionally monitors the relative pole displacement,
    breaks early below ``1e-8`` and warns
    (:class:`VectorFitConvergenceWarning`) when it is still above
    1 % after ``n_iter`` iterations; the diagnostics are exposed as
    ``VectorFitResult.n_iter_used`` / ``pole_shift``.

    Returns
    -------
    VectorFitResult
    """
    frequencies = np.asarray(frequencies, dtype=float)
    Z = np.asarray(Z_values, dtype=complex)
    if frequencies.shape != Z.shape:
        raise ValueError(
            "frequencies and Z_values must have the same shape"
        )
    if n_poles < 1:
        raise ValueError(
            "vector_fit: n_poles must be at least 1, got "
            f"n_poles={n_poles}. n_poles == 0 produces an s-free fit "
            "whose SymPy conversion yields a constant Z(rho, f) "
            "expression without frequency dependence — that is never "
            "the intended behaviour, so this case is rejected at the "
            "API boundary."
        )
    # The previous trigger
    # ``2 * n_poles >= N`` counted 4 real DOFs per pole regardless of
    # whether the search was constrained to conjugate-symmetric pairs.
    # Under ``complex_poles=True`` (the default) the conjugate-symmetry
    # constraint halves the actual free parameters: a pair carries
    # 4 real DOFs (real pole, imag pole, real residue, imag residue),
    # and a lone unpaired real pole carries 2. The corrected formula
    # counts ``n_independent_poles = n_poles // 2 + (n_poles % 2)`` for
    # the pair-constrained search and ``n_poles`` for the all-real one.
    # The strict ``>`` in the conjugate branch keeps a uniquely
    # determined fit (e.g. ``n_poles=2, N=2``) from emitting a false
    # positive.
    if complex_poles:
        n_independent_poles = n_poles // 2 + (n_poles % 2)
        is_underdetermined = 2 * n_independent_poles > frequencies.size
    else:
        is_underdetermined = 2 * n_poles >= frequencies.size
    if is_underdetermined:
        import warnings as _warnings

        if complex_poles:
            ratio_explanation = (
                "2 * n_independent_poles > len(frequencies) is the "
                "critical ratio under complex-conjugate-pair fitting"
            )
        else:
            ratio_explanation = (
                "2 * n_poles >= len(frequencies) is the critical "
                "ratio under all-real-pole fitting"
            )
        _warnings.warn(
            "vector_fit: n_poles="
            f"{n_poles} with len(frequencies)={frequencies.size} "
            "leads to an under- or exactly-determined least-squares "
            "fit (each pole contributes residue + pole-location "
            f"DOFs; {ratio_explanation}). The numerical solve may "
            "succeed but produces an identity interpolation of the "
            "input samples that carries no information beyond the "
            "raw frequency response. Reduce n_poles to "
            f"{max(1, frequencies.size // 2 - 1)} or fewer, or supply "
            "more frequency samples.",
            VectorFitUnderdeterminedWarning,
            stacklevel=2,
        )
    if weighting not in ("uniform", "inverse_magnitude"):
        raise ValueError(
            "vector_fit: weighting must be 'uniform' or "
            f"'inverse_magnitude', got {weighting!r}."
        )
    s = 2j * math.pi * frequencies
    if weighting == "inverse_magnitude":
        w = 1.0 / np.maximum(np.abs(Z), 1e-12 * float(np.max(np.abs(Z))))
    else:
        w = np.ones(s.size)

    poles = _initial_poles(frequencies, n_poles, complex_conj=complex_poles)
    poles = _enforce_conjugate_symmetry(poles)

    # Outer loop: pole relocation. The system is, schematically,
    #   sum_k r_k / (s_n - p_k) + d + s_n*h - sum_k r'_k Z_n / (s_n - p_k)
    #     = Z_n
    # where {r_k, d, h, r'_k} are the unknowns and the new poles are
    # the zeros of the denominator polynomial 1 + sum r'_k/(s-p_k).
    # Both the numerator residues r_k and the sigma residues r'_k use
    # the conjugate-pair real basis (audit 2026-07-08, WP-B4), so
    # complex pole pairs carry genuinely complex residues.
    extras = (1 if include_R_inf else 0) + (1 if include_L_inf else 0)
    n_iter_used = 0
    pole_shift = float("inf")
    for _ in range(max(1, n_iter)):
        structure = _pole_pair_structure(poles)
        basis = _pair_basis(s, poles, structure)      # (N, n_poles)
        A = np.zeros((s.size, n_poles + extras + n_poles), dtype=complex)
        A[:, :n_poles] = basis
        col = n_poles
        if include_R_inf:
            A[:, col] = 1.0
            col += 1
        if include_L_inf:
            A[:, col] = s
            col += 1
        A[:, col:col + n_poles] = -Z[:, None] * basis

        # Solve with real-imag stacking (least-squares), optionally
        # row-weighted.
        Aw = A * w[:, None]
        M = np.vstack([Aw.real, Aw.imag])
        rhs = np.concatenate([(w * Z).real, (w * Z).imag])
        sol, *_ = np.linalg.lstsq(M, rhs, rcond=None)

        # Reconstruct the complex sigma residues from the pair basis.
        sigma_residues = _coeffs_to_residues(
            sol[n_poles + extras:], poles, structure,
        )
        # New poles: eigenvalues of A_p - b·c^T where A_p = diag(poles),
        # b = ones, c = sigma_residues. (Standard Gustavsen formulation.)
        Ap = np.diag(poles)
        b = np.ones(n_poles)
        new_poles = np.linalg.eigvals(Ap - np.outer(b, sigma_residues))
        new_poles = _enforce_stability(new_poles)
        new_poles = _enforce_conjugate_symmetry(new_poles)
        # Sort for deterministic order (deepest |imag| first).
        sorted_idx = np.argsort(-np.abs(new_poles.imag))
        new_poles = new_poles[sorted_idx]
        n_iter_used += 1
        scale = float(np.max(np.abs(poles))) or 1.0
        pole_shift = float(
            np.max(np.abs(np.sort_complex(new_poles)
                          - np.sort_complex(poles)))
        ) / scale
        poles = new_poles
        if pole_shift < 1e-8:
            break
    if pole_shift > 1e-2:
        import warnings as _warnings

        _warnings.warn(
            f"vector_fit: pole relocation did not converge within "
            f"n_iter={n_iter} iterations (last relative pole shift "
            f"{pole_shift:.2e}). The returned fit may sit in a limit "
            "cycle; inspect rms_error, increase n_iter, or adjust "
            "n_poles / weighting.",
            VectorFitConvergenceWarning,
            stacklevel=2,
        )

    # Final residue solve with the converged poles (pair basis —
    # residues at conjugate poles are exact conjugates by
    # construction, no post-hoc repair required).
    structure = _pole_pair_structure(poles)
    basis = _pair_basis(s, poles, structure)
    A = np.zeros((s.size, n_poles + extras), dtype=complex)
    A[:, :n_poles] = basis
    col = n_poles
    if include_R_inf:
        A[:, col] = 1.0
        col += 1
    if include_L_inf:
        A[:, col] = s
        col += 1
    Aw = A * w[:, None]
    M = np.vstack([Aw.real, Aw.imag])
    rhs = np.concatenate([(w * Z).real, (w * Z).imag])
    sol, *_ = np.linalg.lstsq(M, rhs, rcond=None)
    residues = _coeffs_to_residues(sol[:n_poles], poles, structure)
    R_inf = float(sol[n_poles]) if include_R_inf else 0.0
    L_inf = float(sol[n_poles + (1 if include_R_inf else 0)]) if include_L_inf else 0.0

    Z_fit = (
        np.full_like(s, R_inf, dtype=complex)
        + L_inf * s
        + np.sum(residues[:, None] / (s[None, :] - poles[:, None]), axis=0)
    )
    rms = float(np.sqrt(np.mean(np.abs(Z_fit - Z) ** 2)))

    return VectorFitResult(
        n_iter_used=n_iter_used,
        pole_shift=pole_shift,
        poles=poles,
        residues=residues,
        R_inf=R_inf,
        L_inf=L_inf,
        rms_error=rms,
        fit_frequencies=frequencies,
        fit_values=Z,
    )

Standard rho-f form

The five-coefficient ansatz

\[ Z(\rho, f) \;=\; k_1 \rho \;+\; (k_2 + j k_3)\,f \;+\; (k_4 + j k_5)\,f\,\rho \]

is fitted across a \(\rho\)-sweep and is already in the canonical (rho, f) symbol set used by groundinsight.BusType.impedance_formula. See :mod:groundfield.postprocess.rho_f_standard and the IO reference for the export path.

rho_f_standard

Standard-form rho-f model for the groundinsight bridge.

The reduced grounding-cluster impedance uses the physically-motivated 5-coefficient form

\[ Z(\rho, f) \;=\; k_1 \rho \;+\; (k_2 + j k_3)\,f \;+\; (k_4 + j k_5)\,f\,\rho, \]

with \(\rho\) a soil-resistivity parameter (typically the upper-layer \(\rho_1\) in a fixed-structure 2-layer setup) and \(f\) the frequency. The five real coefficients have direct physical interpretation:

• \(k_1\,\rho\) — DC spreading resistance (Dwight-class scaling with the dominant local soil resistivity).
• \(k_2\,f\) — purely-inductive coupling that does not depend on the soil (e.g. a metallic-cable loop-inductance term).
• \(k_3\,f\) — purely-resistive frequency-dependent term (negligible in most quasi-static typical studies).
• \(k_4\,f\,\rho\) — Carson-type earth-return resistance: scales with both frequency and soil resistivity.
• \(k_5\,f\,\rho\) — Carson-type earth-return reactance.

This is not a general rational function — it is a fixed parametric ansatz that captures the leading orders for production-grade grounding-cluster impedances. It is fitted from a parametric family of groundfield runs that span both \(\rho\) and \(f\), and exported as a SymPy expression with two free symbols \(\rho\) (rho) and \(f\) (f) — the BusType.impedance_formula convention used by groundinsight.

Mathematically the fit is a linear least-squares problem in the five real unknowns:

• Real part: \(\Re Z = k_1\rho + k_2 f + k_4 f\rho\) → 3-feature regression in (\(\rho\), \(f\), \(f\rho\)).
• Imaginary part: \(\Im Z = k_3 f + k_5 f\rho\) → 2-feature regression in (\(f\), \(f\rho\)).

Both halves are decoupled in the coefficients, so the fit is unique whenever the sample set spans at least two distinct \(\rho\) values and at least two distinct frequencies.

References

• The rho-f model is the reduced grey-box representation handed to groundinsight.
• groundinsight.BusType.impedance_formula: the consumer of the SymPy expression returned by :func:fit_to_sympy_standard.

RhoFStandardFit dataclass

RhoFStandardFit(
    k1: float,
    k2: float,
    k3: float,
    k4: float,
    k5: float,
    rms_error: float,
    rms_relative: float,
    sample_rho: np.ndarray,
    sample_f: np.ndarray,
    sample_Z: np.ndarray,
)

Result of a standard-form rho-f fit.

Attributes:

Name	Type	Description
k1, k2, k3, k4, k5		The five real coefficients of the formula \(Z = k_1\rho + (k_2 + j k_3)f + (k_4 + j k_5)f\rho\).
rms_error	float	Root-mean-square residual error in \(\Omega\) over the input samples.
rms_relative	float	rms_error divided by the mean of \(|Z|\) over the samples — a dimensionless quality figure.
sample_rho, sample_f, sample_Z		Original samples used for the fit (1-D arrays of length \(N\)), kept for diagnostics.

coefficients property

coefficients: dict[str, float]

Return the five coefficients as a dict.

evaluate

evaluate(
    rho: float | np.ndarray, f: float | np.ndarray
) -> complex | np.ndarray

Evaluate the fitted \(Z(\rho, f)\) at arbitrary points.

Source code in src/groundfield/postprocess/rho_f_standard.py

def evaluate(
    self, rho: float | np.ndarray, f: float | np.ndarray,
) -> complex | np.ndarray:
    """Evaluate the fitted $Z(\\rho, f)$ at arbitrary points."""
    rho = np.asarray(rho, dtype=float)
    f = np.asarray(f, dtype=float)
    return (
        self.k1 * rho
        + (self.k2 + 1j * self.k3) * f
        + (self.k4 + 1j * self.k5) * f * rho
    )

fit_rho_f_standard

fit_rho_f_standard(
    rho_samples: Sequence[float],
    f_samples: Sequence[float],
    Z_samples: Sequence[complex],
) -> RhoFStandardFit

Fit the 5-coefficient standard rho-f form via linear LSQ.

Parameters:

Name	Type	Description	Default
rho_samples	Sequence[float]	Soil-resistivity values \(\rho\) in \(\Omega\,\mathrm{m}\) for each sample (length \(N\)).	required
f_samples	Sequence[float]	Frequencies \(f\) in Hz for each sample (length \(N\)).	required
Z_samples	Sequence[complex]	Complex driving-point impedances \(Z\) in \(\Omega\) at each \((\rho, f)\) sample (length \(N\)).	required

Returns:

Type	Description
RhoFStandardFit

Raises:

Type	Description
ValueError	If sample arrays have inconsistent lengths, or if there are fewer than two distinct \(\rho\) values or fewer than two distinct \(f\) values (the LSQ is then under-determined).

Source code in src/groundfield/postprocess/rho_f_standard.py

def fit_rho_f_standard(
    rho_samples: Sequence[float],
    f_samples: Sequence[float],
    Z_samples: Sequence[complex],
) -> RhoFStandardFit:
    """Fit the 5-coefficient standard rho-f form via linear LSQ.

    Parameters
    ----------
    rho_samples
        Soil-resistivity values $\\rho$ in $\\Omega\\,\\mathrm{m}$
        for each sample (length $N$).
    f_samples
        Frequencies $f$ in Hz for each sample (length $N$).
    Z_samples
        Complex driving-point impedances $Z$ in $\\Omega$ at each
        $(\\rho, f)$ sample (length $N$).

    Returns
    -------
    RhoFStandardFit

    Raises
    ------
    ValueError
        If sample arrays have inconsistent lengths, or if there are
        fewer than two distinct $\\rho$ values *or* fewer than two
        distinct $f$ values (the LSQ is then under-determined).
    """
    rho = np.asarray(rho_samples, dtype=float)
    f = np.asarray(f_samples, dtype=float)
    Z = np.asarray(Z_samples, dtype=complex)
    if rho.shape != f.shape or rho.shape != Z.shape:
        raise ValueError(
            "rho_samples, f_samples and Z_samples must have identical shape; "
            f"got rho={rho.shape}, f={f.shape}, Z={Z.shape}."
        )
    if rho.size < 4:
        raise ValueError(
            "fit_rho_f_standard needs at least 4 samples; got "
            f"{rho.size}."
        )
    if np.unique(rho).size < 2:
        raise ValueError(
            "fit_rho_f_standard needs at least two distinct rho values "
            "to identify k1 and k4 separately."
        )
    if np.unique(f).size < 2:
        raise ValueError(
            "fit_rho_f_standard needs at least two distinct f values "
            "to identify k2 and k4 separately."
        )

    # Real part: Re(Z) = k1·ρ + k2·f + k4·f·ρ
    A_real = np.column_stack([rho, f, f * rho])
    sol_real, *_ = np.linalg.lstsq(A_real, Z.real, rcond=None)
    k1, k2, k4 = float(sol_real[0]), float(sol_real[1]), float(sol_real[2])

    # Imag part: Im(Z) = k3·f + k5·f·ρ
    A_imag = np.column_stack([f, f * rho])
    sol_imag, *_ = np.linalg.lstsq(A_imag, Z.imag, rcond=None)
    k3, k5 = float(sol_imag[0]), float(sol_imag[1])

    Z_fit = (
        k1 * rho
        + (k2 + 1j * k3) * f
        + (k4 + 1j * k5) * f * rho
    )
    rms = float(np.sqrt(np.mean(np.abs(Z_fit - Z) ** 2)))
    rms_rel = rms / max(float(np.mean(np.abs(Z))), 1e-12)

    return RhoFStandardFit(
        k1=k1, k2=k2, k3=k3, k4=k4, k5=k5,
        rms_error=rms, rms_relative=rms_rel,
        sample_rho=rho, sample_f=f, sample_Z=Z,
    )

fit_to_sympy_standard

fit_to_sympy_standard(
    fit: RhoFStandardFit, *, decimals: int = 6
)

Convert a :class:RhoFStandardFit to a SymPy expression.

Returns a :class:sympy.Expr in two free symbols rho and f (both real), in the canonical typical form

\[ Z(\rho, f) \;=\; k_1\rho \;+\; (k_2 + j k_3)\,f \;+\; (k_4 + j k_5)\,f\,\rho. \]

The expression is suitable for direct insertion into groundinsight.BusType.impedance_formula.

Source code in src/groundfield/postprocess/rho_f_standard.py

def fit_to_sympy_standard(fit: RhoFStandardFit, *, decimals: int = 6):
    """Convert a :class:`RhoFStandardFit` to a SymPy expression.

    Returns a :class:`sympy.Expr` in two free symbols ``rho`` and
    ``f`` (both real), in the canonical typical form

    $$
    Z(\\rho, f) \\;=\\; k_1\\rho \\;+\\; (k_2 + j k_3)\\,f
                      \\;+\\; (k_4 + j k_5)\\,f\\,\\rho.
    $$

    The expression is suitable for direct insertion into
    ``groundinsight.BusType.impedance_formula``.
    """
    import sympy as sp

    rho = sp.Symbol("rho", real=True, positive=True)
    f = sp.Symbol("f", real=True, positive=True)
    j = sp.I
    k1 = sp.Float(fit.k1, decimals)
    k2 = sp.Float(fit.k2, decimals)
    k3 = sp.Float(fit.k3, decimals)
    k4 = sp.Float(fit.k4, decimals)
    k5 = sp.Float(fit.k5, decimals)
    expr = k1 * rho + (k2 + j * k3) * f + (k4 + j * k5) * f * rho
    return sp.simplify(expr)

rho_f_standard_from_results

rho_f_standard_from_results(
    results: Sequence,
    rhos: Sequence[float],
    electrode_name: str,
) -> RhoFStandardFit

Build the (ρ, f, Z) sample table from a list of FieldResults.

Use case: parametric soil-resistivity sweep. Run one Engine.solve per soil resistivity, collect the FieldResults along with the driving \(\rho\), and pass the lot to this function.

Parameters:

Name	Type	Description	Default
results	Sequence	List of :class:groundfield.FieldResult (one per \(\rho\) value).	required
rhos	Sequence[float]	Soil-resistivity parameter \(\rho\) corresponding to each FieldResult (same length).	required
electrode_name	str	Name of the electrode to extract.	required

Returns:

Type	Description
RhoFStandardFit

Source code in src/groundfield/postprocess/rho_f_standard.py

def rho_f_standard_from_results(
    results: Sequence,
    rhos: Sequence[float],
    electrode_name: str,
) -> RhoFStandardFit:
    """Build the (ρ, f, Z) sample table from a list of FieldResults.

    Use case: parametric soil-resistivity sweep. Run one
    ``Engine.solve`` per soil resistivity, collect the FieldResults
    along with the driving $\\rho$, and pass the lot to this
    function.

    Parameters
    ----------
    results
        List of :class:`groundfield.FieldResult` (one per
        $\\rho$ value).
    rhos
        Soil-resistivity parameter $\\rho$ corresponding to each
        FieldResult (same length).
    electrode_name
        Name of the electrode to extract.

    Returns
    -------
    RhoFStandardFit
    """
    if len(results) != len(rhos):
        raise ValueError(
            f"results and rhos must have same length; "
            f"got {len(results)} vs {len(rhos)}."
        )
    rho_arr: list[float] = []
    f_arr: list[float] = []
    Z_arr: list[complex] = []
    for res, rho_val in zip(results, rhos):
        if electrode_name not in res.electrode_potentials:
            raise KeyError(
                f"electrode '{electrode_name}' not in FieldResult; "
                f"available: {list(res.electrode_potentials)}"
            )
        U = np.asarray(res.electrode_potentials[electrode_name], dtype=complex)
        I = np.asarray(res.electrode_currents[electrode_name], dtype=complex)
        f_local = np.asarray(res.frequencies, dtype=float)
        Z = np.where(np.abs(I) > 0.0, U / np.where(I != 0, I, 1.0), 0.0 + 0.0j)
        rho_arr.extend([float(rho_val)] * len(f_local))
        f_arr.extend(f_local.tolist())
        Z_arr.extend(Z.tolist())
    return fit_rho_f_standard(rho_arr, f_arr, Z_arr)

Touch and step voltages

The :mod:groundfield.postprocess.safety module turns a :class:FieldResult into the engineering safety quantities used in EN 50522:2010 / IEC 61936-1.

Physical background

The earth potential rise of a touched grounding cluster is

\[ U_E \;=\; \varphi_\text{cluster}. \]

Standing on the soil surface 1 m away from the touched part, the person's feet sit at the surface potential \(\varphi(\mathbf r_\text{feet})\) and the voltage appearing across the body is

\[ U_T \;=\; U_E - \varphi(\mathbf r_\text{feet}). \]

The step voltage between two surface points at the typical \(d_\text{step} = 1\,\mathrm{m}\) separation is

\[ U_S \;=\; \varphi(\mathbf r_1) - \varphi(\mathbf r_1 + d_\text{step}\,\hat{\mathbf e}). \]

Both quantities are returned as complex phasors per frequency index, so that inductive- and resistive-coupling effects above DC remain visible in typical studies.

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.5), length=1.5)
gf.create_source(world, attached_to="g1", magnitude=10.0)

result = gf.create_engine(backend="image", segment_length=0.05).solve(world)

# Touch voltage 1 m east of the rod:
U_T = gf.touch_voltage(result, world, electrode="g1", distance=1.0)

# Worst-case touch voltage on a circle around the rod:
angles, voltages = gf.touch_voltage_envelope(
    result, world, electrode="g1", distance=1.0, n_angles=24,
)
U_T_worst = float(abs(voltages).max())

# Step voltage 1 m east of the rod, walking radially outward:
U_S = gf.step_voltage(result, position=(1.0, 0.0, 0.0),
                      direction=(1.0, 0.0, 0.0), step=1.0)

# Reference: EN 50522 Table B.4 permissible touch voltage at t_F = 0.5 s.
U_TP = gf.permissible_touch_voltage_en50522(0.5)   # 225 V

Validity envelope

• Frequency: quasi-static envelope \(f \le 1\,\mathrm{kHz}\), inherited from the underlying Green's function.
• Backends: every solver that populates point_sources (image, image_2layer, image_nlayer, mom, mom_sommerfeld, cim, bem). Stub backends raise inside :meth:FieldResult.potential.
• Coordinate convention: soil surface at \(z = 0\), positive \(z\) points downwards into the soil. The default surface_z = 0.0 models bare-foot contact on the ground surface.
• permissible_touch_voltage_en50522 reproduces EN 50522:2010 Table B.4 verbatim — eight anchors at \(t_f \in \{0.05, 0.10, 0.20, 0.50, 1.00, 2.00, 5.00, 10.00\}\,\mathrm{s}\) with \(U_{TP} \in \{725, 655, 525, 225, 115, 95, 85, 85\}\,\mathrm{V}\), values rounded to 5 V in the standard. Log-log interpolation inside the grid; outside the grid the values are clamped to the table endpoints (no relaxation below 50 ms; the terminal 85 V plateau between 5 s and 10 s is reproduced exactly).

API reference — safety

safety

Touch and step voltages for grounding-system safety assessment.

This module turns a :class:groundfield.solver.result.FieldResult into the engineering quantities used in the grounding-system safety verification according to EN 50522:2010 / IEC 61936-1: the touch voltage :math:U_T and the step voltage :math:U_S. It also provides the corresponding permissible values :math:U_{TP}(t) from EN 50522:2010, Figure B.3.

Mathematical / physical content

A grounding cluster at earth potential rise

.. math::

U_E \;=\; \varphi_\text{cluster},

is connected — through any galvanically bonded metallic part — to whatever a person can touch. Standing on the soil surface 1 m away from the touched part (the conventional test point of EN 50522), the person's feet sit at the surface potential :math:\varphi(\mathbf{r}_\text{feet}). The voltage appearing across the body is

.. math::

U_T \;=\; U_E \;-\; \varphi(\mathbf{r}_\text{feet}).

For a step voltage the body bridges two surface points at the typical step distance :math:d_\text{step} = 1\,\mathrm{m}:

.. math::

U_S \;=\; \varphi(\mathbf{r}_1) \;-\;
          \varphi(\mathbf{r}_1 + d_\text{step}\,\hat{\mathbf e}).

Both quantities are returned as complex phasors per frequency index — consistent with the rest of groundfield — so that inductive- and resistive-coupling effects above DC remain visible.

Validity envelope

• Frequency: quasi-static envelope :math:f \le 1\,\mathrm{kHz}, inherited from the underlying Green's function.
• Geometry: relies on :meth:FieldResult.potential and :meth:FieldResult.electrode_potentials. Works for every backend that populates point_sources (image, image_2layer, image_nlayer, mom, mom_sommerfeld, cim, bem) and is silently inapplicable to stub backends (raises through :meth:FieldResult.potential).
• Coordinate convention: the soil surface is at :math:z = 0, positive :math:z points downwards into the soil (see :mod:groundfield.geometry.electrodes). The default surface_z = 0.0 therefore models bare-foot contact at the ground surface.
• Permissible-voltage helper :func:permissible_touch_voltage_en50522 covers the low-voltage standard curve (EN 50522:2010, Fig. B.3) for fault-clearing times :math:50\,\mathrm{ms} \le t_F \le 10\, \mathrm{s}. Outside that range the value is clamped to the table endpoints, mirroring the stationary plateau the standard prescribes for long-duration faults.

References

• EN 50522:2010 — Earthing of power installations exceeding 1 kV a.c., Annex B (allowable touch voltages).
• IEC 61936-1:2014 — Power installations exceeding 1 kV a.c.
• IEEE Std 80-2013 — Guide for safety in AC substation grounding.

permissible_touch_voltage_en50522

permissible_touch_voltage_en50522(
    t_clear_s: float,
) -> float

Permissible touch voltage :math:U_{TP}(t_F) per EN 50522:2010.

Returns the maximum admissible touch voltage as a function of the fault-clearing time :math:t_F. The reference values are taken verbatim from EN 50522:2010, Table B.4 — the normative tabulation of the curve in Figure B.3 (values rounded to 5 V in the standard).

Mathematically the helper is a piecewise-loglog interpolation over the canonical anchor points

.. math::

\{(t_k, U_{TP,k})\} \;=\;
\{(0.05, 725),\ (0.10, 655),\ (0.20, 525),\
(0.50, 225),\ (1.00, 115),\ (2.00, 95),\
(5.00, 85),\ (10.0, 85)\},

with units :math:\mathrm{s} and :math:\mathrm{V}. Outside the grid the values are clamped to the endpoints — the standard does not extend the curve below 50 ms, and the terminal plateau at 85 V is explicit (the table reports identical values at 5 s and 10 s).

Parameters:

Name	Type	Description	Default
t_clear_s	float	Fault-clearing time :math:t_F in seconds. Must be > 0.	required

Returns:

Type	Description
float	Permissible touch voltage :math:U_{TP} in V.

Raises:

Type	Description
ValueError	If t_clear_s is not strictly positive.

Notes

The EN 50522 curve is intended for engineering design; this helper is used as a reference line on plots of computed touch voltages, not as a regulatory pass/fail. Use :func:touch_voltage to compute the actual :math:U_T and compare against this reference.

Source code in src/groundfield/postprocess/safety.py

def permissible_touch_voltage_en50522(t_clear_s: float) -> float:
    """Permissible touch voltage :math:`U_{TP}(t_F)` per EN 50522:2010.

    Returns the maximum admissible touch voltage as a function of
    the fault-clearing time :math:`t_F`. The reference values are
    taken **verbatim** from EN 50522:2010, Table B.4 — the
    normative tabulation of the curve in Figure B.3 (values
    rounded to 5 V in the standard).

    Mathematically the helper is a piecewise-loglog interpolation
    over the canonical anchor points

    .. math::

        \\{(t_k, U_{TP,k})\\} \\;=\\;
        \\{(0.05, 725),\\ (0.10, 655),\\ (0.20, 525),\\
        (0.50, 225),\\ (1.00, 115),\\ (2.00, 95),\\
        (5.00, 85),\\ (10.0, 85)\\},

    with units :math:`\\mathrm{s}` and :math:`\\mathrm{V}`. Outside
    the grid the values are clamped to the endpoints — the standard
    does not extend the curve below 50 ms, and the terminal plateau
    at 85 V is explicit (the table reports identical values at 5 s
    and 10 s).

    Parameters
    ----------
    t_clear_s
        Fault-clearing time :math:`t_F` in seconds. Must be > 0.

    Returns
    -------
    float
        Permissible touch voltage :math:`U_{TP}` in V.

    Raises
    ------
    ValueError
        If ``t_clear_s`` is not strictly positive.

    Notes
    -----
    The EN 50522 curve is intended for engineering design; this
    helper is used as a *reference line* on plots of computed touch
    voltages, not as a regulatory pass/fail. Use
    :func:`touch_voltage` to compute the actual :math:`U_T` and
    compare against this reference.
    """
    if not math.isfinite(t_clear_s) or t_clear_s <= 0.0:
        raise ValueError(f"t_clear_s must be > 0, got {t_clear_s!r}.")

    ts = np.array([t for t, _ in _EN50522_TP_GRID], dtype=float)
    us = np.array([u for _, u in _EN50522_TP_GRID], dtype=float)
    if t_clear_s <= ts[0]:
        return float(us[0])
    if t_clear_s >= ts[-1]:
        return float(us[-1])

    # Log-log linear interpolation.
    log_t = math.log(t_clear_s)
    log_ts = np.log(ts)
    log_us = np.log(us)
    return float(math.exp(np.interp(log_t, log_ts, log_us)))

step_voltage

step_voltage(
    result: "FieldResult",
    *,
    position: tuple[float, float, float],
    direction: tuple[float, float, float] = (1.0, 0.0, 0.0),
    step: float = 1.0,
    surface_z: float | None = None,
    frequency_index: int = 0
) -> complex

Step voltage between two surface points.

Evaluates :math:U_S = \varphi(\mathbf{r}_1) - \varphi(\mathbf{r}_1 + d_\text{step}\,\hat{\mathbf{e}}) on the soil surface.

Parameters:

Name	Type	Description	Default
result	'FieldResult'	Solver output.	required
position	tuple[float, float, float]	First foot position :math:(x, y, z) in metres. The :math:z coordinate is overridden by surface_z if the latter is given; the default keeps z as supplied so callers can use the helper for buried-step studies if desired.	required
direction	tuple[float, float, float]	Step direction in world coordinates. The horizontal projection is used (the :math:z component is ignored — both feet stay on the same surface).	(1.0, 0.0, 0.0)
step	float	Step length :math:d_\text{step} in metres. EN 50522 uses 1.0.	1.0
surface_z	float | None	Optional explicit surface :math:z. None keeps the :math:z of position (typical for bare-foot contact at the ground surface, position = (x, y, 0.0)).	None
frequency_index	int	Index into :attr:FieldResult.frequencies.	0

Returns:

Type	Description
complex	Phasor :math:U_S(f) in V.

Raises:

Type	Description
ValueError	If step is not strictly positive or direction has no horizontal component.

Source code in src/groundfield/postprocess/safety.py

def step_voltage(
    result: "FieldResult",
    *,
    position: tuple[float, float, float],
    direction: tuple[float, float, float] = (1.0, 0.0, 0.0),
    step: float = 1.0,
    surface_z: float | None = None,
    frequency_index: int = 0,
) -> complex:
    """Step voltage between two surface points.

    Evaluates :math:`U_S = \\varphi(\\mathbf{r}_1) -
    \\varphi(\\mathbf{r}_1 + d_\\text{step}\\,\\hat{\\mathbf{e}})`
    on the soil surface.

    Parameters
    ----------
    result
        Solver output.
    position
        First foot position :math:`(x, y, z)` in metres. The
        :math:`z` coordinate is overridden by ``surface_z`` if the
        latter is given; the default keeps ``z`` as supplied so
        callers can use the helper for buried-step studies if
        desired.
    direction
        Step direction in world coordinates. The horizontal
        projection is used (the :math:`z` component is ignored —
        both feet stay on the same surface).
    step
        Step length :math:`d_\\text{step}` in metres. EN 50522
        uses ``1.0``.
    surface_z
        Optional explicit surface :math:`z`. ``None`` keeps the
        :math:`z` of ``position`` (typical for bare-foot contact
        at the ground surface, ``position = (x, y, 0.0)``).
    frequency_index
        Index into :attr:`FieldResult.frequencies`.

    Returns
    -------
    complex
        Phasor :math:`U_S(f)` in V.

    Raises
    ------
    ValueError
        If ``step`` is not strictly positive or ``direction`` has
        no horizontal component.
    """
    if not math.isfinite(step) or step <= 0.0:
        raise ValueError(f"step must be > 0, got {step!r}.")

    e_hat = _normalised(direction)
    e_h = np.array([e_hat[0], e_hat[1], 0.0])
    if float(np.linalg.norm(e_h)) <= 0.0:
        raise ValueError(
            f"direction has no horizontal component: {tuple(direction)}."
        )
    e_h = e_h / float(np.linalg.norm(e_h))

    z0 = float(position[2]) if surface_z is None else float(surface_z)
    p1 = np.array([float(position[0]), float(position[1]), z0])
    p2 = p1 + step * e_h
    pts = np.stack([p1, p2])

    phi = result.potential(pts, frequency_index=frequency_index)
    return complex(phi[0] - phi[1])

touch_voltage

touch_voltage(
    result: "FieldResult",
    world: "World",
    *,
    electrode: str,
    distance: float = 1.0,
    direction: tuple[float, float, float] = (1.0, 0.0, 0.0),
    surface_z: float = 0.0,
    frequency_index: int = 0
) -> complex

Touch voltage at a single point on the soil surface.

Evaluates :math:U_T = U_E - \varphi(\mathbf{r}_\text{feet}), where :math:U_E is the cluster potential of the touched electrode and :math:\mathbf{r}_\text{feet} is the surface point distance metres away from the electrode's :attr:connection_point in the (horizontal projection of the) given direction.

Parameters:

Name	Type	Description	Default
result	'FieldResult'	Solver output from :meth:Engine.solve.	required
world	'World'	Companion world used to resolve the electrode's connection point.	required
electrode	str	Name of the touched electrode. Its cluster potential is the EPR seen by the hand.	required
distance	float	Horizontal distance between the touched part and the feet, in metres. EN 50522 uses 1.0.	1.0
direction	tuple[float, float, float]	Direction in which the feet sit, expressed in world coordinates. The horizontal projection of this vector is used (the :math:z component is ignored on purpose — feet stay on the surface). Defaults to +x.	(1.0, 0.0, 0.0)
surface_z	float	:math:z coordinate of the soil surface in metres. Default 0.0 (groundfield convention; positive z is below ground).	0.0
frequency_index	int	Index into :attr:FieldResult.frequencies. Default 0.	0

Returns:

Type	Description
complex	Phasor :math:U_T(f) = U_E(f) - \varphi(\mathbf{r}_\text{feet}, f) in V.

Raises:

Type	Description
KeyError	If electrode is unknown to world or result.
ValueError	If distance is not strictly positive or direction is the zero vector after horizontal projection.

Source code in src/groundfield/postprocess/safety.py

def touch_voltage(
    result: "FieldResult",
    world: "World",
    *,
    electrode: str,
    distance: float = 1.0,
    direction: tuple[float, float, float] = (1.0, 0.0, 0.0),
    surface_z: float = 0.0,
    frequency_index: int = 0,
) -> complex:
    """Touch voltage at a single point on the soil surface.

    Evaluates :math:`U_T = U_E - \\varphi(\\mathbf{r}_\\text{feet})`,
    where :math:`U_E` is the cluster potential of the touched
    electrode and :math:`\\mathbf{r}_\\text{feet}` is the surface
    point ``distance`` metres away from the electrode's
    :attr:`connection_point` in the (horizontal projection of the)
    given ``direction``.

    Parameters
    ----------
    result
        Solver output from :meth:`Engine.solve`.
    world
        Companion world used to resolve the electrode's connection
        point.
    electrode
        Name of the touched electrode. Its cluster potential is the
        EPR seen by the hand.
    distance
        Horizontal distance between the touched part and the feet,
        in metres. EN 50522 uses ``1.0``.
    direction
        Direction in which the feet sit, expressed in world
        coordinates. The horizontal projection of this vector is
        used (the :math:`z` component is ignored on purpose — feet
        stay on the surface). Defaults to ``+x``.
    surface_z
        :math:`z` coordinate of the soil surface in metres.
        Default ``0.0`` (groundfield convention; positive ``z`` is
        below ground).
    frequency_index
        Index into :attr:`FieldResult.frequencies`. Default 0.

    Returns
    -------
    complex
        Phasor :math:`U_T(f) = U_E(f) - \\varphi(\\mathbf{r}_\\text{feet}, f)`
        in V.

    Raises
    ------
    KeyError
        If ``electrode`` is unknown to ``world`` or ``result``.
    ValueError
        If ``distance`` is not strictly positive or ``direction``
        is the zero vector after horizontal projection.
    """
    if not math.isfinite(distance) or distance <= 0.0:
        raise ValueError(f"distance must be > 0, got {distance!r}.")

    cp = world.get_electrode(electrode).connection_point
    e_hat = _normalised(direction)
    # Horizontal projection: enforce feet on the soil surface.
    e_h = np.array([e_hat[0], e_hat[1], 0.0])
    if float(np.linalg.norm(e_h)) <= 0.0:
        raise ValueError(
            f"direction has no horizontal component: {tuple(direction)}."
        )
    e_h = e_h / float(np.linalg.norm(e_h))

    feet = np.array([cp[0], cp[1], surface_z]) + distance * e_h
    phi_feet = complex(
        result.potential(feet[None, :], frequency_index=frequency_index)[0]
    )
    U_E = _cluster_potential(result, electrode, frequency_index)
    return U_E - phi_feet

touch_voltage_envelope

touch_voltage_envelope(
    result: "FieldResult",
    world: "World",
    *,
    electrode: str,
    distance: float = 1.0,
    n_angles: int = 24,
    surface_z: float = 0.0,
    frequency_index: int = 0
) -> tuple[np.ndarray, np.ndarray]

Touch-voltage profile around an electrode.

Walks a horizontal circle of radius distance around the electrode's :attr:connection_point in n_angles equal angular steps and returns the touch voltage at every direction. The maximum of |U_T| is the conservative envelope used in safety verification.

Parameters:

Name	Type	Description	Default
result	'FieldResult'	See :func:touch_voltage.	required
world	'FieldResult'	See :func:touch_voltage.	required
electrode	str	Name of the touched electrode.	required
distance	float	Radius of the circle in metres (1.0 by EN 50522).	1.0
n_angles	int	Number of equally spaced sample directions :math:\theta_k = 2\pi k / n with :math:k = 0, \dots, n - 1. Must be >= 3.	24
surface_z	float	See :func:touch_voltage.	0.0
frequency_index	float	See :func:touch_voltage.	0.0

Returns:

Name	Type	Description
angles	(ndarray, shape(n_angles))	Angles :math:\theta_k in radians, measured from the +x axis.
voltages	np.ndarray, shape (n_angles,), dtype complex	Touch-voltage phasor :math:U_T(\theta_k) in V.

Source code in src/groundfield/postprocess/safety.py

def touch_voltage_envelope(
    result: "FieldResult",
    world: "World",
    *,
    electrode: str,
    distance: float = 1.0,
    n_angles: int = 24,
    surface_z: float = 0.0,
    frequency_index: int = 0,
) -> tuple[np.ndarray, np.ndarray]:
    """Touch-voltage profile around an electrode.

    Walks a horizontal circle of radius ``distance`` around the
    electrode's :attr:`connection_point` in ``n_angles`` equal
    angular steps and returns the touch voltage at every direction.
    The maximum of ``|U_T|`` is the conservative envelope used in
    safety verification.

    Parameters
    ----------
    result, world
        See :func:`touch_voltage`.
    electrode
        Name of the touched electrode.
    distance
        Radius of the circle in metres (1.0 by EN 50522).
    n_angles
        Number of equally spaced sample directions
        :math:`\\theta_k = 2\\pi k / n` with :math:`k = 0,
        \\dots, n - 1`. Must be ``>= 3``.
    surface_z, frequency_index
        See :func:`touch_voltage`.

    Returns
    -------
    angles : np.ndarray, shape (n_angles,)
        Angles :math:`\\theta_k` in radians, measured from the
        ``+x`` axis.
    voltages : np.ndarray, shape (n_angles,), dtype complex
        Touch-voltage phasor :math:`U_T(\\theta_k)` in V.
    """
    if n_angles < 3:
        raise ValueError(f"n_angles must be >= 3, got {n_angles}.")
    if not math.isfinite(distance) or distance <= 0.0:
        raise ValueError(f"distance must be > 0, got {distance!r}.")

    cp = world.get_electrode(electrode).connection_point
    angles = np.linspace(0.0, 2.0 * np.pi, n_angles, endpoint=False)
    feet = np.column_stack(
        [
            cp[0] + distance * np.cos(angles),
            cp[1] + distance * np.sin(angles),
            np.full(n_angles, surface_z, dtype=float),
        ]
    )
    phi_feet = result.potential(feet, frequency_index=frequency_index)
    U_E = _cluster_potential(result, electrode, frequency_index)
    return angles, U_E - phi_feet

Current sharing and split factor

The :mod:groundfield.postprocess.current_balance module turns the per-electrode currents in :class:FieldResult into the engineering quantities that answer "where does the injected source current actually return?".

Physical background

For each galvanic cluster \(c\) the net soil leakage is

\[ I_c \;=\; \sum_{e \in c} I_e, \]

with \(I_e\) the per-electrode soil-leakage current stored in :attr:FieldResult.electrode_currents (positive in direction electrode → soil). Because the cluster members share a potential \(U_c\), the cluster impedance is \(Z_c = U_c / I_c\), and the per-electrode share inside one cluster is the complex ratio \(s_{e \mid c} = I_e / I_c\).

The split factor of a current source is

\[ s \;=\; \frac{\sum_{e \in c_\text{src}} I_e}{I_\text{src}}, \]

with \(c_\text{src}\) the cluster of the source's attached_to electrode. By construction, \(s = 1 + 0\,j\) when the entire injected current leaves the source cluster through the soil (no metallic parallel path), and \(s < 1\) when a metallic conductor (PEN trunk, parallel measurement lead, cable shield) carries part of the current as a parallel resistive path. This is the galvanic current division across parallel paths.

Not the same as the Reduktionsfaktor

In the German EVU / Schirmtechnik literature (Oeding & Oswald 2016) the Reduktionsfaktor refers to the additional transformatorische / inductive coupling correction between a current-carrying conductor and a parallel grounding / shield conductor. That quantity is angle-dependent: it vanishes when the two conductors are perpendicular (no flux linkage) but is large for collinear runs.

The split factor here is purely galvanic — present whenever there are parallel resistive paths, irrespective of the geometric angle between conductors. The proper Reduktionsfaktor is on the roadmap; the inductance backends in :mod:groundfield.coupling (Neumann, Carson, Sommerfeld) are already in place and will be picked up by a dedicated future helper.

Example

import groundfield as gf

soil = gf.HomogeneousSoil(resistivity=100.0)
world = gf.create_world(soil=soil)
g1 = gf.create_electrode(
    world, "rod", name="g1", position=(0.0, 0.0, 0.5), length=1.5,
)
g_aux = gf.create_electrode(
    world, "rod", name="g_aux", position=(50.0, 0.0, 0.5), length=1.5,
)
# Finite-impedance metallic feed lead between source and aux cluster.
gf.create_conductor(
    world, name="feed_lead", start=g1, end=g_aux,
    conductor_type="bare_copper", cross_section=50e-6,
)
gf.create_source(
    world, name="src", attached_to="g1", return_to="g_aux", magnitude=10.0,
)
result = gf.create_engine(backend="image", segment_length=0.05).solve(world)

# Per-cluster summary, sorted by descending |ΣI|:
gf.cluster_current_balance(result)
# Per-electrode table (with kind / depth annotations):
gf.electrode_current_table(result, world=world)
# Split factor of the source — < 1 because the metallic feed lead
# carries part of I_src as a parallel resistive path:
s = gf.split_factor(result, world)
# Top-15 bar chart of |I| per electrode:
gf.plot_current_sharing(result, world=world, top_n=15)

Validity envelope

• Frequency: quasi-static envelope \(f \le 1\,\mathrm{kHz}\).
• Backends: every solver that populates electrode_potentials / electrode_currents / clusters (image, image_2layer, image_nlayer, mom, mom_sommerfeld, cim, bem). Stub backends produce empty / NaN cluster tables.
• The image backend treats Source.return_to as informational metadata; the injected current dissipates via the Dirichlet far-field boundary. The split factor still detects metallic parallel paths because the cluster's net soil leakage is reduced by whatever current is shunted into a finite-impedance branch.

API reference — current sharing

current_balance

Current sharing and split-factor analysis for a FieldResult.

This module turns the per-electrode currents stored in :class:groundfield.solver.result.FieldResult into the engineering quantities that answer "where does the injected source current actually return?". In a TN-Ortsnetz with hundreds of houses, dozens of cable cabinets and one or more metallic return paths (PEN trunk, measurement leads, cable shields), the soil-leakage current of every cluster matters and cannot be read off electrode_currents directly without aggregation.

Quantities

• Per-cluster soil leakage :math:I_{c} = \sum_{e \in c} I_e. With ideal galvanic bonds (cross_section=None) the cluster members share a potential :math:U_c; the cluster impedance is then :math:Z_c = U_c / I_c. Returned as a tabular summary by :func:cluster_current_balance.
• Per-electrode share of cluster current — what fraction of the cluster's net soil leakage flows through this specific electrode: :math:s_{e \mid c} = I_e / I_c (complex). Returned alongside the per-electrode potential / impedance by :func:electrode_current_table.
• Split factor :math:s = I_{c_\text{src}} / I_\text{src} = \sum_{e \in c_\text{src}} I_e / I_\text{src}, with :math:c_\text{src} the cluster of the source's attached_to electrode. s = 1 means the injected current leaves the source cluster entirely through the soil (no metallic short-cut to the return-path cluster). s < 1 means a metallic conductor (PEN trunk, parallel measurement lead, cable shield) carries part of the current as a parallel resistive path.

This is the galvanic current-split between parallel paths. It is not the Reduktionsfaktor in the German EVU / Schirmtechnik sense (Oeding & Oswald 2016): that latter quantity is the additional transformatorische / inductive coupling correction between a current-carrying conductor and a parallel grounding / shield conductor. The Reduktionsfaktor vanishes when the two conductors are perpendicular (no flux linkage) but the split factor still applies — the current always splits among parallel resistive paths irrespective of the geometric angle. A future helper may add the proper Reduktionsfaktor based on the Neumann mutual-inductance / Carson / Sommerfeld backends already present in :mod:groundfield.coupling.

Validity envelope

• Frequency: quasi-static envelope :math:f \le 1\,\mathrm{kHz}.
• Conventions: :class:FieldResult.electrode_currents carries the per-electrode soil-leakage current in A (complex phasor), with positive sign in the direction electrode → soil. All quantities here are complex per frequency_index.
• Cluster discovery uses :attr:FieldResult.clusters — the same map produced by every backend.

cluster_current_balance

cluster_current_balance(
    result: "FieldResult", *, frequency_index: int = 0
) -> pd.DataFrame

Per-cluster soil leakage, potential and impedance.

Returns one row per unique galvanic cluster present in the result. The cluster impedance is computed as :math:Z_c = U_c / I_c with :math:I_c = \sum_{e \in c} I_e. Where :math:I_c = 0 (purely passive observer cluster) the impedance columns are NaN.

Parameters:

Name	Type	Description	Default
result	'FieldResult'	Solver output.	required
frequency_index	int	Index into :attr:FieldResult.frequencies. Default 0.	0

Returns:

Type	Description
DataFrame	Columns cluster_root, n_members, members (list of strings), U_re, U_im, abs_U, sum_I_re, sum_I_im, abs_sum_I, Z_re, Z_im, abs_Z, arg_Z_deg. Rows are sorted by descending abs_sum_I so the dominant clusters surface first — the typical debugging order.

Source code in src/groundfield/postprocess/current_balance.py

def cluster_current_balance(
    result: "FieldResult",
    *,
    frequency_index: int = 0,
) -> pd.DataFrame:
    """Per-cluster soil leakage, potential and impedance.

    Returns one row per *unique* galvanic cluster present in the
    result. The cluster impedance is computed as
    :math:`Z_c = U_c / I_c` with :math:`I_c = \\sum_{e \\in c} I_e`.
    Where :math:`I_c = 0` (purely passive observer cluster) the
    impedance columns are ``NaN``.

    Parameters
    ----------
    result
        Solver output.
    frequency_index
        Index into :attr:`FieldResult.frequencies`. Default 0.

    Returns
    -------
    pandas.DataFrame
        Columns
        ``cluster_root``,
        ``n_members``,
        ``members`` (list of strings),
        ``U_re``, ``U_im``, ``abs_U``,
        ``sum_I_re``, ``sum_I_im``, ``abs_sum_I``,
        ``Z_re``, ``Z_im``, ``abs_Z``, ``arg_Z_deg``.
        Rows are sorted by descending ``abs_sum_I`` so the
        dominant clusters surface first — the typical
        debugging order.
    """
    rows = []
    for root, members in _unique_clusters(result).items():
        # Cluster potential = potential of any member (ideal cluster
        # bond -> shared phi). Use the canonical root.
        U = complex(result.electrode_potentials[root][frequency_index])
        sum_I = 0.0 + 0j
        for m in members:
            if m not in result.electrode_currents:
                continue
            sum_I += complex(result.electrode_currents[m][frequency_index])
        if sum_I != 0:
            Z = U / sum_I
        else:
            Z = complex("nan")
        rows.append(
            {
                "cluster_root": root,
                "n_members": len(members),
                "members": list(members),
                "U_re": U.real,
                "U_im": U.imag,
                "abs_U": abs(U),
                "sum_I_re": sum_I.real,
                "sum_I_im": sum_I.imag,
                "abs_sum_I": abs(sum_I),
                "Z_re": Z.real if np.isfinite(Z.real) else float("nan"),
                "Z_im": Z.imag if np.isfinite(Z.imag) else float("nan"),
                "abs_Z": (
                    abs(Z) if (np.isfinite(Z.real) and np.isfinite(Z.imag)) else float("nan")
                ),
                "arg_Z_deg": (
                    math.degrees(cmath.phase(Z))
                    if (np.isfinite(Z.real) and np.isfinite(Z.imag))
                    else float("nan")
                ),
            }
        )
    df = pd.DataFrame(rows)
    if not df.empty:
        df = df.sort_values("abs_sum_I", ascending=False).reset_index(drop=True)
    return df

electrode_current_table

electrode_current_table(
    result: "FieldResult",
    world: "World | None" = None,
    *,
    frequency_index: int = 0
) -> pd.DataFrame

Per-electrode potential, current and share of cluster current.

For each electrode in the result the table reports the potential :math:U_e (= cluster potential), the soil-leakage current :math:I_e, the implied two-terminal impedance :math:Z_e = U_e / I_e, and the fractional share of the cluster total :math:s_{e \mid c} = I_e / I_c. The latter is the engineering key indicator for "of all the soil current leaving cluster c, how much physically flows through this specific electrode" — useful when looking at a single transformer station with a ring + several rods, or a building cluster with multiple foundation parts.

Parameters:

Name	Type	Description	Default
result	'FieldResult'	Solver output.	required
world	'World | None'	Optional companion world. When given, the table includes the electrode kind (rod / ring / mesh / ...) and the connection-point depth in metres — a small but very helpful annotation when scanning a 200-EFH typical run.	None
frequency_index	int	Index into :attr:FieldResult.frequencies. Default 0.	0

Returns:

Type	Description
DataFrame	Columns name, cluster_root, kind (if world), depth_m (if world), U_re, U_im, abs_U, I_re, I_im, abs_I, Z_re, Z_im, abs_Z, arg_Z_deg, share_of_cluster_re, share_of_cluster_im, share_of_cluster_pct. Sorted by descending abs_I.

Source code in src/groundfield/postprocess/current_balance.py

def electrode_current_table(
    result: "FieldResult",
    world: "World | None" = None,
    *,
    frequency_index: int = 0,
) -> pd.DataFrame:
    """Per-electrode potential, current and share of cluster current.

    For each electrode in the result the table reports the
    potential :math:`U_e` (= cluster potential), the soil-leakage
    current :math:`I_e`, the implied two-terminal impedance
    :math:`Z_e = U_e / I_e`, and the **fractional share** of the
    cluster total :math:`s_{e \\mid c} = I_e / I_c`. The latter is
    the engineering key indicator for "of all the soil current
    leaving cluster *c*, how much physically flows through this
    specific electrode" — useful when looking at a single
    transformer station with a ring + several rods, or a building
    cluster with multiple foundation parts.

    Parameters
    ----------
    result
        Solver output.
    world
        Optional companion world. When given, the table includes
        the electrode kind (``rod`` / ``ring`` / ``mesh`` / ...)
        and the connection-point depth in metres — a small but
        very helpful annotation when scanning a 200-EFH typical run.
    frequency_index
        Index into :attr:`FieldResult.frequencies`. Default 0.

    Returns
    -------
    pandas.DataFrame
        Columns
        ``name``, ``cluster_root``, ``kind`` (if ``world``),
        ``depth_m`` (if ``world``),
        ``U_re``, ``U_im``, ``abs_U``,
        ``I_re``, ``I_im``, ``abs_I``,
        ``Z_re``, ``Z_im``, ``abs_Z``, ``arg_Z_deg``,
        ``share_of_cluster_re``, ``share_of_cluster_im``,
        ``share_of_cluster_pct``. Sorted by descending ``abs_I``.
    """
    # Map every electrode to its cluster root and to the cluster's
    # net soil-leakage current at the chosen frequency.
    cluster_root_of: dict[str, str] = {}
    cluster_sum_I: dict[str, complex] = {}
    for root, members in _unique_clusters(result).items():
        s = 0.0 + 0j
        for m in members:
            if m in result.electrode_currents:
                s += complex(result.electrode_currents[m][frequency_index])
        cluster_sum_I[root] = s
        for m in members:
            cluster_root_of[m] = root

    geo_lookup: dict[str, tuple[str, float]] = {}
    if world is not None:
        for e in world.electrodes:
            cp = e.connection_point
            geo_lookup[e.name] = (e.kind, float(cp[2]))

    rows = []
    for name in sorted(result.electrode_potentials.keys()):
        U = complex(result.electrode_potentials[name][frequency_index])
        I = complex(
            result.electrode_currents.get(name, [0.0 + 0j] * (frequency_index + 1))[
                frequency_index
            ]
        )
        if I != 0:
            Z = U / I
        else:
            Z = complex("nan")
        root = cluster_root_of.get(name, name)
        I_c = cluster_sum_I.get(root, 0.0 + 0j)
        if I_c != 0:
            s_e = I / I_c
            s_pct = abs(s_e) * 100.0
        else:
            s_e = complex("nan")
            s_pct = float("nan")

        row: dict[str, object] = {
            "name": name,
            "cluster_root": root,
        }
        if world is not None:
            kind, depth = geo_lookup.get(name, ("?", float("nan")))
            row["kind"] = kind
            row["depth_m"] = depth
        row.update(
            {
                "U_re": U.real,
                "U_im": U.imag,
                "abs_U": abs(U),
                "I_re": I.real,
                "I_im": I.imag,
                "abs_I": abs(I),
                "Z_re": Z.real if np.isfinite(Z.real) else float("nan"),
                "Z_im": Z.imag if np.isfinite(Z.imag) else float("nan"),
                "abs_Z": (
                    abs(Z) if (np.isfinite(Z.real) and np.isfinite(Z.imag)) else float("nan")
                ),
                "arg_Z_deg": (
                    math.degrees(cmath.phase(Z))
                    if (np.isfinite(Z.real) and np.isfinite(Z.imag))
                    else float("nan")
                ),
                "share_of_cluster_re": (
                    s_e.real if np.isfinite(s_e.real) else float("nan")
                ),
                "share_of_cluster_im": (
                    s_e.imag if np.isfinite(s_e.imag) else float("nan")
                ),
                "share_of_cluster_pct": s_pct,
            }
        )
        rows.append(row)
    df = pd.DataFrame(rows)
    if not df.empty:
        df = df.sort_values("abs_I", ascending=False).reset_index(drop=True)
    return df

plot_current_sharing

plot_current_sharing(
    result: "FieldResult",
    world: "World | None" = None,
    *,
    by: Literal["electrode", "cluster"] = "electrode",
    top_n: int = 15,
    frequency_index: int = 0,
    figsize: tuple[float, float] = (8.0, 5.0),
    title: str | None = None
)

Bar chart of the top_n current contributors.

Renders |I| (in A) for either every electrode or every galvanic cluster, sorted descending. The default by = "electrode" is the default — which physical electrode actually carries the test current? The companion option by = "cluster" aggregates over each cluster.

Parameters:

Name	Type	Description	Default
result	'FieldResult'	Solver output.	required
world	'World | None'	Optional companion world. Forwarded to :func:electrode_current_table for the kind / depth annotations; ignored otherwise.	None
by	Literal['electrode', 'cluster']	"electrode" (default) or "cluster".	'electrode'
top_n	int	Maximum number of bars to render. Smaller contributors are truncated. Pass 0 to render all.	15
frequency_index	int	Index into :attr:FieldResult.frequencies.	0
figsize	tuple[float, float]	Matplotlib figure size in inches.	(8.0, 5.0)
title	str | None	Optional title override.	None

Returns:

Type	Description
Figure

Source code in src/groundfield/postprocess/current_balance.py

def plot_current_sharing(
    result: "FieldResult",
    world: "World | None" = None,
    *,
    by: Literal["electrode", "cluster"] = "electrode",
    top_n: int = 15,
    frequency_index: int = 0,
    figsize: tuple[float, float] = (8.0, 5.0),
    title: str | None = None,
):
    """Bar chart of the ``top_n`` current contributors.

    Renders ``|I|`` (in A) for either every electrode or every
    galvanic cluster, sorted descending. The default ``by =
    "electrode"`` is the default — *which physical electrode
    actually carries the test current?* The companion option
    ``by = "cluster"`` aggregates over each cluster.

    Parameters
    ----------
    result
        Solver output.
    world
        Optional companion world. Forwarded to
        :func:`electrode_current_table` for the kind / depth
        annotations; ignored otherwise.
    by
        ``"electrode"`` (default) or ``"cluster"``.
    top_n
        Maximum number of bars to render. Smaller contributors are
        truncated. Pass ``0`` to render all.
    frequency_index
        Index into :attr:`FieldResult.frequencies`.
    figsize
        Matplotlib figure size in inches.
    title
        Optional title override.

    Returns
    -------
    matplotlib.figure.Figure
    """
    import matplotlib.pyplot as plt

    if by == "electrode":
        df = electrode_current_table(
            result, world=world, frequency_index=frequency_index
        )
        labels = df["name"].tolist()
        magnitudes = df["abs_I"].to_numpy()
        ylabel = "|I| in A (per electrode)"
    elif by == "cluster":
        df = cluster_current_balance(result, frequency_index=frequency_index)
        labels = df["cluster_root"].tolist()
        magnitudes = df["abs_sum_I"].to_numpy()
        ylabel = "|ΣI| in A (per cluster)"
    else:
        raise ValueError(f"by must be 'electrode' or 'cluster', got {by!r}.")

    if top_n and top_n > 0:
        labels = labels[:top_n]
        magnitudes = magnitudes[:top_n]

    fig, ax = plt.subplots(figsize=figsize)
    ax.bar(range(len(labels)), magnitudes, color="C0")
    ax.set_xticks(range(len(labels)))
    ax.set_xticklabels(labels, rotation=45, ha="right")
    ax.set_ylabel(ylabel)
    ax.grid(True, axis="y", alpha=0.3)

    if title is None:
        f_hz = result.frequencies[frequency_index]
        suffix = f"top {top_n}" if (top_n and top_n > 0) else "all"
        title = f"Current sharing by {by} — {suffix}, f = {f_hz:g} Hz"
    ax.set_title(title)
    fig.tight_layout()
    return fig

split_factor

split_factor(
    result: "FieldResult",
    world: "World",
    *,
    source_name: str | None = None,
    frequency_index: int = 0
) -> complex

Split factor :math:s = I_{c_\text{src}} / I_\text{src}.

Computes the galvanic split factor of a current source: the fraction of the injected current that leaves the source cluster through the soil rather than through any metallic parallel path (PEN trunk, parallel measurement lead, cable shield).

Mathematically

.. math::

s \;=\; \frac{\sum_{e \in c_\text{src}} I_e}{I_\text{src}},

with :math:c_\text{src} the galvanic cluster of the source's attached_to electrode and :math:I_\text{src} the complex phasor of the source. By construction:

• For a stand-alone source on an electrode without metallic parallel paths, KCL forces :math:s = 1 + 0\,j. The injected current leaves the cluster entirely through the soil.
• For a source whose cluster is connected to another grounding cluster (e.g. the auxiliary electrode of a fall-of-potential measurement) by a metallic conductor, :math:s < 1. The smaller the conductor's series impedance, the smaller :math:s.
• The complex argument :math:\arg s reveals the inductive content of the parallel path — a Carson- or Sommerfeld- corrected lead at frequency >0 produces an imaginary share.

Not to be confused with the Reduktionsfaktor

In the German EVU / Schirmtechnik literature (Oeding & Oswald 2016), the Reduktionsfaktor refers to the transformatorische / inductive coupling correction between a current-carrying conductor and a parallel grounding / shield conductor. That quantity is angle-dependent: it vanishes when the two conductors are perpendicular (no flux linkage) but is large for collinear runs.

The split factor implemented here is purely galvanic — the resistive division of current across parallel paths. It is present whenever there are multiple parallel paths, regardless of their geometric orientation.

Parameters:

Name	Type	Description	Default
result	'FieldResult'	Solver output.	required
world	'World'	Companion world used to resolve the source object.	required
source_name	str | None	Optional name of the current source. None (default) picks the unique current source; raises if the world contains more than one.	None
frequency_index	int	Index into :attr:FieldResult.frequencies. Default 0.	0

Returns:

Type	Description
complex	Phasor :math:s(f) (dimensionless).

Raises:

Type	Description
ValueError	If the world has no current source, multiple current sources without explicit source_name, or the source has zero magnitude.
KeyError	If source_name does not exist among current sources, or the source's attached_to electrode is unknown to the result.

Source code in src/groundfield/postprocess/current_balance.py

def split_factor(
    result: "FieldResult",
    world: "World",
    *,
    source_name: str | None = None,
    frequency_index: int = 0,
) -> complex:
    """Split factor :math:`s = I_{c_\\text{src}} / I_\\text{src}`.

    Computes the **galvanic split factor** of a current source:
    the fraction of the injected current that leaves the source
    cluster through the **soil** rather than through any metallic
    parallel path (PEN trunk, parallel measurement lead, cable
    shield).

    Mathematically

    .. math::

        s \\;=\\; \\frac{\\sum_{e \\in c_\\text{src}} I_e}{I_\\text{src}},

    with :math:`c_\\text{src}` the galvanic cluster of the source's
    ``attached_to`` electrode and :math:`I_\\text{src}` the complex
    phasor of the source. By construction:

    * For a stand-alone source on an electrode without metallic
      parallel paths, KCL forces :math:`s = 1 + 0\\,j`. The
      injected current leaves the cluster entirely through the
      soil.
    * For a source whose cluster is connected to *another*
      grounding cluster (e.g. the auxiliary electrode of a
      fall-of-potential measurement) by a metallic conductor,
      :math:`s < 1`. The smaller the conductor's series
      impedance, the smaller :math:`s`.
    * The complex argument :math:`\\arg s` reveals the inductive
      content of the parallel path — a Carson- or Sommerfeld-
      corrected lead at frequency >0 produces an imaginary share.

    Not to be confused with the *Reduktionsfaktor*
    -----------------------------------------------
    In the German EVU / Schirmtechnik literature (Oeding & Oswald
    2016), the *Reduktionsfaktor* refers to the **transformatorische
    / inductive coupling correction** between a current-carrying
    conductor and a parallel grounding / shield conductor. That
    quantity is angle-dependent: it vanishes when the two conductors
    are perpendicular (no flux linkage) but is large for collinear
    runs.

    The split factor implemented here is **purely galvanic** — the
    resistive division of current across parallel paths. It is
    present whenever there are multiple parallel paths, regardless
    of their geometric orientation.

    Parameters
    ----------
    result
        Solver output.
    world
        Companion world used to resolve the source object.
    source_name
        Optional name of the current source. ``None`` (default)
        picks the unique current source; raises if the world
        contains more than one.
    frequency_index
        Index into :attr:`FieldResult.frequencies`. Default 0.

    Returns
    -------
    complex
        Phasor :math:`s(f)` (dimensionless).

    Raises
    ------
    ValueError
        If the world has no current source, multiple current
        sources without explicit ``source_name``, or the source
        has zero magnitude.
    KeyError
        If ``source_name`` does not exist among current sources,
        or the source's ``attached_to`` electrode is unknown to
        the result.
    """
    source = _resolve_source(world, source_name)
    I_src = _source_current(source)
    if abs(I_src) == 0.0:
        raise ValueError(
            f"Source '{source.name}' has zero magnitude — split factor "
            "is undefined."
        )

    attached = source.attached_to
    members = result.clusters.get(attached, [attached])
    if not members:
        raise KeyError(
            f"Source '{source.name}' is attached to '{attached}', "
            "which has no cluster entry in the result."
        )

    sum_I = 0.0 + 0j
    for m in members:
        if m in result.electrode_currents:
            sum_I += complex(result.electrode_currents[m][frequency_index])
    return sum_I / I_src

Parameter sweeps and convergence studies

The :mod:groundfield.postprocess.sweep and :mod:groundfield.postprocess.convergence modules turn the typical parameter axes into a single tabular response. Both produce long-format :class:pandas.DataFrame objects that feed naturally into the :math:\rho-:math:f fit and the vector-fitting pipelines.

sweep — Cartesian product over named axes

import groundfield as gf

def world_factory(*, rho, h_1):
    soil = gf.TwoLayerSoil(rho_1=rho, rho_2=10*rho, h_1=h_1)
    world = gf.create_world(soil=soil)
    gf.create_electrode(world, "rod", name="g1",
                        position=(0, 0, 0.5), length=1.5)
    gf.create_source(world, attached_to="g1", magnitude=1.0)
    return world

eng = gf.create_engine(backend="image_2layer", segment_length=0.1,
                       frequencies=[50.0, 200.0])
df = gf.sweep(
    world_factory,
    eng,
    axes={"rho": [50, 100, 500, 1000], "h_1": [1.0, 2.0, 5.0]},
)
gf.plot_sweep_lines(df, x="rho", y="abs_Z", color="h_1",
                    log_x=True, log_y=True)
gf.plot_sweep_heatmap(df, x="rho", y="h_1", response="abs_Z",
                      frequency_Hz=50.0)

The default response captures the cluster impedance and EPR at the source's cluster per frequency (Z_re, Z_im, abs_Z, arg_Z_deg, U_E_re, U_E_im, abs_U_E, I_re, I_im, abs_I). Pass response=... to extract custom scalars.

convergence_study — refinement over segment_length

df = gf.convergence_study(
    world, eng, segment_lengths=[0.5, 0.2, 0.1, 0.05, 0.02],
)
gf.plot_convergence(df, response="abs_Z", reference=R_sunde)

The helper clones the engine via :meth:Engine.model_copy, so the original engine is not mutated. The plot's x-axis is inverted so finer segment_length sits on the right (the convergence-asymptote direction). Pass an analytical reference (Sunde, Dwight, IEEE 80) via reference=... for a clean asymptote line.

API reference — sweep

sweep

Cartesian-product parameter sweeps for parameter axis studies.

This module turns the typical study axes — soil resistivity :math:\rho_1 / :math:\rho_2, layer thickness :math:h_1, electrode geometry, frequency — into a single tabular response that is cheap to plot and easy to feed into vector-fitting, :math:\rho-:math:f regression, or downstream groundinsight.

API surface

:func:sweep Walk the Cartesian product of an arbitrary number of named axes, build a fresh :class:World (and optionally a fresh :class:Engine) per combination, solve, and extract a scalar response per frequency. Returns a long-format :class:pandas.DataFrame with one row per (axis values × frequency). :func:plot_sweep_lines Line plot of one response column against a chosen axis, optionally with one curve per value of a second axis. :func:plot_sweep_heatmap Pivot-table heatmap of one response column over a (x_axis, y_axis) pair (e.g. :math:\rho_1 vs. :math:h_1).

Mathematical / physical content

The default response extractor picks up the cluster impedance :math:Z_c = U_c / I_c at the source's galvanic cluster, with :math:U_c the cluster potential and :math:I_c = \sum_{e \in c_\text{src}} I_e the net soil leakage. Both numerator and denominator are read from :class:FieldResult per frequency_index, so the response is the engineering :math:Z(\rho_1, h_1, f) curve that :mod:groundfield.postprocess.rho_f_standard and :mod:groundfield.postprocess.vector_fitting consume.

Validity envelope

• Frequency: quasi-static envelope :math:f \le 1\,\mathrm{kHz}.
• Linearity: the sweep does not interpolate between samples — for a smooth :math:Z(\rho, f) surface, hand the resulting DataFrame to :func:fit_rho_f_standard / :func:fit_to_sympy_standard.
• Cost: every Cartesian combination triggers one full :meth:Engine.solve. For 50 ρ-values × 10 h-values × 1 frequency list with N segments each, expect 500 dense system solves; budget accordingly with :func:expected_segments.

plot_sweep_heatmap

plot_sweep_heatmap(
    df: pd.DataFrame,
    *,
    x: str,
    y: str,
    response: str = "abs_Z",
    frequency_Hz: float | None = None,
    agg: str = "mean",
    cmap: str = "viridis",
    figsize: tuple[float, float] = (7.5, 5.5),
    title: str | None = None
)

Heatmap of response over the (x, y) axis pair.

Pivot-tables df by y (rows) and x (columns), aggregating response with the chosen agg (default "mean" — useful when the sweep contains additional axes that should collapse). Selects a single frequency slice when frequency_Hz is given.

Parameters:

Name	Type	Description	Default
df	DataFrame	Long-format DataFrame from :func:sweep.	required
x	str	Column names for the heatmap axes. Both must be numeric.	required
y	str	Column names for the heatmap axes. Both must be numeric.	required
response	str	Column to colour-code. Default "abs_Z".	'abs_Z'
frequency_Hz	float | None	If given and the DataFrame has a frequency_Hz column, keep only rows that match this frequency. Raises if no rows match.	None
agg	str	Aggregation passed to :meth:pd.DataFrame.pivot_table. Default "mean".	'mean'
cmap	str	Standard matplotlib options.	'viridis'
figsize	str	Standard matplotlib options.	'viridis'
title	str	Standard matplotlib options.	'viridis'

Returns:

Type	Description
Figure

Source code in src/groundfield/postprocess/sweep.py

def plot_sweep_heatmap(
    df: pd.DataFrame,
    *,
    x: str,
    y: str,
    response: str = "abs_Z",
    frequency_Hz: float | None = None,
    agg: str = "mean",
    cmap: str = "viridis",
    figsize: tuple[float, float] = (7.5, 5.5),
    title: str | None = None,
):
    """Heatmap of ``response`` over the ``(x, y)`` axis pair.

    Pivot-tables ``df`` by ``y`` (rows) and ``x`` (columns),
    aggregating ``response`` with the chosen ``agg`` (default
    ``"mean"`` — useful when the sweep contains additional axes
    that should collapse). Selects a single frequency slice
    when ``frequency_Hz`` is given.

    Parameters
    ----------
    df
        Long-format DataFrame from :func:`sweep`.
    x, y
        Column names for the heatmap axes. Both must be numeric.
    response
        Column to colour-code. Default ``"abs_Z"``.
    frequency_Hz
        If given and the DataFrame has a ``frequency_Hz`` column,
        keep only rows that match this frequency. Raises if no
        rows match.
    agg
        Aggregation passed to :meth:`pd.DataFrame.pivot_table`.
        Default ``"mean"``.
    cmap, figsize, title
        Standard matplotlib options.

    Returns
    -------
    matplotlib.figure.Figure
    """
    import matplotlib.pyplot as plt

    for col in (x, y, response):
        if col not in df.columns:
            raise KeyError(
                f"column '{col}' not in DataFrame; have: {list(df.columns)}"
            )

    sub = df
    if frequency_Hz is not None and "frequency_Hz" in df.columns:
        sub = df[df["frequency_Hz"] == frequency_Hz]
        if sub.empty:
            raise ValueError(
                f"No rows match frequency_Hz={frequency_Hz}; "
                f"available: {sorted(df['frequency_Hz'].unique())}"
            )

    pivot = sub.pivot_table(index=y, columns=x, values=response, aggfunc=agg)
    fig, ax = plt.subplots(figsize=figsize)
    im = ax.imshow(
        pivot.values,
        origin="lower",
        aspect="auto",
        extent=(
            float(pivot.columns.min()), float(pivot.columns.max()),
            float(pivot.index.min()), float(pivot.index.max()),
        ),
        cmap=cmap,
        interpolation="nearest",
    )
    fig.colorbar(im, ax=ax, label=response)
    ax.set_xlabel(x)
    ax.set_ylabel(y)
    if title is None:
        title = f"{response}({x}, {y})"
        if frequency_Hz is not None:
            title += f"  @ f = {frequency_Hz:g} Hz"
    ax.set_title(title)
    fig.tight_layout()
    return fig

plot_sweep_lines

plot_sweep_lines(
    df: pd.DataFrame,
    *,
    x: str,
    y: str = "abs_Z",
    color: str | None = None,
    figsize: tuple[float, float] = (8.0, 5.0),
    log_x: bool = False,
    log_y: bool = False,
    title: str | None = None
)

Line plot of y versus x, one curve per color value.

Parameters:

Name	Type	Description	Default
df	DataFrame	Long-format DataFrame as produced by :func:sweep.	required
x	str	Column name to use on the x-axis.	required
y	str	Column name of the response. Default "abs_Z".	'abs_Z'
color	str | None	Optional second column name; one line per distinct value is drawn. None (default) collapses to a single line.	None
figsize	tuple[float, float]	Standard matplotlib options.	(8.0, 5.0)
log_x	tuple[float, float]	Standard matplotlib options.	(8.0, 5.0)
log_y	tuple[float, float]	Standard matplotlib options.	(8.0, 5.0)
title	tuple[float, float]	Standard matplotlib options.	(8.0, 5.0)

Returns:

Type	Description
Figure

Source code in src/groundfield/postprocess/sweep.py

def plot_sweep_lines(
    df: pd.DataFrame,
    *,
    x: str,
    y: str = "abs_Z",
    color: str | None = None,
    figsize: tuple[float, float] = (8.0, 5.0),
    log_x: bool = False,
    log_y: bool = False,
    title: str | None = None,
):
    """Line plot of ``y`` versus ``x``, one curve per ``color`` value.

    Parameters
    ----------
    df
        Long-format DataFrame as produced by :func:`sweep`.
    x
        Column name to use on the x-axis.
    y
        Column name of the response. Default ``"abs_Z"``.
    color
        Optional second column name; one line per distinct value
        is drawn. ``None`` (default) collapses to a single line.
    figsize, log_x, log_y, title
        Standard matplotlib options.

    Returns
    -------
    matplotlib.figure.Figure
    """
    import matplotlib.pyplot as plt

    if x not in df.columns:
        raise KeyError(f"column '{x}' not in DataFrame; have: {list(df.columns)}")
    if y not in df.columns:
        raise KeyError(f"column '{y}' not in DataFrame; have: {list(df.columns)}")

    fig, ax = plt.subplots(figsize=figsize)
    if color is None:
        sub = df.sort_values(x)
        ax.plot(sub[x], sub[y], marker="o")
    else:
        if color not in df.columns:
            raise KeyError(
                f"color column '{color}' not in DataFrame; have: {list(df.columns)}"
            )
        for val, sub in df.groupby(color):
            sub = sub.sort_values(x)
            ax.plot(sub[x], sub[y], marker="o", label=f"{color} = {val}")
        ax.legend()

    if log_x:
        ax.set_xscale("log")
    if log_y:
        ax.set_yscale("log")
    ax.set_xlabel(x)
    ax.set_ylabel(y)
    ax.grid(True, which="both", alpha=0.3)
    if title is None:
        title = f"{y} vs. {x}" + (f"  (color: {color})" if color else "")
    ax.set_title(title)
    fig.tight_layout()
    return fig

sweep

sweep(
    world_factory: Callable[..., "World"],
    engine: "Engine | Callable[..., Engine]",
    *,
    axes: dict[str, Sequence[Any]],
    response: (
        Callable[
            ["FieldResult", "World", int], dict[str, float]
        ]
        | None
    ) = None
) -> pd.DataFrame

Cartesian-product parameter sweep across user-defined axes.

For every combination (a_1 = v_1, a_2 = v_2, ...) in the Cartesian product of axes, the function

1. builds a fresh world via world_factory(**combination),
2. resolves the engine — either the static :class:Engine passed in or a per-combination engine via engine(**combination) if a callable is given,
3. solves with :meth:World.solve,
4. iterates over every frequency in :attr:FieldResult.frequencies and extracts a row via response(result, world, frequency_index).

The axis values and frequency_Hz are added to every row automatically.

Parameters:

Name	Type	Description	Default
world_factory	Callable[..., 'World']	Callable that builds a fresh :class:World from the per-combination kwargs. Must return a fully populated world (soil, electrodes, sources, optional conductors).	required
engine	'Engine | Callable[..., Engine]'	Either a static :class:Engine (reused for every combination) or a callable engine(**combination) -> Engine (rebuilt per combination — useful when the segment length depends on the geometry).	required
axes	dict[str, Sequence[Any]]	Mapping axis_name -> Sequence of values. Must contain at least one axis. Empty sequences raise immediately.	required
response	Callable[['FieldResult', 'World', int], dict[str, float]] | None	Optional response extractor. Receives (result, world, frequency_index) and must return a dict[str, float|complex]. Defaults to :func:_default_response (cluster impedance + EPR at the source's cluster).	None

Returns:

Type	Description
DataFrame	Long-format. Columns are the axis names, frequency_Hz, and every key returned by response. One row per Cartesian-product point per frequency.

Raises:

Type	Description
ValueError	If axes is empty or any axis is empty.

Source code in src/groundfield/postprocess/sweep.py

def sweep(
    world_factory: Callable[..., "World"],
    engine: "Engine | Callable[..., Engine]",
    *,
    axes: dict[str, Sequence[Any]],
    response: Callable[["FieldResult", "World", int], dict[str, float]] | None = None,
) -> pd.DataFrame:
    """Cartesian-product parameter sweep across user-defined axes.

    For every combination ``(a_1 = v_1, a_2 = v_2, ...)`` in the
    Cartesian product of ``axes``, the function

    1. builds a fresh world via ``world_factory(**combination)``,
    2. resolves the engine — either the static :class:`Engine`
       passed in or a per-combination engine via
       ``engine(**combination)`` if a callable is given,
    3. solves with :meth:`World.solve`,
    4. iterates over every frequency in
       :attr:`FieldResult.frequencies` and extracts a row via
       ``response(result, world, frequency_index)``.

    The axis values and ``frequency_Hz`` are added to every row
    automatically.

    Parameters
    ----------
    world_factory
        Callable that builds a fresh :class:`World` from the
        per-combination kwargs. Must return a fully populated
        world (soil, electrodes, sources, optional conductors).
    engine
        Either a static :class:`Engine` (reused for every
        combination) or a callable
        ``engine(**combination) -> Engine`` (rebuilt per
        combination — useful when the segment length depends on
        the geometry).
    axes
        Mapping ``axis_name -> Sequence`` of values. Must contain
        at least one axis. Empty sequences raise immediately.
    response
        Optional response extractor. Receives
        ``(result, world, frequency_index)`` and must return a
        ``dict[str, float|complex]``. Defaults to
        :func:`_default_response` (cluster impedance + EPR at
        the source's cluster).

    Returns
    -------
    pandas.DataFrame
        Long-format. Columns are the axis names, ``frequency_Hz``,
        and every key returned by ``response``. One row per
        Cartesian-product point per frequency.

    Raises
    ------
    ValueError
        If ``axes`` is empty or any axis is empty.
    """
    from groundfield.solver.engine import Engine

    if not axes:
        raise ValueError("axes must contain at least one axis.")
    for name, values in axes.items():
        if len(list(values)) == 0:
            raise ValueError(f"axis '{name}' is empty.")

    if response is None:
        response = _default_response

    keys = list(axes.keys())
    values_list = [list(axes[k]) for k in keys]

    rows: list[dict[str, Any]] = []
    for combo in itertools.product(*values_list):
        params = dict(zip(keys, combo))
        world = world_factory(**params)
        eng = engine if isinstance(engine, Engine) else engine(**params)
        result = eng.solve(world)
        for f_idx, f_hz in enumerate(result.frequencies):
            row = dict(params)
            row["frequency_Hz"] = float(f_hz)
            row.update(response(result, world, f_idx))
            rows.append(row)
    return pd.DataFrame(rows)

API reference — convergence

convergence

Convergence study over the engine's segment_length.

The PDE / field model in groundfield is a reference computation (see CLAUDE.md); to honour that role every non-trivial result should be backed up by a refinement study. This module turns the canonical "halve the segment length, watch what happens" experiment into one function call.

Mathematical / physical content

The image-family discretiser splits each electrode into segments of length :math:\Delta s. As :math:\Delta s \to 0 the multi-port grounding matrix approaches the continuous integral operator and the cluster impedance :math:Z_c converges to the PDE-grade reference. The convergence is monotone in :math:\Delta s for the average-potential method (cf. Sunde 1968; Tagg 1964). A practical refinement plot therefore answers two questions at once:

• "Has my chosen :math:\Delta s already converged within X %?" — the curve flattens out.
• "What is the asymptotic (PDE-grade) value?" — extrapolation / Richardson if needed.

Validity envelope

• Frequency: quasi-static envelope :math:f \le 1\,\mathrm{kHz}.
• Backends: any image-family solver (image / image_2layer / image_nlayer / mom / mom_sommerfeld / cim / bem). FEM is supported but its mesh is generated independently — the segment_length knob does not directly control its accuracy.
• The function clones the engine via :meth:Engine.model_copy, so the original engine is not mutated.

convergence_study

convergence_study(
    world: "World",
    engine: "Engine",
    *,
    segment_lengths: Sequence[float],
    response: (
        Callable[
            ["FieldResult", "World", int], dict[str, float]
        ]
        | None
    ) = None
) -> pd.DataFrame

Solve the same world repeatedly with refining segment_length.

For every :math:\Delta s_k in segment_lengths the function builds a clone of the engine (only segment_length differs), solves the world, and extracts a scalar response per frequency. The returned DataFrame is sorted descending in segment_length_m so finer resolutions land on the right in the default plot.

Parameters:

Name	Type	Description	Default
world	'World'	World to solve. Not modified.	required
engine	'Engine'	Base engine. Cloned per refinement step via :meth:Engine.model_copy; the original is left untouched.	required
segment_lengths	Sequence[float]	Sequence of segment lengths in metres. Must be strictly positive and contain at least two distinct values.	required
response	Callable[['FieldResult', 'World', int], dict[str, float]] | None	Optional response extractor. Defaults to :func:groundfield.postprocess.sweep._default_response (cluster impedance and EPR at the source's cluster).	None

Returns:

Type	Description
DataFrame	Long-format. Columns segment_length_m, frequency_Hz, n_segments and every key of the response (e.g. Z_re, Z_im, abs_Z, arg_Z_deg, U_E_re, ...).

Raises:

Type	Description
ValueError	If segment_lengths is empty, contains non-positive values, or has fewer than two distinct values (a convergence study below two refinement steps is meaningless).

Source code in src/groundfield/postprocess/convergence.py

def convergence_study(
    world: "World",
    engine: "Engine",
    *,
    segment_lengths: Sequence[float],
    response: Callable[["FieldResult", "World", int], dict[str, float]] | None = None,
) -> pd.DataFrame:
    """Solve the same world repeatedly with refining ``segment_length``.

    For every :math:`\\Delta s_k` in ``segment_lengths`` the
    function builds a clone of the engine (only ``segment_length``
    differs), solves the world, and extracts a scalar response
    per frequency. The returned DataFrame is sorted **descending**
    in ``segment_length_m`` so finer resolutions land on the right
    in the default plot.

    Parameters
    ----------
    world
        World to solve. **Not modified.**
    engine
        Base engine. Cloned per refinement step via
        :meth:`Engine.model_copy`; the original is left untouched.
    segment_lengths
        Sequence of segment lengths in metres. Must be strictly
        positive and contain at least two distinct values.
    response
        Optional response extractor. Defaults to
        :func:`groundfield.postprocess.sweep._default_response`
        (cluster impedance and EPR at the source's cluster).

    Returns
    -------
    pandas.DataFrame
        Long-format. Columns ``segment_length_m``,
        ``frequency_Hz``, ``n_segments`` and every key of the
        response (e.g. ``Z_re``, ``Z_im``, ``abs_Z``,
        ``arg_Z_deg``, ``U_E_re``, ...).

    Raises
    ------
    ValueError
        If ``segment_lengths`` is empty, contains non-positive
        values, or has fewer than two distinct values (a
        convergence study below two refinement steps is
        meaningless).
    """
    ds_list = [float(ds) for ds in segment_lengths]
    if len(ds_list) < 2:
        raise ValueError(
            f"segment_lengths must contain at least 2 values, got {len(ds_list)}."
        )
    if any(not (ds > 0) for ds in ds_list):
        raise ValueError(
            f"segment_lengths must be strictly positive, got {ds_list}."
        )
    if len(set(ds_list)) < 2:
        raise ValueError(
            "segment_lengths must contain at least 2 distinct values."
        )

    if response is None:
        response = _default_response

    rows: list[dict[str, Any]] = []
    for ds in ds_list:
        eng = engine.model_copy(update={"segment_length": ds})
        result = world.solve(eng)
        n_seg = len(result.point_sources)
        for f_idx, f_hz in enumerate(result.frequencies):
            row: dict[str, Any] = {
                "segment_length_m": ds,
                "frequency_Hz": float(f_hz),
                "n_segments": n_seg,
            }
            row.update(response(result, world, f_idx))
            rows.append(row)
    df = pd.DataFrame(rows).sort_values(
        ["frequency_Hz", "segment_length_m"], ascending=[True, False]
    ).reset_index(drop=True)
    return df

plot_convergence

plot_convergence(
    df: pd.DataFrame,
    *,
    response: str = "abs_Z",
    reference: float | None = None,
    figsize: tuple[float, float] = (7.5, 4.5),
    title: str | None = None
)

Plot response versus segment_length (one line per frequency).

The x-axis is logarithmic and inverted so finer resolutions sit on the right (the convergence "asymptote direction"). When the DataFrame contains a single frequency only, the legend is suppressed.

Parameters:

Name	Type	Description	Default
df	DataFrame	DataFrame returned by :func:convergence_study.	required
response	str	Column name to plot. Default "abs_Z".	'abs_Z'
reference	float | None	Optional asymptotic value to draw as a horizontal dashed reference line — useful when an analytical reference (Sunde, Dwight, IEEE Std 80) is known.	None
figsize	tuple[float, float]	Standard matplotlib options.	(7.5, 4.5)
title	tuple[float, float]	Standard matplotlib options.	(7.5, 4.5)

Returns:

Type	Description
Figure

Source code in src/groundfield/postprocess/convergence.py

def plot_convergence(
    df: pd.DataFrame,
    *,
    response: str = "abs_Z",
    reference: float | None = None,
    figsize: tuple[float, float] = (7.5, 4.5),
    title: str | None = None,
):
    """Plot ``response`` versus ``segment_length`` (one line per frequency).

    The x-axis is logarithmic and **inverted** so finer
    resolutions sit on the right (the convergence "asymptote
    direction"). When the DataFrame contains a single frequency
    only, the legend is suppressed.

    Parameters
    ----------
    df
        DataFrame returned by :func:`convergence_study`.
    response
        Column name to plot. Default ``"abs_Z"``.
    reference
        Optional asymptotic value to draw as a horizontal dashed
        reference line — useful when an analytical reference
        (Sunde, Dwight, IEEE Std 80) is known.
    figsize, title
        Standard matplotlib options.

    Returns
    -------
    matplotlib.figure.Figure
    """
    import matplotlib.pyplot as plt

    if response not in df.columns:
        raise KeyError(
            f"column '{response}' not in DataFrame; have: {list(df.columns)}"
        )
    if "segment_length_m" not in df.columns:
        raise KeyError(
            "DataFrame must have a 'segment_length_m' column "
            "(produced by convergence_study)."
        )

    fig, ax = plt.subplots(figsize=figsize)
    if "frequency_Hz" in df.columns and df["frequency_Hz"].nunique() > 1:
        for f_hz, sub in df.groupby("frequency_Hz"):
            sub = sub.sort_values("segment_length_m", ascending=False)
            ax.plot(
                sub["segment_length_m"], sub[response],
                marker="o", label=f"f = {f_hz:g} Hz",
            )
        ax.legend()
    else:
        sub = df.sort_values("segment_length_m", ascending=False)
        ax.plot(sub["segment_length_m"], sub[response], marker="o")

    if reference is not None:
        ax.axhline(
            float(reference), ls="--", color="k", alpha=0.6,
            label=f"reference = {reference:g}",
        )
        # Re-render legend so the reference line appears even in
        # the single-frequency case.
        ax.legend()

    ax.set_xscale("log")
    ax.invert_xaxis()  # finer ds (smaller value) on the right
    ax.set_xlabel("segment_length in m")
    ax.set_ylabel(response)
    ax.grid(True, which="both", alpha=0.3)
    if title is None:
        title = f"Convergence of {response} vs. segment_length"
    ax.set_title(title)
    fig.tight_layout()
    return fig

World-geometry plots (no solve required)

The :mod:groundfield.postprocess.geometry_plot module renders the physical world — electrodes, conductors and current sources — before the solver runs. production-grade networks with several hundred electrodes benefit from a quick sanity check that catches typos in coordinates, missing conductors, or sources attached to the wrong electrode without paying for a full solve.

Conductor colour scheme

The conductor colour follows :data:groundfield.conductors.ConductorType:

conductor_type	Colour
pen	green (#2c7a2c)
bare_copper	orange (#d97300)
cable_shield	grey (#888888)
overhead	steel blue (#1f77b4)
generic	dark grey (#444444)

The line style flags the soil-coupling mode: coupling_to_soil="galvanic" is solid, "isolated" is dashed.

Example

import groundfield as gf

world = gf.create_world(soil=gf.HomogeneousSoil(resistivity=100.0))
g_rod = gf.create_electrode(
    world, "rod", name="g1", position=(0.0, 0.0, 0.5), length=1.5,
)
g_ring = gf.create_electrode(
    world, "ring", name="g2", center=(10.0, 0.0, 0.8), radius=2.5,
)
gf.create_conductor(world, name="bond", start=g_rod, end=g_ring,
                    conductor_type="bare_copper")
gf.create_source(world, name="src", attached_to=g_rod, magnitude=10.0,
                 return_to=g_ring)

# Top-down 2-D view (default).
gf.plot_world(world, plane="xy")
# Vertical slice with the soil surface as a dotted grey line.
gf.plot_world(world, plane="xz")
# 3-D wireframe with inverted z (depth grows downwards on screen).
gf.plot_world_3d(world)

# Pure helper — bounding box of the geometry, useful for custom
# extents on the field plots.
x_min, x_max, y_min, y_max, z_min, z_max = gf.world_bounds_3d(world)

API reference — geometry plots

geometry_plot

World-geometry plots (no solve required).

This module provides quick visualisations of the physical world — electrodes, conductors and current sources — before the solver runs. In default setups with several hundred electrodes (200 EFH plus KVS plus substation plus measurement aux/probe) it is very useful to inspect the geometry first to catch typos in positions, accidental clusters, missing conductors, or sources attached to the wrong electrode, without paying the cost of a field solve.

Functions:

Name	Description
world_bounds_3d	Smallest axis-aligned bounding box in :math:(x, y, z) of the world's electrodes plus conductor endpoints. Extension of :func:groundfield.postprocess.plotting.world_bounds_xy to three dimensions.
plot_world	2-D top-down (plane="xy") or vertical (plane="xz") geometry plot. Electrodes drawn via the existing :func:_draw_electrodes helper of :mod:groundfield.postprocess.plotting; conductors as colour-coded line segments; sources as red star markers with an optional arrow to return_to.
plot_world_3d	3-D wireframe with the soil-surface plane shown in light grey and the :math:z-axis pointing downwards (groundfield convention; positive depth is below ground).

Conductor colour scheme

The conductor colour follows :data:groundfield.conductors.ConductorType:

================== ========================= ===================== conductor_type Colour Default style ================== ========================= ===================== pen #2c7a2c (green) solid bare_copper #d97300 (orange) solid cable_shield #888888 (grey) solid overhead #1f77b4 (steel blue) solid generic #444444 (dark grey) solid ================== ========================= =====================

The line style flags the soil-coupling mode of the conductor: coupling_to_soil = "galvanic" is drawn solid, "isolated" is drawn dashed.

Validity

Pure geometry — no solver result is required. The plot does not visualise any field quantity; for that, see :func:groundfield.postprocess.plotting.plot_potential_contour and friends. Best used as a debugging step before world.solve(...).

plot_world

plot_world(
    world: "World",
    *,
    plane: Literal["xy", "xz"] = "xy",
    extent: tuple[float, float, float, float] | None = None,
    padding_m: float = 5.0,
    show_conductors: bool = True,
    show_sources: bool = True,
    annotate_electrodes: bool = False,
    figsize: tuple[float, float] = (8.0, 6.0),
    ax=None,
    title: str | None = None
)

Top-down or vertical 2-D geometry plot of a :class:World.

Pure-geometry visualisation; no field quantity is evaluated and no solver is invoked. Useful as a sanity check before :meth:World.solve on a large typical network.

Parameters:

Name	Type	Description	Default
world	'World'	World to draw.	required
plane	Literal['xy', 'xz']	"xy" (default — top-down) or "xz" (side / vertical slice).	'xy'
extent	tuple[float, float, float, float] | None	Optional explicit (a_min, a_max, b_min, b_max) of the plotted area in metres. None (default) derives the extent from :func:world_bounds_3d plus padding_m.	None
padding_m	float	Extra padding on each side of the bounding box in metres (used only when extent is None).	5.0
show_conductors	bool	If True (default), draw every entry of :attr:World.conductors as a colour-coded line.	True
show_sources	bool	If True (default), draw every entry of :attr:World.sources as a red star at its anchor electrode, plus an arrow to return_to if set.	True
annotate_electrodes	bool	If True, attach a small text label with the electrode name to each electrode. Default False because production-grade worlds with > 50 electrodes look cluttered.	False
figsize	tuple[float, float]	Matplotlib figure size in inches (used when ax is None).	(8.0, 6.0)
ax		Optional pre-existing :class:matplotlib.axes.Axes. When passed, figsize is ignored and the host figure is returned.	None
title	str | None	Optional title override; default World 'name' — plane geometry.	None

Returns:

Type	Description
Figure

Source code in src/groundfield/postprocess/geometry_plot.py

def plot_world(
    world: "World",
    *,
    plane: Literal["xy", "xz"] = "xy",
    extent: tuple[float, float, float, float] | None = None,
    padding_m: float = 5.0,
    show_conductors: bool = True,
    show_sources: bool = True,
    annotate_electrodes: bool = False,
    figsize: tuple[float, float] = (8.0, 6.0),
    ax=None,
    title: str | None = None,
):
    """Top-down or vertical 2-D geometry plot of a :class:`World`.

    Pure-geometry visualisation; no field quantity is evaluated
    and no solver is invoked. Useful as a sanity check before
    :meth:`World.solve` on a large typical network.

    Parameters
    ----------
    world
        World to draw.
    plane
        ``"xy"`` (default — top-down) or ``"xz"`` (side / vertical
        slice).
    extent
        Optional explicit ``(a_min, a_max, b_min, b_max)`` of the
        plotted area in metres. ``None`` (default) derives the
        extent from :func:`world_bounds_3d` plus ``padding_m``.
    padding_m
        Extra padding on each side of the bounding box in metres
        (used only when ``extent`` is ``None``).
    show_conductors
        If ``True`` (default), draw every entry of
        :attr:`World.conductors` as a colour-coded line.
    show_sources
        If ``True`` (default), draw every entry of
        :attr:`World.sources` as a red star at its anchor
        electrode, plus an arrow to ``return_to`` if set.
    annotate_electrodes
        If ``True``, attach a small text label with the electrode
        name to each electrode. Default ``False`` because
        production-grade worlds with > 50 electrodes look cluttered.
    figsize
        Matplotlib figure size in inches (used when ``ax`` is
        ``None``).
    ax
        Optional pre-existing :class:`matplotlib.axes.Axes`. When
        passed, ``figsize`` is ignored and the host figure is
        returned.
    title
        Optional title override; default ``World 'name' — plane
        geometry``.

    Returns
    -------
    matplotlib.figure.Figure
    """
    import matplotlib.pyplot as plt

    from groundfield.postprocess.plotting import _draw_electrodes

    if extent is None:
        bounds = world_bounds_3d(world)
        x_min, x_max, y_min, y_max, z_min, z_max = bounds
        if plane == "xy":
            extent = (
                x_min - padding_m, x_max + padding_m,
                y_min - padding_m, y_max + padding_m,
            )
        elif plane == "xz":
            # In z, "padding" should not push the surface (z=0)
            # below the soil — but z<0 (overhead) is legitimate.
            # We pad symmetrically and rely on plot inversion.
            extent = (
                x_min - padding_m, x_max + padding_m,
                z_min - padding_m, z_max + padding_m,
            )
        else:
            raise ValueError(f"plane must be 'xy' or 'xz', got {plane!r}.")

    if ax is None:
        fig, ax = plt.subplots(figsize=figsize)
    else:
        fig = ax.figure

    # 1. Surface line for the xz view — emphasises the soil
    #    boundary at z = 0.
    if plane == "xz":
        ax.axhline(0.0, color="#666666", linestyle=":", linewidth=0.8,
                   alpha=0.7, zorder=1)

    # 2. Conductors first (so electrodes draw on top).
    if show_conductors:
        for c in world.conductors:
            _draw_conductor_2d(ax, c, plane)

    # 3. Electrodes — reuse the existing helper for consistency
    #    with the field plots.
    _draw_electrodes(ax, world, plane)

    # 4. Sources on top.
    if show_sources:
        for s in world.sources:
            _draw_source_2d(ax, s, world, plane)

    # 5. Optional annotations.
    if annotate_electrodes:
        for e in world.electrodes:
            cp = e.connection_point
            x, y = _project_2d(cp, plane)
            ax.annotate(
                e.name, (x, y),
                textcoords="offset points", xytext=(5, 5),
                fontsize=7, alpha=0.85, zorder=6,
            )

    ax.set_xlim(extent[0], extent[1])
    ax.set_ylim(extent[2], extent[3])
    ax.set_aspect("equal")
    ax.grid(True, alpha=0.3)
    if plane == "xy":
        ax.set_xlabel("x in m")
        ax.set_ylabel("y in m")
    else:
        ax.set_xlabel("x in m")
        ax.set_ylabel("z in m (depth)")
        ax.invert_yaxis()  # soil downwards

    if title is None:
        n_e = len(world.electrodes)
        n_c = len(world.conductors)
        n_s = len(world.sources)
        title = (
            f"World '{world.name}' — {plane} geometry "
            f"({n_e} electrodes, {n_c} conductors, {n_s} sources)"
        )
    ax.set_title(title)

    _build_legend(ax, world, show_conductors, show_sources)
    fig.tight_layout()
    return fig

plot_world_3d

plot_world_3d(
    world: "World",
    *,
    show_conductors: bool = True,
    show_sources: bool = True,
    show_surface: bool = True,
    figsize: tuple[float, float] = (9.0, 7.0),
    elev: float = 22.0,
    azim: float = -55.0,
    title: str | None = None
)

3-D wireframe of a :class:World using mpl_toolkits.mplot3d.

The :math:z axis is inverted so that depth points downwards on screen (groundfield convention: positive :math:z is into the soil). A faint grey surface plane at :math:z = 0 marks the soil surface.

Parameters:

Name	Type	Description	Default
world	'World'	World to draw.	required
show_conductors	bool	See :func:plot_world.	True
show_sources	bool	See :func:plot_world.	True
show_surface	bool	If True (default), render a translucent grey square at :math:z = 0 over the world's :math:(x, y) bounding box plus 5 m padding.	True
figsize	tuple[float, float]	Matplotlib figure size in inches.	(9.0, 7.0)
elev	float	Initial viewing angles (passed to :meth:Axes3D.view_init).	22.0
azim	float	Initial viewing angles (passed to :meth:Axes3D.view_init).	22.0
title	str | None	Optional title override.	None

Returns:

Type	Description
Figure

Source code in src/groundfield/postprocess/geometry_plot.py

def plot_world_3d(
    world: "World",
    *,
    show_conductors: bool = True,
    show_sources: bool = True,
    show_surface: bool = True,
    figsize: tuple[float, float] = (9.0, 7.0),
    elev: float = 22.0,
    azim: float = -55.0,
    title: str | None = None,
):
    """3-D wireframe of a :class:`World` using ``mpl_toolkits.mplot3d``.

    The :math:`z` axis is **inverted** so that depth points
    downwards on screen (groundfield convention: positive
    :math:`z` is into the soil). A faint grey surface plane at
    :math:`z = 0` marks the soil surface.

    Parameters
    ----------
    world
        World to draw.
    show_conductors, show_sources
        See :func:`plot_world`.
    show_surface
        If ``True`` (default), render a translucent grey square
        at :math:`z = 0` over the world's :math:`(x, y)` bounding
        box plus 5 m padding.
    figsize
        Matplotlib figure size in inches.
    elev, azim
        Initial viewing angles (passed to
        :meth:`Axes3D.view_init`).
    title
        Optional title override.

    Returns
    -------
    matplotlib.figure.Figure
    """
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D  # noqa: F401  (registers 3d projection)

    fig = plt.figure(figsize=figsize)
    ax = fig.add_subplot(111, projection="3d")

    bounds = world_bounds_3d(world)
    x_min, x_max, y_min, y_max, z_min, z_max = bounds

    # Optional surface plane at z = 0.
    if show_surface and (world.electrodes or world.conductors):
        pad = 5.0
        xs = np.array([[x_min - pad, x_max + pad],
                       [x_min - pad, x_max + pad]])
        ys = np.array([[y_min - pad, y_min - pad],
                       [y_max + pad, y_max + pad]])
        zs = np.zeros_like(xs)
        ax.plot_surface(xs, ys, zs, color="#cccccc", alpha=0.25,
                        edgecolor="none", linewidth=0)

    if show_conductors:
        for c in world.conductors:
            _draw_conductor_3d(ax, c)

    for e in world.electrodes:
        _draw_electrode_3d(ax, e)

    if show_sources:
        for s in world.sources:
            _draw_source_3d(ax, s, world)

    ax.set_xlabel("x in m")
    ax.set_ylabel("y in m")
    ax.set_zlabel("z in m (depth)")
    # Invert z so that depth points downwards on screen.
    ax.invert_zaxis()
    ax.view_init(elev=elev, azim=azim)
    if title is None:
        title = (
            f"World '{world.name}' — 3-D geometry "
            f"({len(world.electrodes)} electrodes, "
            f"{len(world.conductors)} conductors, "
            f"{len(world.sources)} sources)"
        )
    ax.set_title(title)
    fig.tight_layout()
    return fig

world_bounds_3d

world_bounds_3d(
    world: "World",
) -> tuple[float, float, float, float, float, float]

Smallest axis-aligned :math:(x, y, z) bounding box of the world.

Inspects every electrode in :attr:World.electrodes and every conductor endpoint in :attr:World.conductors, returning the six bounds (x_min, x_max, y_min, y_max, z_min, z_max) in metres.

The :math:(x, y) portion is consistent with :func:groundfield.postprocess.plotting.world_bounds_xy but additionally includes conductor endpoints — important for overhead lines or measurement leads that extend well beyond the electrode footprint. The :math:z portion uses

• :class:RodElectrode: [position[2], position[2] + length] (head + foot of the rod).
• :class:RingElectrode, :class:StripElectrode, :class:PolylineElectrode, :class:MeshElectrode, :class:GridMeshElectrode: the :math:z coordinate(s) of the electrode (single buried depth).
• Conductors: min/max over both endpoint :math:z values (catches overhead :math:z<0 and buried :math:z>0).

Returns:

Type	Description
tuple	(x_min, x_max, y_min, y_max, z_min, z_max) in metres. For an empty world the result is the trivial (0, 0, 0, 0, 0, 0); callers should add positive padding before plotting.

Source code in src/groundfield/postprocess/geometry_plot.py

def world_bounds_3d(
    world: "World",
) -> tuple[float, float, float, float, float, float]:
    """Smallest axis-aligned :math:`(x, y, z)` bounding box of the world.

    Inspects every electrode in :attr:`World.electrodes` and every
    conductor endpoint in :attr:`World.conductors`, returning the
    six bounds ``(x_min, x_max, y_min, y_max, z_min, z_max)`` in
    metres.

    The :math:`(x, y)` portion is consistent with
    :func:`groundfield.postprocess.plotting.world_bounds_xy` but
    additionally includes conductor endpoints — important for
    overhead lines or measurement leads that extend well beyond
    the electrode footprint. The :math:`z` portion uses

    * :class:`RodElectrode`: ``[position[2], position[2] + length]``
      (head + foot of the rod).
    * :class:`RingElectrode`, :class:`StripElectrode`,
      :class:`PolylineElectrode`, :class:`MeshElectrode`,
      :class:`GridMeshElectrode`: the :math:`z` coordinate(s) of the
      electrode (single buried depth).
    * Conductors: ``min/max`` over both endpoint :math:`z` values
      (catches overhead :math:`z<0` and buried :math:`z>0`).

    Returns
    -------
    tuple
        ``(x_min, x_max, y_min, y_max, z_min, z_max)`` in metres.
        For an empty world the result is the trivial
        ``(0, 0, 0, 0, 0, 0)``; callers should add positive
        padding before plotting.
    """
    from groundfield.geometry.electrodes import (
        GridMeshElectrode,
        MeshElectrode,
        PolylineElectrode,
        RingElectrode,
        RodElectrode,
        StarElectrode,
        StripElectrode,
    )

    if not world.electrodes and not world.conductors:
        return (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

    xs: list[float] = []
    ys: list[float] = []
    zs: list[float] = []

    for e in world.electrodes:
        if isinstance(e, RodElectrode):
            x, y, z = e.position
            xs.append(x)
            ys.append(y)
            zs.extend([z, z + e.length])
        elif isinstance(e, RingElectrode):
            cx, cy, cz = e.center
            xs.extend([cx - e.radius, cx + e.radius])
            ys.extend([cy - e.radius, cy + e.radius])
            zs.append(cz)
        elif isinstance(e, StripElectrode):
            xs.extend([e.start[0], e.end[0]])
            ys.extend([e.start[1], e.end[1]])
            zs.extend([e.start[2], e.end[2]])
        elif isinstance(e, (MeshElectrode, GridMeshElectrode)):
            cx, cy, cz = e.corner
            dx, dy = e.size
            xs.extend([cx, cx + dx])
            ys.extend([cy, cy + dy])
            zs.append(cz)
        elif isinstance(e, PolylineElectrode):
            xs.extend(v[0] for v in e.vertices)
            ys.extend(v[1] for v in e.vertices)
            zs.extend(v[2] for v in e.vertices)
        elif isinstance(e, StarElectrode):
            for start, end in e.edges:
                xs.extend([start[0], end[0]])
                ys.extend([start[1], end[1]])
                zs.extend([start[2], end[2]])
        else:  # pragma: no cover — defensive
            cp = e.connection_point
            xs.append(cp[0])
            ys.append(cp[1])
            zs.append(cp[2])

    for c in world.conductors:
        for endpoint in (c.start, c.end):
            xs.append(endpoint[0])
            ys.append(endpoint[1])
            zs.append(endpoint[2])

    return (min(xs), max(xs), min(ys), max(ys), min(zs), max(zs))

Plotting API reference

plotting

Matplotlib plots for potential distributions and profiles.

Functions:

Name	Description
plot_potential_contour	Contour / pseudo-colour plot of the potential in a slice plane (xy at fixed \(z\) or xz at fixed \(y\)).
plot_potential_profile	Line plot of the potential along an arbitrary direction for one or several \(z\) depths.
plot_potential_radial	Radial profile \(\varphi(r)\) starting from an electrode for several depths — the standard "how far does the trumpet reach?" plot.
All functions return the ``matplotlib.figure.Figure`` they produce so
that notebooks can post-process them

plot_potential_contour

plot_potential_contour(
    result: "FieldResult",
    *,
    world: "World | None" = None,
    plane: Literal["xy", "xz"] = "xy",
    z: float = 0.0,
    y: float = 0.0,
    extent: tuple[float, float, float, float] | None = None,
    n: int = 120,
    frequency_index: int = 0,
    levels: int = 20,
    log: bool = False,
    cmap: str = "viridis"
)

Contour plot of the potential in a slice plane.

Parameters:

Name	Type	Description	Default
result	'FieldResult'	Result object from :meth:Engine.solve.	required
world	'World | None'	Optional companion world; if given, electrodes are drawn into the plot.	None
plane	Literal['xy', 'xz']	"xy" (horizontal slice at depth \(z\)) or "xz" (vertical slice at \(y\)).	'xy'
z	float	Depth or \(y\) value of the slice plane in metres.	0.0
y	float	Depth or \(y\) value of the slice plane in metres.	0.0
extent	tuple[float, float, float, float] | None	(a_min, a_max, b_min, b_max) of the slice plane in metres. Default is derived from the point-source distribution.	None
n	int	Resolution per axis.	120
frequency_index	int	Index into :attr:FieldResult.frequencies.	0
levels	int	Number of contour levels.	20
log	bool	Logarithmic colour scale (better for fields that span several decades).	False
cmap	str	Matplotlib colormap name.	'viridis'

Returns:

Name	Type	Description
fig	Figure

Source code in src/groundfield/postprocess/plotting.py

def plot_potential_contour(
    result: "FieldResult",
    *,
    world: "World | None" = None,
    plane: Literal["xy", "xz"] = "xy",
    z: float = 0.0,
    y: float = 0.0,
    extent: tuple[float, float, float, float] | None = None,
    n: int = 120,
    frequency_index: int = 0,
    levels: int = 20,
    log: bool = False,
    cmap: str = "viridis",
):
    """Contour plot of the potential in a slice plane.

    Parameters
    ----------
    result
        Result object from :meth:`Engine.solve`.
    world
        Optional companion world; if given, electrodes are drawn into
        the plot.
    plane
        ``"xy"`` (horizontal slice at depth $z$) or ``"xz"``
        (vertical slice at $y$).
    z, y
        Depth or $y$ value of the slice plane in metres.
    extent
        ``(a_min, a_max, b_min, b_max)`` of the slice plane in metres.
        Default is derived from the point-source distribution.
    n
        Resolution per axis.
    frequency_index
        Index into :attr:`FieldResult.frequencies`.
    levels
        Number of contour levels.
    log
        Logarithmic colour scale (better for fields that span several
        decades).
    cmap
        Matplotlib colormap name.

    Returns
    -------
    fig : matplotlib.figure.Figure
    """
    import matplotlib.pyplot as plt

    if extent is None:
        srcs = np.array([ps.position for ps in result.point_sources])
        if plane == "xy":
            ax = srcs[:, 0]
            bx = srcs[:, 1]
        else:
            ax = srcs[:, 0]
            bx = srcs[:, 2]
        pad = max(5.0, 0.5 * (ax.ptp() + bx.ptp()))
        extent = (ax.min() - pad, ax.max() + pad,
                  bx.min() - pad, bx.max() + pad)

    fixed = z if plane == "xy" else y
    A, B, flat = _make_grid(plane, extent, fixed, n)

    phi = result.potential(flat, frequency_index=frequency_index).real
    phi = phi.reshape(A.shape)

    fig, axx = plt.subplots(figsize=(7, 5))
    if log:
        # Drop non-positive values, otherwise the log scale breaks.
        phi_plot = np.where(phi > 0, phi, np.nan)
        cs = axx.contourf(A, B, phi_plot, levels=levels, cmap=cmap,
                          locator=__import__("matplotlib.ticker",
                                             fromlist=["LogLocator"]).LogLocator())
    else:
        cs = axx.contourf(A, B, phi, levels=levels, cmap=cmap)
    contours = axx.contour(A, B, phi, levels=levels, colors="k",
                           linewidths=0.4, alpha=0.5)
    axx.clabel(contours, inline=True, fontsize=7, fmt="%.0f")

    cbar = fig.colorbar(cs, ax=axx)
    cbar.set_label("Potential φ in V")

    if plane == "xy":
        axx.set_xlabel("x in m")
        axx.set_ylabel("y in m")
        axx.set_title(f"Potential φ(x, y, z={z:g} m), "
                      f"f={result.frequencies[frequency_index]} Hz")
        axx.set_aspect("equal")
    else:
        axx.set_xlabel("x in m")
        axx.set_ylabel("z in m (depth)")
        axx.invert_yaxis()  # soil downwards
        axx.set_title(f"Potential φ(x, y={y:g} m, z), "
                      f"f={result.frequencies[frequency_index]} Hz")

    if world is not None:
        _draw_electrodes(axx, world, plane)

    fig.tight_layout()
    return fig

plot_potential_profile

plot_potential_profile(
    result: "FieldResult",
    *,
    start: tuple[float, float, float],
    direction: tuple[float, float, float] = (1.0, 0.0, 0.0),
    distance: float = 30.0,
    n: int = 200,
    depths: Iterable[float] | None = None,
    frequency_index: int = 0
)

Potential along an arbitrary line, optionally for several depths.

Parameters:

Name	Type	Description	Default
result	'FieldResult'	Result object.	required
start	tuple[float, float, float]	Start point (x, y, z) of the line in metres. z is overridden by depths if given.	required
direction	tuple[float, float, float]	Direction vector (will be normalised).	(1.0, 0.0, 0.0)
distance	float	Length of the line in metres.	30.0
n	int	Number of evaluation points.	200
depths	Iterable[float] | None	List of depths \(z\) in metres. One curve per depth.	None
frequency_index	int	Index into :attr:FieldResult.frequencies.	0

Returns:

Name	Type	Description
fig	Figure

Source code in src/groundfield/postprocess/plotting.py

def plot_potential_profile(
    result: "FieldResult",
    *,
    start: tuple[float, float, float],
    direction: tuple[float, float, float] = (1.0, 0.0, 0.0),
    distance: float = 30.0,
    n: int = 200,
    depths: Iterable[float] | None = None,
    frequency_index: int = 0,
):
    """Potential along an arbitrary line, optionally for several depths.

    Parameters
    ----------
    result
        Result object.
    start
        Start point ``(x, y, z)`` of the line in metres. ``z`` is
        overridden by ``depths`` if given.
    direction
        Direction vector (will be normalised).
    distance
        Length of the line in metres.
    n
        Number of evaluation points.
    depths
        List of depths $z$ in metres. One curve per depth.
    frequency_index
        Index into :attr:`FieldResult.frequencies`.

    Returns
    -------
    fig : matplotlib.figure.Figure
    """
    import matplotlib.pyplot as plt

    if depths is None:
        depths = [0.0]
    direction_arr = np.asarray(direction, dtype=float)
    direction_arr = direction_arr / np.linalg.norm(direction_arr)
    s = np.linspace(0.0, distance, n)

    fig, ax = plt.subplots(figsize=(7, 4))
    for z_d in depths:
        pts = np.column_stack(
            [start[0] + s * direction_arr[0],
             start[1] + s * direction_arr[1],
             np.full_like(s, z_d)]
        )
        phi = result.potential(pts, frequency_index=frequency_index).real
        ax.plot(s, phi, label=f"z = {z_d:g} m")

    ax.set_xlabel("distance s along profile in m")
    ax.set_ylabel("Potential φ in V")
    ax.set_title(
        f"Potential profile from {start}, "
        f"f = {result.frequencies[frequency_index]} Hz"
    )
    ax.grid(True, alpha=0.3)
    ax.legend()
    fig.tight_layout()
    return fig

plot_potential_radial

plot_potential_radial(
    result: "FieldResult",
    *,
    around: str | tuple[float, float, float],
    world: "World | None" = None,
    r_max: float = 30.0,
    n: int = 200,
    depths: Iterable[float] = (0.0, 0.5, 1.0),
    frequency_index: int = 0,
    log_x: bool = False
)

Trumpet-shape decay around an electrode (or fixed point).

Evaluates \(\varphi(r)\) along the +x direction starting at the connection point of the given electrode (or at a fixed point). Multiple depths produce multiple curves — the typical "surface vs. deeper soil" comparison.

Parameters:

Name	Type	Description	Default
around	str | tuple[float, float, float]	Electrode name (resolved via world) or (x, y, z) tuple.	required
world	'World | None'	Required when around is an electrode name.	None
r_max	float	Maximum radius in metres.	30.0
n	int	Number of evaluation points.	200
depths	Iterable[float]	\(z\) depths for the comparison curves.	(0.0, 0.5, 1.0)
frequency_index	int	Index into :attr:FieldResult.frequencies.	0
log_x	bool	If True use a logarithmic x-axis (better visibility of the trumpet shape).	False

Returns:

Name	Type	Description
fig	Figure

Source code in src/groundfield/postprocess/plotting.py

def plot_potential_radial(
    result: "FieldResult",
    *,
    around: str | tuple[float, float, float],
    world: "World | None" = None,
    r_max: float = 30.0,
    n: int = 200,
    depths: Iterable[float] = (0.0, 0.5, 1.0),
    frequency_index: int = 0,
    log_x: bool = False,
):
    """Trumpet-shape decay around an electrode (or fixed point).

    Evaluates $\\varphi(r)$ along the ``+x`` direction starting
    at the connection point of the given electrode (or at a fixed
    point). Multiple ``depths`` produce multiple curves — the typical
    "surface vs. deeper soil" comparison.

    Parameters
    ----------
    around
        Electrode name (resolved via ``world``) or ``(x, y, z)`` tuple.
    world
        Required when ``around`` is an electrode name.
    r_max
        Maximum radius in metres.
    n
        Number of evaluation points.
    depths
        $z$ depths for the comparison curves.
    frequency_index
        Index into :attr:`FieldResult.frequencies`.
    log_x
        If ``True`` use a logarithmic x-axis (better visibility of the
        trumpet shape).

    Returns
    -------
    fig : matplotlib.figure.Figure
    """
    import matplotlib.pyplot as plt

    if isinstance(around, str):
        if world is None:
            raise ValueError(
                "When 'around' is an electrode name, 'world' is required."
            )
        cp = world.get_electrode(around).connection_point
        x0, y0 = cp[0], cp[1]
        label_origin = around
    else:
        x0, y0 = around[0], around[1]
        label_origin = f"{around}"

    rs = np.linspace(0.1 if not log_x else 0.05, r_max, n)
    fig, ax = plt.subplots(figsize=(7, 4))
    for z_d in depths:
        pts = np.column_stack(
            [x0 + rs, np.full_like(rs, y0), np.full_like(rs, z_d)]
        )
        phi = result.potential(pts, frequency_index=frequency_index).real
        ax.plot(rs, phi, label=f"z = {z_d:g} m")

    ax.set_xlabel(f"distance r from {label_origin} in m")
    ax.set_ylabel("Potential φ in V")
    ax.set_title(
        f"Radial profile around {label_origin}, "
        f"f = {result.frequencies[frequency_index]} Hz"
    )
    if log_x:
        ax.set_xscale("log")
    ax.grid(True, which="both", alpha=0.3)
    ax.legend()
    fig.tight_layout()
    return fig

plot_surface_potential

plot_surface_potential(
    result: "FieldResult",
    world: "World",
    *,
    z: float = 0.0,
    padding_m: float = 15.0,
    n: int = 200,
    extent: tuple[float, float, float, float] | None = None,
    frequency_index: int = 0,
    levels: int = 25,
    log: bool = False,
    symmetric: bool = False,
    cmap: str = "viridis",
    show_electrodes: bool = True,
    show_contour_lines: bool = True,
    figsize: tuple[float, float] = (8.5, 7.0),
    title: str | None = None
)

Surface-potential pseudo-colour plot over the entire world.

Differs from :func:plot_potential_contour in two ways:

1. The default extent is derived from the world's electrode bounding box (via :func:world_bounds_xy) plus padding_m, not from the discretised current point sources. This makes the plot naturally cover all buildings, cable cabinets and the substation in a TN network — including the "boundary regions" where the potential decays back to remote earth.
2. The plot is locked to a horizontal slice (always plane="xy"); this is what typical calls the surface potential. Use :func:plot_potential_contour for vertical slices.

Mathematical content

For each grid point \((x_i, y_j, z)\) the function evaluates

.. math::

\varphi(x_i, y_j, z) \;=\; \frac{1}{4\pi\sigma_0}
\sum_k \frac{I_k}{|\mathbf{r}_{ij} - \mathbf{r}_k|}
\;+\; \text{image-charge series}

via :meth:FieldResult.potential, which dispatches to the correct Green's function for the soil model the result was computed with (homogeneous / 2-layer / multi-layer). The real part is plotted; for log=True the colour scale uses \(\log_{10} |\varphi|\) so the boundary decay is visible across several decades.

Parameters:

Name	Type	Description	Default
result	'FieldResult'	Solver output from :meth:Engine.solve.	required
world	'World'	Companion world used both to derive the plot extent and (optionally) to overlay electrodes.	required
z	float	Depth of the slice in metres. 0.0 is the ground surface; positive values go into the soil.	0.0
padding_m	float	Extra space added on each side of the electrode bounding box, in metres. Larger values let you see how the potential decays towards remote earth.	15.0
n	int	Resolution per axis (the grid is n × n).	200
extent	tuple[float, float, float, float] | None	Optional explicit (x_min, x_max, y_min, y_max) override. None (default) uses world_bounds_xy(world) plus padding_m.	None
frequency_index	int	Index into :attr:FieldResult.frequencies.	0
levels	int	Number of contour fill levels.	25
log	bool	Logarithmic colour scale based on |φ| (better for the boundary decay).	False
symmetric	bool	Use a symmetric [-φ_max, +φ_max] colour scale around 0. Useful when the potential takes both signs (multi-source or net-injection-zero studies).	False
cmap	str	Matplotlib colormap name.	'viridis'
show_electrodes	bool	Draw the electrode geometry on top of the contour.	True
show_contour_lines	bool	Draw black iso-potential lines on top of the fill.	True
figsize	tuple[float, float]	Matplotlib figure size in inches.	(8.5, 7.0)
title	str | None	Optional title override; if None a sensible default is constructed from z and the chosen frequency.	None

Returns:

Name	Type	Description
fig	Figure

Source code in src/groundfield/postprocess/plotting.py

def plot_surface_potential(
    result: "FieldResult",
    world: "World",
    *,
    z: float = 0.0,
    padding_m: float = 15.0,
    n: int = 200,
    extent: tuple[float, float, float, float] | None = None,
    frequency_index: int = 0,
    levels: int = 25,
    log: bool = False,
    symmetric: bool = False,
    cmap: str = "viridis",
    show_electrodes: bool = True,
    show_contour_lines: bool = True,
    figsize: tuple[float, float] = (8.5, 7.0),
    title: str | None = None,
):
    """Surface-potential pseudo-colour plot over the entire world.

    Differs from :func:`plot_potential_contour` in two ways:

    1. The default ``extent`` is derived from the *world's
       electrode bounding box* (via :func:`world_bounds_xy`) plus
       ``padding_m``, not from the discretised current point
       sources. This makes the plot naturally cover all
       buildings, cable cabinets and the substation in a TN
       network — including the "boundary regions" where the
       potential decays back to remote earth.
    2. The plot is locked to a horizontal slice (always
       ``plane="xy"``); this is what typical calls the *surface
       potential*. Use :func:`plot_potential_contour` for vertical
       slices.

    Mathematical content
    --------------------
    For each grid point $(x_i, y_j, z)$ the function evaluates

    .. math::

        \\varphi(x_i, y_j, z) \\;=\\; \\frac{1}{4\\pi\\sigma_0}
        \\sum_k \\frac{I_k}{|\\mathbf{r}_{ij} - \\mathbf{r}_k|}
        \\;+\\; \\text{image-charge series}

    via :meth:`FieldResult.potential`, which dispatches to the
    correct Green's function for the soil model the result was
    computed with (homogeneous / 2-layer / multi-layer). The
    real part is plotted; for ``log=True`` the colour scale uses
    $\\log_{10} |\\varphi|$ so the boundary decay is visible across
    several decades.

    Parameters
    ----------
    result
        Solver output from :meth:`Engine.solve`.
    world
        Companion world used both to derive the plot extent and
        (optionally) to overlay electrodes.
    z
        Depth of the slice in metres. ``0.0`` is the ground
        surface; positive values go into the soil.
    padding_m
        Extra space added on each side of the electrode bounding
        box, in metres. Larger values let you see how the
        potential decays towards remote earth.
    n
        Resolution per axis (the grid is ``n × n``).
    extent
        Optional explicit ``(x_min, x_max, y_min, y_max)``
        override. ``None`` (default) uses
        ``world_bounds_xy(world)`` plus ``padding_m``.
    frequency_index
        Index into :attr:`FieldResult.frequencies`.
    levels
        Number of contour fill levels.
    log
        Logarithmic colour scale based on ``|φ|`` (better for the
        boundary decay).
    symmetric
        Use a symmetric ``[-φ_max, +φ_max]`` colour scale around 0.
        Useful when the potential takes both signs (multi-source
        or net-injection-zero studies).
    cmap
        Matplotlib colormap name.
    show_electrodes
        Draw the electrode geometry on top of the contour.
    show_contour_lines
        Draw black iso-potential lines on top of the fill.
    figsize
        Matplotlib figure size in inches.
    title
        Optional title override; if ``None`` a sensible default is
        constructed from ``z`` and the chosen frequency.

    Returns
    -------
    fig : matplotlib.figure.Figure
    """
    import matplotlib.pyplot as plt
    from matplotlib.colors import LogNorm, Normalize

    if extent is None:
        x_min, x_max, y_min, y_max = world_bounds_xy(world)
        extent = (
            x_min - padding_m, x_max + padding_m,
            y_min - padding_m, y_max + padding_m,
        )

    A, B, flat = _make_grid("xy", extent, z, n)
    phi = result.potential(flat, frequency_index=frequency_index).real
    phi = phi.reshape(A.shape)

    fig, ax = plt.subplots(figsize=figsize)

    if log:
        phi_abs = np.abs(phi)
        # Avoid log(0) by clipping at the smallest positive value;
        # mask-zero so the contourf renders the under-range as bg.
        positive = phi_abs[phi_abs > 0]
        v_min = float(positive.min()) if positive.size else 1e-9
        v_max = float(phi_abs.max()) if phi_abs.size else 1.0
        norm = LogNorm(vmin=max(v_min, 1e-9), vmax=max(v_max, 10 * v_min))
        cs = ax.contourf(A, B, phi_abs, levels=levels, cmap=cmap, norm=norm)
        cbar_label = "|φ| in V (log scale)"
    elif symmetric:
        v = float(np.max(np.abs(phi)))
        norm = Normalize(vmin=-v, vmax=v)
        cs = ax.contourf(A, B, phi, levels=levels, cmap=cmap, norm=norm)
        cbar_label = "Potential φ in V"
    else:
        cs = ax.contourf(A, B, phi, levels=levels, cmap=cmap)
        cbar_label = "Potential φ in V"

    if show_contour_lines and not log:
        contours = ax.contour(A, B, phi, levels=levels, colors="k",
                              linewidths=0.4, alpha=0.4)
        ax.clabel(contours, inline=True, fontsize=7, fmt="%.0f")
    elif show_contour_lines and log:
        # Log mode: contour lines on the absolute value, fewer levels
        # so labels stay legible.
        n_lines = max(5, levels // 3)
        contours = ax.contour(A, B, np.abs(phi), levels=n_lines,
                              colors="k", linewidths=0.4, alpha=0.4)

    cbar = fig.colorbar(cs, ax=ax)
    cbar.set_label(cbar_label)

    ax.set_xlabel("x in m")
    ax.set_ylabel("y in m")
    ax.set_aspect("equal")
    if title is None:
        f_hz = result.frequencies[frequency_index]
        title = (
            f"Surface potential φ(x, y, z = {z:g} m), "
            f"f = {f_hz:g} Hz"
        )
    ax.set_title(title)

    if show_electrodes:
        _draw_electrodes(ax, world, "xy")

    fig.tight_layout()
    return fig

world_bounds_xy

world_bounds_xy(
    world: "World",
) -> tuple[float, float, float, float]

Compute the horizontal bounding box of a world's electrodes.

Inspects every electrode in world.electrodes and returns the smallest axis-aligned rectangle in the \((x, y)\) plane that contains the geometry. Used by :func:plot_surface_potential as the natural "the whole world" extent for a surface contour plot.

The footprint of each electrode kind is treated explicitly:

• :class:RodElectrode — point at position[:2].
• :class:RingElectrode — square of side \(2r\) around the centre.
• :class:StripElectrode — bounding box of start and end.
• :class:MeshElectrode / :class:GridMeshElectrode — corner to corner + size.

Returns:

Type	Description
tuple	(x_min, x_max, y_min, y_max) in metres. For a world with no electrodes the result is the trivial (0, 0, 0, 0); callers should add positive padding before plotting.

Source code in src/groundfield/postprocess/plotting.py

def world_bounds_xy(world: "World") -> tuple[float, float, float, float]:
    """Compute the horizontal bounding box of a world's electrodes.

    Inspects every electrode in ``world.electrodes`` and returns
    the smallest axis-aligned rectangle in the $(x, y)$ plane that
    contains the geometry. Used by :func:`plot_surface_potential` as
    the natural "the whole world" extent for a surface contour plot.

    The footprint of each electrode kind is treated explicitly:

    * :class:`RodElectrode` — point at ``position[:2]``.
    * :class:`RingElectrode` — square of side $2r$ around the centre.
    * :class:`StripElectrode` — bounding box of ``start`` and ``end``.
    * :class:`MeshElectrode` / :class:`GridMeshElectrode` —
      ``corner`` to ``corner + size``.

    Returns
    -------
    tuple
        ``(x_min, x_max, y_min, y_max)`` in metres. For a world with
        no electrodes the result is the trivial ``(0, 0, 0, 0)``;
        callers should add positive padding before plotting.
    """
    from groundfield.geometry.electrodes import (
        GridMeshElectrode,
        MeshElectrode,
        PolylineElectrode,
        RingElectrode,
        RodElectrode,
        StarElectrode,
        StripElectrode,
    )

    if not world.electrodes:
        return (0.0, 0.0, 0.0, 0.0)

    xs: list[float] = []
    ys: list[float] = []
    for e in world.electrodes:
        if isinstance(e, RodElectrode):
            xs.append(e.position[0])
            ys.append(e.position[1])
        elif isinstance(e, RingElectrode):
            cx, cy, _ = e.center
            xs.extend([cx - e.radius, cx + e.radius])
            ys.extend([cy - e.radius, cy + e.radius])
        elif isinstance(e, StripElectrode):
            xs.extend([e.start[0], e.end[0]])
            ys.extend([e.start[1], e.end[1]])
        elif isinstance(e, (MeshElectrode, GridMeshElectrode)):
            cx, cy, _ = e.corner
            dx, dy = e.size
            xs.extend([cx, cx + dx])
            ys.extend([cy, cy + dy])
        elif isinstance(e, PolylineElectrode):
            xs.extend(v[0] for v in e.vertices)
            ys.extend(v[1] for v in e.vertices)
        elif isinstance(e, StarElectrode):
            for start, end in e.edges:
                xs.extend([start[0], end[0]])
                ys.extend([start[1], end[1]])
        else:  # pragma: no cover — defensive against future kinds
            cp = e.connection_point
            xs.append(cp[0])
            ys.append(cp[1])
    return (min(xs), max(xs), min(ys), max(ys))