`cim` — Complex Image Method¶

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 (the implementation takes the real part to suppress numerical residue).

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:

Sample $g_j = \Gamma_1(\lambda_j)$ 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.

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.

Sample grid choice¶

The grid is uniform on $\lambda \in [\lambda_{\min}, \lambda_{\max}]$ with bounds expressed as multiples of $1 / h_{\min}$:

$\lambda_{\min} = 10^{-3} / h_{\min}$ (well below the band where $\Gamma_1$ has structure),
$\lambda_{\max} = 50 / h_{\min}$ (well above the decay scale of every interface).

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 (sub-percent residual); the segmented variant becomes attractive when $n \gtrsim 5$ or contrasts are extreme.

Self-action strategy¶

For $n = 1$ the fit returns $P = 0$ and the self-action falls back exactly on the homogeneous self-kernel.

For $n = 2$ the engine deliberately bypasses the matrix-pencil fit and reuses the closed-form Tagg / Sunde self-kernel of image_2layer. The reason is numerical conditioning: $\Gamma_1 \equiv K_1$ is a constant in $\lambda$, so the matrix pencil sees zero variation in the Hankel block and would have to fit a single pole on the unit circle. This is ill-conditioned and would produce $P = 1$ with a near-singular pole at $\beta = 0$. Falling back to the geometric series gives a bit-exact match at zero numerical risk.

For $n \ge 3$ the genuine matrix-pencil fit runs. The result metadata exposes cim_n_images (the effective $P$) and cim_rms (the residual on the sample grid) for diagnostics.

Reaction matrix¶

After the fit, the homogeneous self-kernel handles the direct + air-mirror part; the layered correction is a sum over the $P$ complex images, evaluated with point-source kernels at every (field, source) segment pair. Cluster constraints are enforced through the same multi-port system as in image.

Validity envelope¶

Property	Range / value
Soil model	`HomogeneousSoil`, `TwoLayerSoil`, or `MultiLayerSoil`
Frequency	quasi-static, $f < 1\,\text{kHz}$
Electrode placement	every segment must lie inside the upper layer
Number of images $P$	4–12 typical; 8 default
Layer-contrast range	$
Number of layers	tested up to $n = 8$ (Dan et al. report)

Convergence and cost¶

Fit accuracy. With $P = 8$ on a 64-sample grid, the $\Gamma_1$ residual stays below $10^{-3}$ for typical contrasts. The cluster-impedance error in the final result is typically smaller than the fit residual by an order of magnitude (the spatial integration averages out the spectral residual).
Per-segment cost. $O(N^2 \cdot P)$ kernel evaluations, plus $O(N_s P^2)$ for the fit (one-shot, dominated by the SVD). Total cost is a small constant multiple of the homogeneous engine, and is independent of $n$.
Reduction. At $n = 1$ the fit gives $P = 0$ and the engine collapses to image exactly. At $n = 2$ the engine reuses the exact Tagg / Sunde kernel — bit-exact match with image_2layer.

Cross-validation notes¶

Counterpart	Expected agreement	What is checked
`image` ($n = 1$)	bit-exact	fit returns $P = 0$
`image_2layer` ($n = 2$)	bit-exact	engine reuses the Tagg / Sunde kernel
`mom_sommerfeld` ($n \ge 3$)	$\le 5\,\%$	quadrature is the reference
`bem` (any $n$)	$\le 5\,\%$	shares the CIM kernel for $n \ge 3$
Three-layer collapse $\rho_1 = \rho_2 = \rho_3$	$\le 2\,\%$	engine reduces to the 1-layer case via the empty-fit shortcut

The engine's cross-validation role is the fast layered solver: it covers the $n \ge 3$ regime where image_nlayer refuses, and is checked against mom_sommerfeld for absolute correctness.

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¶

import groundfield as gf

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)

engine = gf.create_engine(backend="cim",
                          segment_length=0.1,
                          frequencies=[50.0])
result = world.solve(engine)
print(result.cluster_impedance("g1")[0])
print(result.metadata.get("cim_n_images"))  # effective P

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. The kernel of every backend that consumes the layered Green's function (this engine and bem) becomes: $$ \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 is therefore the same as the homogeneous backend, multiplied by 2 * P — independent of the layer count.

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 log-spaced 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).

Validity

Quasi-static, $f < 1\,\mathrm{kHz}$.
For n_layers == 1 the fit returns P = 0 and the engine collapses to the homogeneous image-charge sum exactly.
For n_layers == 2 with a moderate $P$ (typically 6–10) the engine reproduces image_2layer to better than 0.1 %.
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
)

Result of :func:fit_complex_images.

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 on the sample grid (diagnostic only).

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

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

Steps:

Build a uniform sampling grid $\lambda_j = j \Delta$ on $(0, \lambda_{\max}]$. The grid step $\Delta$ is chosen from the smallest layer thickness $h_{\min}$ so that $\Delta = (1/h_{\min}) \cdot (\lambda_{\max,\text{factor}} / N)$ covers the regime where Γ_1 has any structure (it decays for $\lambda \gg 1/h_{\min}$).
Sample $g_j = \Gamma_1(\lambda_j)$.
Apply the matrix-pencil method to g_j to recover the $P$ poles $p_k = e^{-2 \Delta \beta_k}$ of the sum of exponentials.
Solve the least-squares system for the coefficients $a_k$.

For stack.n_layers <= 1 returns an empty fit (Γ_1 ≡ 0).

Parameters:

Name	Type	Description	Default
`stack`	`LayerStack`	Layer stack to fit.	required
`n_images`	`int`	Target number of complex images $P$. Sensible range 4–12 for two-/three-layer stacks; small numbers are fine for soft contrasts, harder contrasts need more.	`8`
`n_samples`	`int`	Number of samples drawn from Γ_1(λ).	`64`
`lambda_min_factor`	`float`	Bounds of the λ-grid expressed as multiples of $1/h_{\min}$.	`0.001`
`lambda_max_factor`	`float`	Bounds of the λ-grid expressed as multiples of $1/h_{\min}$.	`0.001`

Returns:

Type	Description
`ComplexImageFit`

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,
) -> ComplexImageFit:
    """Fit Γ_1(λ) of an n-layer stack by ``n_images`` complex exponentials.

    Steps:

    1. Build a uniform sampling grid $\\lambda_j = j \\Delta$
       on $(0, \\lambda_{\\max}]$. The grid step
       $\\Delta$ is chosen from the smallest layer thickness
       $h_{\\min}$ so that
       $\\Delta = (1/h_{\\min}) \\cdot
       (\\lambda_{\\max,\\text{factor}} / N)$ covers the regime where
       Γ_1 has any structure (it decays for
       $\\lambda \\gg 1/h_{\\min}$).
    2. Sample $g_j = \\Gamma_1(\\lambda_j)$.
    3. Apply the matrix-pencil method to ``g_j`` to recover the
       $P$ poles $p_k = e^{-2 \\Delta \\beta_k}$ of the
       sum of exponentials.
    4. Solve the least-squares system for the coefficients $a_k$.

    For ``stack.n_layers <= 1`` returns an empty fit (Γ_1 ≡ 0).

    Parameters
    ----------
    stack
        Layer stack to fit.
    n_images
        Target number of complex images $P$. Sensible range
        4–12 for two-/three-layer stacks; small numbers are fine for
        soft contrasts, harder contrasts need more.
    n_samples
        Number of samples drawn from Γ_1(λ).
    lambda_min_factor, lambda_max_factor
        Bounds of the λ-grid expressed as multiples of
        $1/h_{\\min}$.

    Returns
    -------
    ComplexImageFit
    """
    if stack.n_layers <= 1:
        return ComplexImageFit(
            a=np.zeros(0, dtype=complex),
            beta=np.zeros(0, dtype=complex),
            rms=0.0,
        )

    h_min = float(np.min(stack.h))
    lam_min = lambda_min_factor / h_min
    lam_max = lambda_max_factor / h_min

    # Uniform spacing required by the matrix-pencil method.
    delta = (lam_max - lam_min) / (n_samples - 1)
    lam = lam_min + delta * np.arange(n_samples)
    g = reflection_gamma(stack, lam)
    g = g.astype(complex)

    # If Γ_1(λ) is essentially zero everywhere (degenerate case
    # ρ_1 = ρ_2 = … = ρ_n; matrix-pencil would receive a zero Hankel
    # block and divide by zero), short-circuit with an empty fit.
    if np.max(np.abs(g)) < 1e-12:
        return ComplexImageFit(
            a=np.zeros(0, dtype=complex),
            beta=np.zeros(0, dtype=complex),
            rms=0.0,
        )

    # Matrix-pencil method: form the Hankel pencil [Y0, Y1] of size
    # (n_samples - L) x L, with pencil parameter L ≈ n_samples/3.
    L = max(n_images + 1, n_samples // 3)
    L = min(L, n_samples - n_images - 1)
    if L < n_images + 1:
        # Not enough samples — fall back to fewer images.
        n_images = max(1, L - 1)

    rows = n_samples - L
    Y = np.zeros((rows, L + 1), dtype=complex)
    for i in range(rows):
        Y[i, :] = g[i:i + L + 1]
    Y0 = Y[:, :L]
    Y1 = Y[:, 1:L + 1]

    # SVD-based filtering: project both blocks onto the dominant
    # singular subspace of order n_images. Filter near-zero singular
    # values to avoid 1/0.
    U, S, Vh = np.linalg.svd(Y0, full_matrices=False)
    s_max = float(S.max()) if S.size else 0.0
    if s_max == 0.0:
        return ComplexImageFit(
            a=np.zeros(0, dtype=complex),
            beta=np.zeros(0, dtype=complex),
            rms=0.0,
        )
    keep_mask = S > s_max * 1e-10
    P = int(min(n_images, np.sum(keep_mask)))
    if P < 1:
        return ComplexImageFit(
            a=np.zeros(0, dtype=complex),
            beta=np.zeros(0, dtype=complex),
            rms=0.0,
        )
    Up = U[:, :P]
    Sp = S[:P]
    Vp = Vh[:P, :]

    # Y1 ≈ Up · diag(Sp) · Vp · Z, with Z holding the poles in its
    # eigenvalues. Equivalently the poles solve
    #   eig( pinv(Y0) · Y1 ) ≈ pinv(Sp) · Up^H · Y1 · Vp^H.
    Z = np.diag(1.0 / Sp) @ Up.conj().T @ Y1 @ Vp.conj().T
    poles = np.linalg.eigvals(Z)

    # Map poles back to the β coefficients: p_k = exp(-2 Δ β_k)
    #   → β_k = -log(p_k) / (2 Δ).
    # Drop poles outside the unit disc (non-physical, would blow up).
    mask = np.abs(poles) < 0.999
    poles = poles[mask]
    if poles.size == 0:
        return ComplexImageFit(
            a=np.zeros(0, dtype=complex),
            beta=np.zeros(0, dtype=complex),
            rms=float("nan"),
        )
    beta = -np.log(poles) / (2.0 * delta)

    # Coefficients a_k by linear least squares on the original samples.
    # g(λ_j) ≈ Σ a_k · exp(-2 λ_j β_k)
    A = np.exp(-2.0 * lam[:, None] * beta[None, :])  # (n_samples, P)
    a, *_ = np.linalg.lstsq(A, g, rcond=None)

    # Discard images with negative real part of β (non-decaying).
    keep = beta.real > 0.0
    beta = beta[keep]
    a = a[keep]
    if beta.size == 0:
        return ComplexImageFit(
            a=np.zeros(0, dtype=complex),
            beta=np.zeros(0, dtype=complex),
            rms=float("nan"),
        )
    A_keep = np.exp(-2.0 * lam[:, None] * beta[None, :])
    rms = float(np.sqrt(np.mean(np.abs(A_keep @ a - g) ** 2)))
    _log.info(
        "cim.fit: n_images_used=%d, rms=%.2e, n_layers=%d",
        beta.size, rms, stack.n_layers,
    )
    return ComplexImageFit(a=a, beta=beta, rms=rms)

solve_cim ¶

solve_cim(
    world: "World",
    engine: "Engine",
    *,
    n_images: int = 8,
    n_samples: int = 64
) -> FieldResult

Complex-Image-Method solver for layered soil.

Accepts :class:HomogeneousSoil, :class:TwoLayerSoil, and :class:MultiLayerSoil. The fit reduces to P = 0 on homogeneous soil (no extra images needed); for layered soils the closed-form CIM Green's function is used.