`cim` — Complex Image Method¶

Status: reduced to closed-form kernels ($n \le 2$ only)

Since 0.11.0 cim rejects $n \ge 3$ soils with NotImplementedError (the historic complex-image kernel was structurally incomplete — audit 2026-07-08, WP-E), and since 0.15.0 it no longer computes a complex-image fit at all. What runs is:

Soil	Kernel	Relation to other backends
$n = 1$	homogeneous image-charge self-kernel	bit-identical to `image`
$n = 2$	exact Tagg / Sunde series self-kernel	bit-identical to `image_2layer`
$n \ge 3$	—	`NotImplementedError`; use `mom_sommerfeld` or `fem`

cim is therefore not an independent cross-check of the image-charge family: in a cross-validation table it must be counted as the same computation as image / image_2layer, not as a second engine (ADR-0002 amendment 2026-07-09). The sections below document the complex-image theory and fit_complex_images, which remains available as a standalone spectral helper for a future complete $n \ge 3$ kernel, but is not on any solve path.

Physical context¶

The Complex Image Method (CIM) is the modern workhorse for layered Green's functions in the grounding and antenna literature. It trades a small, controllable approximation error for two important properties:

Closed-form spatial Green's function of the same shape as the homogeneous image-charge sum.
Cost independent of the layer count once the fit is done.

The trick is that the only object in the layered Sommerfeld representation that changes with the layer count is the recursive reflection coefficient $\Gamma_1(\lambda)$ — every other piece is the same exponential / Bessel kernel as in the homogeneous case. If $\Gamma_1(\lambda)$ can be approximated by a finite sum of complex exponentials,

\[ \Gamma_1(\lambda) \;\approx\; \sum_{k=1}^{P} a_k\, e^{-2\lambda \beta_k}, \qquad a_k \in \mathbb{C}, \quad \Re\{\beta_k\} > 0, \]

then substituting that approximation into the Sommerfeld integral gives a closed-form spatial form via the Sommerfeld identity

\[ \int_0^{\infty} e^{-\lambda d}\, J_0(\lambda s)\, d\lambda \;=\; \frac{1}{\sqrt{s^2 + d^2}}, \qquad \Re\{d\} > 0. \]

The result is a $1/r$-type kernel, but with complex image positions $z = -(z_s + 2\beta_k)$. Each pole $(a_k, \beta_k)$ contributes one image; the per-evaluation cost is the same as the homogeneous backend multiplied by $P$, and is independent of $n$.

Governing equation: complex-image kernel¶

Substituting the fit into the layered Sommerfeld integral gives

\[ \varphi(s, z) \;=\; \frac{\rho_1\, I}{4\pi} \Biggl[ \frac{1}{r} + \frac{1}{r_{\text{air}}} + \sum_{k=1}^{P} a_k\, \frac{1}{\sqrt{s^2 + (z + z_s + 2\beta_k)^2}} \Biggr], \]

with $r = \sqrt{s^2 + (z-z_s)^2}$, $r_{\text{air}} = \sqrt{s^2 + (z+z_s)^2}$. The first two terms are the homogeneous direct + air-mirror; the third is the closed-form layered correction. With complex $\beta_k$, the square-root denominator is complex too — but the imaginary parts cancel by symmetry of the fit, so the final potential is real (an implementation would take the real part to suppress numerical residue).

Two caveats, both of which is why this kernel is not wired into a solver today: the form above carries only the $(z + z_s)$ image family, and it lacks the surface-interface multiple-reflection denominator $1 / (1 - \Gamma_1 e^{-2\lambda h_1})$ of the full layered Green's function (see _layered and mom_sommerfeld). Completing it is the prerequisite for letting cim accept $n \ge 3$ again.

Numerical strategy¶

Matrix-pencil fit¶

We approximate $\Gamma_1(\lambda)$ on a uniform sample grid $\lambda_j = \lambda_{\min} + j \Delta$, $j = 0, \dots, N_s - 1$ by $P$ complex exponentials using the matrix-pencil method (Sarkar & Pereira 1995). The procedure:

Split off the asymptote (see below) and sample the decaying remainder $g_j = \Gamma_1(\lambda_j) - K_1$ on the uniform grid.
Form the rectangular Hankel pencil $[Y_0, Y_1]$ with $Y_0[i, k] = g_{i+k}$, $Y_1[i, k] = g_{i+k+1}$, $i = 0, \dots, N_s - L - 1$, $k = 0, \dots, L - 1$, with pencil parameter $L \approx N_s/3$.
SVD $Y_0 = U \Sigma V^*$ and project both blocks onto the $P$-dimensional dominant singular subspace ($U_P, \Sigma_P, V_P$).
The poles $p_k = e^{-2\Delta \beta_k}$ of the sum of exponentials are the eigenvalues of $\Sigma_P^{-1} U_P^* Y_1 V_P^*$.
Recover $\beta_k = -\ln p_k / (2\Delta)$. Discard poles with $|p_k| \ge 1$ (non-decaying / non-physical) or with $\Re\{\beta_k\} \le 0$.
Solve the linear least-squares system for the coefficients $a_k$ on the original samples of $\Gamma_1$, over the surviving set $\{\beta_k\} \cup \{0\}$ — i.e. after the pole filter of step 5, so that the reported weights belong to the reported poles.

The fit is adaptive in the sense that the SVD truncates poles whose singular value falls below $10^{-10}$ relative to the dominant one, so the effective $P$ may be smaller than the requested target. This handles the degenerate case $\Gamma_1 \equiv 0$ (homogeneous soil) gracefully — the fit returns $P = 0$ and the engine collapses to the homogeneous closed form exactly.

The $\lambda \to \infty$ asymptote must be split off¶

A sum of decaying exponentials vanishes for $\lambda \to \infty$, but $\Gamma_1$ does not:

\[ \lim_{\lambda \to \infty} \Gamma_1(\lambda) \;=\; K_1 \;=\; \frac{\rho_2 - \rho_1}{\rho_2 + \rho_1} \;\neq\; 0, \]

because every $e^{-2\lambda h_i}$ in the recursion dies and leaves the top-interface Fresnel coefficient. The constant is carried by an image at $\beta_0 = 0$ — i.e. at the air-mirror position $z = -z_s$, with weight $K_1$, still integrable through the Sommerfeld identity since $d = z + z_s > 0$ — and only $\Gamma_1(\lambda) - K_1$ is handed to the matrix pencil. This follows Li et al. 2006 / Dan et al. 2021.

Up to 0.14.1 groundfield fitted $\Gamma_1$ itself, and the pole filters of step 5 removed precisely the constant term (pole at $p = 1$, $\beta = 0$) that would have carried the asymptote. The residual was therefore of order $|K_1|$ for every stack — e.g. $\rho = [100, 400, 50]\ \Omega\text{m}$, $h = [2, 3]\ \text{m}$: $|K_1| = 0.600$, reported rms $= 0.585$ — and a two-layer stack (where the constant is the whole function) returned $P = 0$ with rms = nan. With the split, the same stacks give rms $\approx 2 \cdot 10^{-13}$ and $0$ exactly (review pass 9, F39).

Failure reporting¶

A fit that cannot represent $\Gamma_1$ no longer passes silently. fit_complex_images reports

rms — residual on the uniform fitting grid, always finite,
rms_extrapolated — residual on an independent 512-point log-spaced verification grid spanning the full requested $\lambda$ range (catches a window that stopped short of the decay band, which the fitting-grid residual cannot see),
converged — both residuals below rms_tol (default $10^{-3}$) and the pencil produced at least one usable pole,
k_inf — the asymptote $K_1$ carried by the $\beta = 0$ image,

and emits a ComplexImageFitWarning whenever converged is False.

Sample grid choice¶

Only the thicknesses $h_2, \dots, h_{n-1}$ enter $\Gamma_1$ — the recursion never uses $h_1$. The grid is uniform on $\lambda \in [\lambda_{\min}, \lambda_{\max}]$ with

$\lambda_{\min} = 10^{-3} / h_{\max}$, $h_{\max} = \max_{i \ge 2} h_i$ (well below the band where $\Gamma_1$ has structure),
$\Delta \le 1 / (8\, h_{\max})$, so the slowest exponential (scale $1/(2 h_{\max})$) is resolved by at least eight samples,
$\lambda_{\max} = 50 / h_{\min}$, $h_{\min} = \min_{i \ge 2} h_i$ (above the decay scale of the thinnest interface); with $\Delta$ capped as above the grid may stop earlier, which is exactly what rms_extrapolated monitors.

Up to 0.14.1 both bounds were derived from $h_{\min}$, so for a stack like $h = [0.5, 10]\ \text{m}$ the step was $\Delta = 1.59$ against a structure scale of $0.05$ — the grid stepped straight over the transition and the pencil collapsed to $P = 1$.

This is a single-segment fit. Dan et al. 2021 propose a segmented-sampling variant in which the $\lambda$-axis is split into pieces and a separate fit is run on each. For the typical contrast and layer-count range a single-segment fit with $P = 8$ is sufficient; the segmented variant becomes attractive when $n \gtrsim 5$ or contrasts are extreme.

Self-action strategy¶

For $n = 1$ the self-action is the homogeneous self-kernel ($\Gamma_1 \equiv 0$, no extra images).

For $n = 2$ the engine uses the closed-form Tagg / Sunde self-kernel of image_2layer, with allow_cross_layer=True so that a rod crossing the interface dispatches to the rigorous cross-layer path (ADR-0007). Here $\Gamma_1 \equiv K_1$ is constant in $\lambda$, so the geometric image series is the complex-image representation (the single $\beta = 0$ image) — evaluating it in closed form is exact and cheaper than any fit.

For $n \ge 3$ the engine raises NotImplementedError. No complex-image fit is evaluated on any solve path. Up to 0.14.1 the fit was computed on every solve and its diagnostics were published as cim_n_images / cim_rms although neither reachable path consumed them; readers were invited to judge "how good is the complex-image approximation here" from a NaN produced by a fit with no influence on the answer (review pass 9, F34). The metadata now carries

Key	Value	Meaning
`cim_fit_used`	`False`	no complex-image fit entered this result
`cim_n_images`	`0`	consequently, no images
`cim_rms`	`None`	no fit, hence no residual (was `nan`)
`reduces_to`	`"image"` / `"image_2layer"`	backend this result is bit-identical to

Reaction matrix¶

The closed-form self-kernel handles direct + air-mirror + layered correction in one call. Cluster constraints are enforced through the same multi-port system as in image.

Validity envelope¶

Property	Range / value
Soil model	`HomogeneousSoil`, `TwoLayerSoil` (a `MultiLayerSoil` is accepted only while it reduces to $n \le 2$)
Frequency	quasi-static, $f < 1\,\text{kHz}$
Electrode placement	free — $n = 2$ interface crossings dispatch to the ADR-0007 cross-layer kernel
Number of layers	$n \le 2$; $n \ge 3$ raises `NotImplementedError`
Number of images $P$	not applicable — no fit runs (`fit_complex_images` itself: 4–12 typical, 8 default)
Layer-contrast range	as `image_2layer` (exact series, truncation controlled by `Engine.image_max_terms` / `image_series_tol`)

Convergence and cost¶

Accuracy. Identical to image ($n = 1$) / image_2layer ($n = 2$) by construction — the same kernels are called.
Fit accuracy (of the standalone fit_complex_images, not of the engine). With $P = 8$ decaying images on a 64-sample grid plus the analytic $\beta = 0$ asymptote image, the $\Gamma_1$ residual measured on three- and four-layer stacks with $|K_i|$ up to $0.9$ is $\le 5 \cdot 10^{-8}$ on the fitting grid and $\le 3 \cdot 10^{-4}$ on the wide verification grid. Anything worse raises a ComplexImageFitWarning.
Per-segment cost. $O(N^2)$ kernel evaluations plus the multi-port solve — the same as image / image_2layer.
Reduction. At $n = 1$ bit-exact match with image; at $n = 2$ bit-exact match with image_2layer.

Cross-validation notes¶

Counterpart	Expected agreement	What is checked
`image` ($n = 1$)	bit-exact	same computation — the homogeneous self-kernel; no independent information
`image_2layer` ($n = 2$)	bit-exact	same computation — the Tagg / Sunde kernel; no independent information
`mom` / `bem` (any $n \le 2$)	$\le 5\,\%$	uniform-current vs. Galerkin / collocation on the same kernel
`mom_sommerfeld` ($n = 2$)	$\le 5\,\%$	direct quadrature of the full layered Green's function — the only genuinely independent kernel
`fem` ($n \le 2$)	$\le 10\,\%$	volume PDE, independent problem form (equivalent-hemisphere bias)

cim is not an independent cross-validation engine. For every soil it accepts it evaluates the same closed-form kernel as the image-charge family, so a table that lists image_2layer and cim as two agreeing engines is reporting one computation twice. Real independence comes from mom_sommerfeld (quadrature) and fem (volume PDE); see the ADR-0002 amendment.

References¶

Sarkar, T. K. & Pereira, O. (1995). Using the matrix pencil method to estimate the parameters of a sum of complex exponentials. IEEE Antennas & Propagation Magazine 37(1). The primary reference for the fit algorithm.
Li, Z.-X., Chen, W., Fan, J.-B. & Lu, J. (2006). A novel mathematical modeling of grounding system buried in multilayer earth. IEEE PWRD 21(3). Quasi-static CIM (QCIM) for grounding systems; the closed-form Green's function in the form used here.
Dan, Y. et al. (2021). Segmented sampling least squares algorithm for Green's function of arbitrary layered soil. IEEE PWRD 36(3). The segmented-sampling refinement and a systematic study of the accuracy versus $P$.
Hua, Y. & Sarkar, T. K. (1990). Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise. IEEE Trans. ASSP 38(5). The original matrix-pencil derivation.

Example¶

A two-layer solve — the deepest stack cim accepts:

import groundfield as gf

soil = gf.TwoLayerSoil(rho_1=100.0, rho_2=400.0, h_1=2.0)
world = gf.create_world(soil=soil)
gf.create_electrode(world, "rod", name="g1",
                    position=(0.0, 0.0, 0.0), length=3.0)
gf.create_source(world, attached_to="g1", magnitude=1.0)

engine = gf.create_engine(backend="cim",
                          segment_length=0.1,
                          frequencies=[50.0])
result = world.solve(engine)
print(result.cluster_impedance("g1")[0])   # (67.6864+0j) Ohm
print(result.metadata["cim_fit_used"])     # False — no fit runs
print(result.metadata["reduces_to"])       # 'image_2layer'

Three or more layers raise, and point at the engines that do implement the full layered Green's function:

soil = gf.MultiLayerSoil(layers=[
    gf.SoilLayer(resistivity=80.0, thickness=0.5),
    gf.SoilLayer(resistivity=300.0, thickness=2.0),
    gf.SoilLayer(resistivity=50.0),  # semi-infinite bottom
])
world = gf.create_world(soil=soil)
gf.create_electrode(world, "rod", name="g1",
                    position=(0.0, 0.0, 0.0), length=0.4)
gf.create_source(world, attached_to="g1", magnitude=1.0)

# gf.create_engine(backend="cim", ...).solve(world)
#   -> NotImplementedError: cim: n_layers = 3 >= 3 is not supported ...
engine = gf.create_engine(backend="mom_sommerfeld", segment_length=0.1)
result = world.solve(engine)
print(result.cluster_impedance("g1")[0])

The spectral helper can still be called directly — for instance to check how well a stack's $\Gamma_1$ would be represented:

import numpy as np
from groundfield.solver._layered import LayerStack
from groundfield.solver.cim import fit_complex_images

stack = LayerStack(rhos=np.array([100.0, 400.0, 50.0]),
                   h=np.array([2.0, 3.0]))
fit = fit_complex_images(stack)
print(fit.a.size, fit.k_inf, fit.rms, fit.converged)
# 8 0.6 2.1e-13 True   (index 0 is the beta = 0 asymptote image)

API reference¶

cim ¶

Complex Image Method backend (cim).

Mathematical / physical model

The :mod:image_nlayer backend expands the upward-looking reflection $\Gamma_1(\lambda)$ as a power series in the per-layer attenuation factors $e^{-2\lambda h_i}$ — a representation that converges fast for two layers but generates many terms for $n \ge 3$ and high contrasts.

The Complex Image Method (CIM) instead approximates $\Gamma_1(\lambda)$ directly by a finite sum of complex exponentials, $$ \Gamma_1(\lambda) \;\approx\; \sum_{k=1}^{P} a_k\, e^{-2\lambda \beta_k}, $$ with complex coefficients $a_k \in \mathbb{C}$ and complex "image depths" $\beta_k \in \mathbb{C}$ (with $\Re\{\beta_k\} > 0$ to keep the integrals convergent). Substituting this approximation into the Sommerfeld integral and using the closed-form $$ \int_0^{\infty} e^{-\lambda d} J_0(\lambda s)\, d\lambda \;=\; \frac{1}{\sqrt{s^2 + d^2}} $$ immediately gives a closed-form spatial Green's function of the same shape as the homogeneous image-charge sum, but with complex image positions. A kernel built from such a fit would read $$ \varphi(s, z) \;=\; \frac{\rho_1\, I}{4\pi}\, \Bigl(\frac{1}{r} + \frac{1}{r_{\text{air}}} + \sum_{k=1}^{P} a_k \Bigl(\frac{1}{\sqrt{s^2 + (z + z_s + 2\beta_k)^2}} + \frac{1}{\sqrt{s^2 + (z - z_s + 2\beta_k)^2}}\Bigr) \Bigr). $$ The cost of one potential evaluation would therefore be the same as the homogeneous backend, multiplied by 2 * P — independent of the layer count. No backend evaluates this kernel today; see "Solver status of the fit" below, and note that the form above still needs the $2 h_1$ image families and the multiple-reflection denominator of :mod:groundfield.solver._layered before it represents the layered Green's function completely.

Numerical fit strategy

We use the matrix-pencil method (a numerically stable variant of Prony's algorithm) to fit a sample of $\Gamma_1(\lambda)$ on a uniform grid in $\lambda$ to P complex exponentials. This is a faithful Python re-implementation of the segmented-sampling least-squares idea of Dan et al. 2021 (without the segmentation heuristic, which is needed mainly for very many layers; for the typical two-/three-layer use cases a single segment with a moderately oversampled grid is enough).

One structural detail matters: a pure sum of decaying exponentials vanishes for $\lambda \to \infty$, whereas $$ \lim_{\lambda \to \infty} \Gamma_1(\lambda) \;=\; K_1 \;=\; \frac{\rho_2 - \rho_1}{\rho_2 + \rho_1} \;\neq\; 0 $$ (every $e^{-2 \lambda h_i}$ in the recursion dies, leaving the top-interface Fresnel coefficient). The asymptote is therefore split off analytically — it is carried by the single $\beta = 0$ image, i.e. an image at the air-mirror position $z = -z_s$ with weight $K_1$ — and only the decaying remainder $\Gamma_1(\lambda) - K_1$ is handed to the matrix pencil. Fitting $\Gamma_1$ itself (as this module did up to 0.14.1) leaves a residual of order $|K_1|$ for every stack, because the pole filters discard precisely the constant term that would carry the asymptote (review pass 9, F39).

Solver status of the fit (0.15.0)

:func:fit_complex_images is a standalone spectral helper: no solver path evaluates it any more. cim (and bem, which shares these kernels) rejects $n \ge 3$ soils, and for $n \le 2$ both engines use exact closed-form self-kernels — the homogeneous image-charge sum for $n = 1$ and the Tagg/Sunde series of :mod:groundfield.solver.image_2layer for $n = 2$ — which do not involve a fit at all. Up to 0.14.1 the fit was nevertheless computed on every solve and its (failed) diagnostics were reported as metadata['cim_n_images'] / ['cim_rms'], which invited readers to judge an approximation that had no influence on the answer (review pass 9, F34). The call is gone; the metadata now carries cim_fit_used: False and a reduces_to key naming the backend this result is bit-identical to.

Validity

Quasi-static, $f < 1\,\mathrm{kHz}$.
n_layers == 1 → homogeneous image-charge sum, bit-identical to image ($\Gamma_1 \equiv 0$).
n_layers == 2 → exact Tagg/Sunde series, bit-identical to image_2layer. cim is therefore not an independent cross-check of image_2layer in this regime (ADR-0002 amendment 2026-07-09).
n_layers ≥ 3 is rejected (audit 2026-07-08, WP-E): the historic complex-image expansion approximated an incomplete Green's function (single image family, no multiple reflections — see :mod:groundfield.solver._layered) and systematically underestimated the layered correction. Use mom_sommerfeld or fem for three and more layers.

References

Sarkar, T. K. & Pereira, O. (1995). Using the matrix pencil method to estimate the parameters of a sum of complex exponentials, IEEE Antennas & Propagation Magazine 37(1).
Li, Z.-X. et al. (2006). A novel mathematical modeling of grounding system buried in multilayer earth, IEEE PWRD 21(3).
Dan, Y. et al. (2021). Segmented sampling least squares algorithm for Green's function of arbitrary layered soil, IEEE PWRD 36(3).

ComplexImageFit `dataclass` ¶

ComplexImageFit(
    a: np.ndarray,
    beta: np.ndarray,
    rms: float,
    k_inf: float = 0.0,
    rms_extrapolated: float = 0.0,
    converged: bool = True,
)

Result of :func:fit_complex_images.

The fit represents $\Gamma_1(\lambda) \approx \sum_k a_k e^{-2 \lambda \beta_k}$. Index 0 is the analytic asymptote image ($\beta_0 = 0$, $a_0 \approx K_1$) whenever the stack has one; the remaining entries are the decaying images recovered by the matrix pencil.

Attributes:

Name	Type	Description
`a`	`(ndarray, shape(P), complex)`	Image weights $a_k$.
`beta`	`(ndarray, shape(P), complex)`	Image depths $\beta_k$ (units of metres).
`rms`	`float`	RMS of the residual $\|\text{fit} - \Gamma_1\|$ on the sample grid. Always finite: a failed fit reports the residual it actually achieved (and warns) instead of `nan`.
`k_inf`	`float`	The $\lambda \to \infty$ asymptote $K_1 = (\rho_2 - \rho_1)/(\rho_2 + \rho_1)$ that the $\beta = 0$ image carries. `0.0` for a homogeneous stack.
`rms_extrapolated`	`float`	RMS of the same residual on an independent log-spaced verification grid covering the full requested $\lambda$ range. Larger than `rms` when the uniform fitting window did not reach the decay band of the thinnest layer (large thickness spread) or when the pencil overfitted.
`converged`	`bool`	Both residuals are `<= rms_tol` and the decaying remainder was either numerically flat or successfully fitted. `False` means a :class:`ComplexImageFitWarning` was emitted.

ComplexImageFitWarning ¶

Bases: UserWarning

The complex-image fit of $\Gamma_1(\lambda)$ is unreliable.

Emitted by :func:fit_complex_images when the RMS residual on the sample grid or on the independent verification grid exceeds rms_tol, or when the matrix pencil returned no usable pole for a remainder that is not numerically flat. All three situations mean the returned images do not represent $\Gamma_1$ and must not be fed into a Green's function.

Up to 0.14.1 the same situations were silent: the pole filters (|p| < 0.999, Re β > 0) discarded the constant term that carries the $\lambda \to \infty$ asymptote $K_1$, so the residual was of order $|K_1|$ for every stack, and a two-layer stack (where that constant is the whole function) returned P = 0 with rms = nan (review pass 9, F39).

A dedicated category lets callers opt in / out with warnings.simplefilter("error", ComplexImageFitWarning).

fit_complex_images ¶

fit_complex_images(
    stack: LayerStack,
    *,
    n_images: int = 8,
    n_samples: int = 64,
    lambda_min_factor: float = 0.001,
    lambda_max_factor: float = 50.0,
    rms_tol: float = _DEFAULT_RMS_TOL,
    warn: bool = True
) -> ComplexImageFit

Fit Γ_1(λ) of an n-layer stack by n_images complex exponentials.

Physical / mathematical context

The upward-looking reflection coefficient of a stratified half-space (see :func:groundfield.solver._layered.reflection_gamma) is a rational function of the interface attenuations $e^{-2 \lambda h_i}$. It is not a decaying function: it tends to the top-interface Fresnel coefficient

\[ \lim_{\lambda \to \infty} \Gamma_1(\lambda) \;=\; K_1 \;=\; \frac{\rho_2 - \rho_1}{\rho_2 + \rho_1}, \]

while a sum $\sum_k a_k e^{-2\lambda\beta_k}$ with $\Re\{\beta_k\} > 0$ tends to 0. The asymptote is hence split off analytically and carried by an image at $\beta_0 = 0$ (weight $K_1$; in the spatial domain an image at the air-mirror position $z = -z_s$, integrable through the Sommerfeld identity because $d = z + z_s > 0$), and the matrix pencil only sees the genuinely decaying remainder $\Gamma_1(\lambda) - K_1$.

Steps:

Build a uniform sampling grid $\lambda_j = \lambda_{\min} + j \Delta$. Only the thicknesses $h_2, \dots, h_{n-1}$ enter $\Gamma_1$ (the recursion never uses $h_1$), so the structure scale is $h_{\max} = \max_i h_{i\ge2}$ and the decay scale is $h_{\min} = \min_i h_{i\ge2}$: the step is capped at $\Delta \le 1/(8 h_{\max})$ so the slowest exponential is resolved, and the window reaches up to $\lambda_{\max} = \lambda_{\max,\text{factor}}/h_{\min}$ where the fastest one has died.
Sample $g_j = \Gamma_1(\lambda_j)$ and subtract $K_1$.
Apply the matrix-pencil method to the remainder to recover the poles $p_k = e^{-2 \Delta \beta_k}$, filtered to $|p_k| < 0.999$ and $\Re\{\beta_k\} > 0$.
Re-solve the least-squares coefficients $a_k$ for the surviving $\{\beta_k\} \cup \{0\}$ on the original samples of $\Gamma_1$ — the filtered model is therefore self-consistent (up to 0.14.1 the coefficients were solved before the filter and then merely sliced, so the reported weights belonged to a different model).

Special cases:

stack.n_layers <= 1 → empty fit ($\Gamma_1 \equiv 0$).
stack.n_layers == 2 → $\Gamma_1 \equiv K_1$ exactly; the single $\beta = 0$ image reproduces it with rms = 0 (up to 0.14.1 this returned P = 0 and rms = nan).
$\rho_1 = \dots = \rho_n$ → empty fit.

Notes

No solver path calls this function: cim/bem reject $n \ge 3$ and use exact closed-form kernels for $n \le 2$ (see the module docstring). It is kept public as the spectral building block for a future complete $n \ge 3$ complex-image kernel, and it now reports failure instead of hiding it.

Parameters:

Name	Type	Description	Default
`stack`	`LayerStack`	Layer stack to fit.	required
`n_images`	`int`	Target number of decaying complex images. Sensible range 4–12; the returned $P$ is one larger when the asymptote image is present, and may be smaller when the SVD truncates poles.	`8`
`n_samples`	`int`	Number of samples drawn from Γ_1(λ).	`64`
`lambda_min_factor`	`float`	Lower grid bound as a multiple of $1/h_{\max}$.	`0.001`
`lambda_max_factor`	`float`	Upper grid bound as a multiple of $1/h_{\min}$ (an upper bound only — the step cap of item 1 above may stop the grid earlier).	`50.0`
`rms_tol`	`float`	Residual (on the fitting and on the verification grid) above which the fit counts as failed (`converged=False` plus a :class:`ComplexImageFitWarning`).	`_DEFAULT_RMS_TOL`
`warn`	`bool`	Set to `False` to suppress the warning in bulk sweeps; `converged` / `rms` / `rms_extrapolated` still report the failure.	`True`

Returns:

Type	Description
`ComplexImageFit`

Warns:

Type	Description
`ComplexImageFitWarning`	If `rms` or `rms_extrapolated` exceeds `rms_tol`, or the pencil produced no usable pole for a non-flat remainder.

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

def fit_complex_images(
    stack: LayerStack,
    *,
    n_images: int = 8,
    n_samples: int = 64,
    lambda_min_factor: float = 1e-3,
    lambda_max_factor: float = 50.0,
    rms_tol: float = _DEFAULT_RMS_TOL,
    warn: bool = True,
) -> ComplexImageFit:
    """Fit Γ_1(λ) of an n-layer stack by ``n_images`` complex exponentials.

    Physical / mathematical context
    ------------------------------
    The upward-looking reflection coefficient of a stratified
    half-space (see :func:`groundfield.solver._layered.reflection_gamma`)
    is a rational function of the interface attenuations
    $e^{-2 \\lambda h_i}$. It is **not** a decaying function: it
    tends to the top-interface Fresnel coefficient

    $$
    \\lim_{\\lambda \\to \\infty} \\Gamma_1(\\lambda)
    \\;=\\; K_1 \\;=\\; \\frac{\\rho_2 - \\rho_1}{\\rho_2 + \\rho_1},
    $$

    while a sum $\\sum_k a_k e^{-2\\lambda\\beta_k}$ with
    $\\Re\\{\\beta_k\\} > 0$ tends to 0. The asymptote is hence
    split off analytically and carried by an image at
    $\\beta_0 = 0$ (weight $K_1$; in the spatial domain an image
    at the air-mirror position $z = -z_s$, integrable through the
    Sommerfeld identity because $d = z + z_s > 0$), and the matrix
    pencil only sees the genuinely decaying remainder
    $\\Gamma_1(\\lambda) - K_1$.

    Steps:

    1. Build a uniform sampling grid
       $\\lambda_j = \\lambda_{\\min} + j \\Delta$. Only the
       thicknesses $h_2, \\dots, h_{n-1}$ enter $\\Gamma_1$
       (the recursion never uses $h_1$), so the *structure* scale is
       $h_{\\max} = \\max_i h_{i\\ge2}$ and the *decay* scale is
       $h_{\\min} = \\min_i h_{i\\ge2}$: the step is capped at
       $\\Delta \\le 1/(8 h_{\\max})$ so the slowest exponential is
       resolved, and the window reaches up to
       $\\lambda_{\\max} = \\lambda_{\\max,\\text{factor}}/h_{\\min}$
       where the fastest one has died.
    2. Sample $g_j = \\Gamma_1(\\lambda_j)$ and subtract $K_1$.
    3. Apply the matrix-pencil method to the remainder to recover the
       poles $p_k = e^{-2 \\Delta \\beta_k}$, filtered to
       $|p_k| < 0.999$ and $\\Re\\{\\beta_k\\} > 0$.
    4. Re-solve the least-squares coefficients $a_k$ for the
       *surviving* $\\{\\beta_k\\} \\cup \\{0\\}$ on the original
       samples of $\\Gamma_1$ — the filtered model is therefore
       self-consistent (up to 0.14.1 the coefficients were solved
       before the filter and then merely sliced, so the reported
       weights belonged to a different model).

    Special cases:

    - ``stack.n_layers <= 1`` → empty fit ($\\Gamma_1 \\equiv 0$).
    - ``stack.n_layers == 2`` → $\\Gamma_1 \\equiv K_1$ exactly;
      the single $\\beta = 0$ image reproduces it with ``rms = 0``
      (up to 0.14.1 this returned ``P = 0`` and ``rms = nan``).
    - $\\rho_1 = \\dots = \\rho_n$ → empty fit.

    Notes
    -----
    No solver path calls this function: ``cim``/``bem`` reject
    $n \\ge 3$ and use exact closed-form kernels for $n \\le 2$
    (see the module docstring). It is kept public as the spectral
    building block for a future complete $n \\ge 3$ complex-image
    kernel, and it now reports failure instead of hiding it.

    Parameters
    ----------
    stack
        Layer stack to fit.
    n_images
        Target number of *decaying* complex images. Sensible range
        4–12; the returned $P$ is one larger when the asymptote
        image is present, and may be smaller when the SVD truncates
        poles.
    n_samples
        Number of samples drawn from Γ_1(λ).
    lambda_min_factor
        Lower grid bound as a multiple of $1/h_{\\max}$.
    lambda_max_factor
        Upper grid bound as a multiple of $1/h_{\\min}$ (an upper
        bound only — the step cap of item 1 above may stop the grid
        earlier).
    rms_tol
        Residual (on the fitting *and* on the verification grid) above
        which the fit counts as failed (``converged=False`` plus a
        :class:`ComplexImageFitWarning`).
    warn
        Set to ``False`` to suppress the warning in bulk sweeps;
        ``converged`` / ``rms`` / ``rms_extrapolated`` still report the
        failure.

    Returns
    -------
    ComplexImageFit

    Warns
    -----
    ComplexImageFitWarning
        If ``rms`` or ``rms_extrapolated`` exceeds ``rms_tol``, or the
        pencil produced no usable pole for a non-flat remainder.
    """
    if stack.n_layers <= 1:
        return _empty_fit()

    k_inf = float(stack.K[0])

    # Only h_2 … h_{n-1} enter Γ_1 — the recursion in
    # `reflection_gamma` never touches h_1.
    h_rel = np.asarray(stack.h[1:], dtype=float)

    if h_rel.size == 0:
        # n == 2: Γ_1 ≡ K_1, a constant. Exactly one image, at β = 0.
        if abs(k_inf) < _FLAT_TOL:
            return _empty_fit()
        return ComplexImageFit(
            a=np.array([k_inf + 0j]),
            beta=np.zeros(1, dtype=complex),
            rms=0.0,
            k_inf=k_inf,
            converged=True,
        )

    h_struct = float(np.max(h_rel))
    h_decay = float(np.min(h_rel))
    lam_min = lambda_min_factor / h_struct
    lam_max = lambda_max_factor / h_decay

    # Uniform spacing required by the matrix-pencil method; capped so
    # the slowest exponential (scale 1/(2 h_struct)) is resolved.
    delta = min(
        (lam_max - lam_min) / (n_samples - 1),
        1.0 / (_SAMPLES_PER_STRUCTURE_SCALE * h_struct),
    )
    lam = lam_min + delta * np.arange(n_samples)
    g = reflection_gamma(stack, lam).astype(complex)

    # Degenerate case ρ_1 = ρ_2 = … = ρ_n: nothing to fit.
    if np.max(np.abs(g)) < _FLAT_TOL:
        return _empty_fit()

    resid = g - k_inf
    if np.max(np.abs(resid)) < _FLAT_TOL:
        # Γ_1 is constant on the whole window (e.g. ρ_2 = … = ρ_n).
        return ComplexImageFit(
            a=np.array([k_inf + 0j]),
            beta=np.zeros(1, dtype=complex),
            rms=float(np.sqrt(np.mean(np.abs(resid) ** 2))),
            k_inf=k_inf,
            converged=True,
        )

    beta_decay = _matrix_pencil_beta(resid, delta, n_images)
    pencil_failed = beta_decay.size == 0

    # β = 0 (the analytic asymptote) plus the decaying images.
    beta = np.concatenate([np.zeros(1, dtype=complex), beta_decay])

    # Coefficients a_k by linear least squares on the original samples
    # of Γ_1 — solved *after* all pole filtering.
    A = np.exp(-2.0 * lam[:, None] * beta[None, :])  # (n_samples, P)
    a, *_ = np.linalg.lstsq(A, g, rcond=None)
    rms = float(np.sqrt(np.mean(np.abs(A @ a - g) ** 2)))

    # Independent verification on a log-spaced grid over the *full*
    # requested λ range: the uniform fitting window may stop well
    # below λ_max when the thickness spread is large, and a residual
    # measured on the fitting grid alone cannot see that.
    lam_v = np.logspace(np.log10(lam_min), np.log10(lam_max), _N_VERIFY)
    g_v = reflection_gamma(stack, lam_v).astype(complex)
    model_v = np.exp(-2.0 * lam_v[:, None] * beta[None, :]) @ a
    rms_wide = float(np.sqrt(np.mean(np.abs(model_v - g_v) ** 2)))

    converged = (
        bool(rms <= rms_tol)
        and bool(rms_wide <= rms_tol)
        and not pencil_failed
    )

    _log.info(
        "cim.fit: n_images_used=%d (P_decay=%d), rms=%.2e, "
        "rms_extrapolated=%.2e, K_1=%.4f, n_layers=%d, "
        "lam=[%.3g, %.3g], delta=%.3g",
        beta.size, beta_decay.size, rms, rms_wide, k_inf,
        stack.n_layers, lam[0], lam[-1], delta,
    )
    if not converged and warn:
        if pencil_failed:
            reason = (
                "the matrix pencil returned no usable pole for a "
                "non-flat remainder Γ_1 - K_1"
            )
        elif rms > rms_tol:
            reason = f"RMS residual {rms:.3e} > rms_tol {rms_tol:.3e}"
        else:
            reason = (
                f"RMS residual on the verification grid "
                f"{rms_wide:.3e} > rms_tol {rms_tol:.3e}"
            )
        warnings.warn(
            f"fit_complex_images: unreliable fit for an "
            f"{stack.n_layers}-layer stack ({reason}). Grid: "
            f"lambda in [{lam[0]:.3g}, {lam[-1]:.3g}], delta="
            f"{delta:.3g}, structure scale 1/(2*h_max)="
            f"{1.0 / (2.0 * h_struct):.3g}, decay scale "
            f"1/(2*h_min)={1.0 / (2.0 * h_decay):.3g}. Raise "
            "n_samples / n_images, or widen lambda_max_factor when the "
            "thickness spread h_max/h_min is large. Do not use these "
            "images in a Green's function.",
            ComplexImageFitWarning,
            stacklevel=2,
        )
    return ComplexImageFit(
        a=a,
        beta=beta,
        rms=rms,
        k_inf=k_inf,
        rms_extrapolated=rms_wide,
        converged=converged,
    )

solve_cim ¶

solve_cim(world: 'World', engine: 'Engine') -> FieldResult

Complex-Image-Method solver for layered soil.

Accepts :class:HomogeneousSoil and :class:TwoLayerSoil; :class:MultiLayerSoil is accepted only while it reduces to $n \le 2$ (see below).

What actually runs

$n = 1$ → homogeneous image-charge self-kernel. Bit-identical to image.
$n = 2$ → exact Tagg/Sunde series self-kernel. Bit-identical to image_2layer (which is why cim is not an independent cross-check of it — ADR-0002 amendment 2026-07-09).
$n \ge 3$ → :class:NotImplementedError (audit 2026-07-08, WP-E). Use mom_sommerfeld or fem.

No complex-image fit is evaluated. Up to 0.14.1 this solver called :func:fit_complex_images on every solve and published its diagnostics even though neither reachable path consumed them (review pass 9, F34); the metadata now says so explicitly with cim_fit_used = False instead of reporting a NaN residual.