"""Pole and branch-cut utilities for N-layer Sommerfeld Green functions.
The routines in this module are intentionally independent of the concrete
``NLayerGreenFunction`` class. They work with callables such as
``denominator(q)`` or ``kernel(q)`` so they can be reused by planar multilayer,
Sommerfeld-interface, or transmission-line Green-function implementations.
The most useful entry points are
* :func:`find_poles_by_winding` - locate zeros of an Airy/modal denominator by
recursively applying the argument principle on rectangular boxes.
* :func:`residue_vector_by_contour` - estimate vector residues of a spectral
kernel after a pole location has been found.
* :func:`branch_cut_integral_kz` - evaluate a diagnostic branch-cut contour in a
vertical-wavenumber variable. This is useful for validating branch-cut DCIM
ideas before making them the production solver.
All lengths are SI, so a pole in the q plane has units of 1/m.
"""
from __future__ import annotations
import warnings
from dataclasses import dataclass, field
from math import pi
from numbers import Integral
from typing import Callable, Sequence, Union
import numpy as np
from scipy.integrate import quad_vec
from scipy.optimize import root
from scipy.special import hankel1, jv
ArrayLikeComplex = Union[complex, np.ndarray]
[docs]
@dataclass(frozen=True)
class ComplexBox:
"""Axis-aligned rectangular search box in the complex q plane.
Args:
real_min: Lower bound for ``Re(q)`` in SI units.
real_max: Upper bound for ``Re(q)`` in SI units.
imag_min: Lower bound for ``Im(q)`` in SI units.
imag_max: Upper bound for ``Im(q)`` in SI units.
depth: Recursion depth used by the winding-number subdivision search.
"""
real_min: float
real_max: float
imag_min: float
imag_max: float
depth: int = 0
@property
def width(self) -> float:
return float(self.real_max - self.real_min)
@property
def height(self) -> float:
return float(self.imag_max - self.imag_min)
@property
def center(self) -> complex:
return complex(0.5 * (self.real_min + self.real_max), 0.5 * (self.imag_min + self.imag_max))
@property
def diameter(self) -> float:
return float(np.hypot(self.width, self.height))
def subdivide(self) -> tuple["ComplexBox", "ComplexBox", "ComplexBox", "ComplexBox"]:
mid_r = 0.5 * (self.real_min + self.real_max)
mid_i = 0.5 * (self.imag_min + self.imag_max)
d = self.depth + 1
return (
ComplexBox(self.real_min, mid_r, self.imag_min, mid_i, d),
ComplexBox(mid_r, self.real_max, self.imag_min, mid_i, d),
ComplexBox(self.real_min, mid_r, mid_i, self.imag_max, d),
ComplexBox(mid_r, self.real_max, mid_i, self.imag_max, d),
)
[docs]
@dataclass(frozen=True)
class PoleSearchConfig:
r"""Controls recursive pole search by the argument principle.
The search locates zeros of a scalar denominator ``D(q)`` by sampling boxes
in the complex ``q`` plane and applying
.. math::
N_Z-N_P = \frac{1}{2\pi i}\oint_{\partial B}\frac{D'(q)}{D(q)}dq
= \frac{\Delta \arg D(q)}{2\pi}.
For the Airy denominators used by :class:`NLayerGreenFunction`, the relevant
singularities are zeros rather than poles of ``D``. Nonzero winding boxes
are subdivided until they are small enough for a local complex root solve.
The default region is intended for outgoing-wave poles with positive
``Re(q)`` and negative ``Im(q)``. Use :meth:`from_k0` to specify the region
in dimensionless ``q/k0`` units.
Args:
real_min: Lower search bound for ``Re(q)`` in SI units.
real_max: Upper search bound for ``Re(q)`` in SI units.
imag_min: Lower search bound for ``Im(q)`` in SI units.
imag_max: Upper search bound for ``Im(q)`` in SI units.
contour_points_per_side: Number of denominator samples on each box side.
max_depth: Maximum recursive subdivision depth.
min_box_size: Minimum box width/height before accepting a candidate.
winding_tol: Absolute winding below this threshold is treated as zero.
denominator_floor: Near-zero denominator threshold that forces subdivision.
root_tol: Tolerance passed to the local root solver.
residual_tol: Maximum accepted ``abs(D(q_p))`` after refinement.
dedup_tol: Distance below which two refined roots are treated as one pole.
max_boxes: Safety bound on the number of boxes tested.
skip_branch_points: Branch points to guard against false pole detections.
branch_point_guard: Exclusion radius around branch points in SI units.
"""
real_min: float
real_max: float
imag_min: float
imag_max: float
contour_points_per_side: int = 24
max_depth: int = 10
min_box_size: float = 1e-6
winding_tol: float = 0.35
denominator_floor: float = 1e-12
root_tol: float = 1e-10
residual_tol: float = 1e-6
dedup_tol: float = 1e-5
max_boxes: int = 20000
skip_branch_points: Sequence[complex] = field(default_factory=tuple)
branch_point_guard: float = 0.0
def __post_init__(self) -> None:
finite_fields = {
"real_min": self.real_min,
"real_max": self.real_max,
"imag_min": self.imag_min,
"imag_max": self.imag_max,
"min_box_size": self.min_box_size,
"winding_tol": self.winding_tol,
"denominator_floor": self.denominator_floor,
"root_tol": self.root_tol,
"residual_tol": self.residual_tol,
"dedup_tol": self.dedup_tol,
"branch_point_guard": self.branch_point_guard,
}
if not all(np.isfinite(value) for value in finite_fields.values()):
raise ValueError("Pole-search bounds and tolerances must be finite.")
if self.real_min >= self.real_max or self.imag_min >= self.imag_max:
raise ValueError("Pole-search bounds must have positive width and height.")
integer_fields = {
"contour_points_per_side": self.contour_points_per_side,
"max_depth": self.max_depth,
"max_boxes": self.max_boxes,
}
if not all(isinstance(value, Integral) for value in integer_fields.values()):
raise ValueError("Pole-search counts must be integers.")
if not 2 <= self.contour_points_per_side <= 4096:
raise ValueError("contour_points_per_side must be between 2 and 4096.")
if self.max_depth < 0 or self.max_boxes <= 0:
raise ValueError("max_depth must be non-negative and max_boxes must be positive.")
if any(value < 0 for value in (
self.min_box_size,
self.winding_tol,
self.denominator_floor,
self.root_tol,
self.residual_tol,
self.dedup_tol,
self.branch_point_guard,
)):
raise ValueError("Pole-search sizes and tolerances must be non-negative.")
if not all(np.isfinite(point) for point in self.skip_branch_points):
raise ValueError("skip_branch_points must contain only finite values.")
@classmethod
def from_k0(
cls,
k0: float,
real_min_factor: float = 0.0,
real_max_factor: float = 6.0,
imag_min_factor: float = -2.0,
imag_max_factor: float = 1e-3,
**kwargs,
) -> "PoleSearchConfig":
k0_abs = float(abs(k0))
if k0_abs <= 0:
raise ValueError("k0 must be nonzero.")
return cls(
real_min=real_min_factor * k0_abs,
real_max=real_max_factor * k0_abs,
imag_min=imag_min_factor * k0_abs,
imag_max=imag_max_factor * k0_abs,
min_box_size=float(kwargs.pop("min_box_size_factor", 1e-6)) * k0_abs,
dedup_tol=float(kwargs.pop("dedup_tol_factor", 1e-5)) * k0_abs,
branch_point_guard=float(kwargs.pop("branch_point_guard_factor", 1e-4)) * k0_abs,
**kwargs,
)
[docs]
@dataclass(frozen=True)
class SommerfeldPole:
"""One located Sommerfeld pole/root.
Args:
q: Complex in-plane pole wave number in SI units.
polarization: Polarization label, normally ``"s"`` or ``"p"``.
residual: Absolute search-characteristic residual at ``q``.
winding_number: Winding number of the terminal box that produced the pole.
box: Terminal search box associated with this pole, when available.
derivative: Diagnostic numerical slope of the search characteristic at
``q``. It is not an analytic mode-normalization derivative.
"""
q: complex
polarization: str
residual: float
winding_number: int = 1
box: ComplexBox | None = None
derivative: complex | None = None
[docs]
@dataclass(frozen=True)
class PoleResidue:
r"""Residue of the seven Bessel-free kernels at one pole.
If a Bessel-free kernel has a simple pole at ``q_p``,
.. math::
F_m(q)=\frac{A_m}{q-q_p}+F_m^{\mathrm{reg}}(q),
then ``residues[m]`` stores ``A_m`` before multiplication by Bessel or Hankel
functions.
"""
pole: SommerfeldPole
residues: np.ndarray
contour_radius: float
method: str = "contour"
[docs]
@dataclass(frozen=True)
class BranchCutConfig:
r"""Controls branch-cut diagnostic integration in a ``kz`` variable.
``branch_layer`` selects the layer whose light-line branch point
``sqrt(eps_l) k0`` is used. The path is parameterized by the vertical
wavenumber
.. math::
k_{z,l}=t \pm i\eta,\qquad
q_\pm(t)=\sqrt{k_l^2-k_{z,l}^2},\qquad
\frac{dq_\pm}{dt}=-\frac{k_{z,l}}{q_\pm}.
For lossy stacks and arbitrary branch-cut conventions this should be treated
as a diagnostic/validation tool rather than a universal production solver.
Args:
branch_layer: Layer index whose light-line branch point is sampled.
t_limit: Integration half-width for ``t`` in SI units.
side_offset: Offset ``eta`` separating the two branch-cut sides.
epsabs: Absolute quadrature tolerance.
epsrel: Relative quadrature tolerance.
limit: Maximum number of ``quad_vec`` subintervals.
use_hankel: Use outgoing Hankel functions instead of Bessel functions.
include_two_sides: Return the two-side jump when true; otherwise sample
only the ``+`` side.
"""
branch_layer: int
t_limit: float
side_offset: float = 0.0
epsabs: float = 1e-8
epsrel: float = 1e-8
limit: int = 200
use_hankel: bool = True
include_two_sides: bool = True
@classmethod
def from_k0(
cls,
k0: float,
branch_layer: int,
t_limit_factor: float = 8.0,
side_offset_factor: float = 1e-6,
**kwargs,
) -> "BranchCutConfig":
k0_abs = float(abs(k0))
if k0_abs <= 0:
raise ValueError("k0 must be nonzero.")
return cls(
branch_layer=branch_layer,
t_limit=t_limit_factor * k0_abs,
side_offset=side_offset_factor * k0_abs,
**kwargs,
)
def _as_1d_complex(value: ArrayLikeComplex) -> np.ndarray:
arr = np.asarray(value, dtype=complex)
if arr.ndim == 0:
return arr.reshape(1)
return arr.reshape(-1)
[docs]
def contour_points(box: ComplexBox, points_per_side: int) -> np.ndarray:
"""Counter-clockwise contour points around ``box`` without duplicate corners.
Args:
box: Rectangular complex-q domain.
points_per_side: Number of samples per side, excluding the next corner.
Returns:
Closed contour samples with the starting point repeated once at the end.
"""
n = int(points_per_side)
if n < 2:
raise ValueError("points_per_side must be at least 2.")
r0, r1 = box.real_min, box.real_max
i0, i1 = box.imag_min, box.imag_max
bottom = np.linspace(r0 + 1j * i0, r1 + 1j * i0, n, endpoint=False)
right = np.linspace(r1 + 1j * i0, r1 + 1j * i1, n, endpoint=False)
top = np.linspace(r1 + 1j * i1, r0 + 1j * i1, n, endpoint=False)
left = np.linspace(r0 + 1j * i1, r0 + 1j * i0, n, endpoint=False)
return np.concatenate((bottom, right, top, left, np.array([r0 + 1j * i0], dtype=complex)))
[docs]
def winding_number(values: np.ndarray, denominator_floor: float = 1e-12) -> int:
r"""Integer winding number of complex samples around the origin.
The discrete estimate is
.. math::
w \approx \operatorname{round}\left[\frac{\Delta\arg D}{2\pi}\right].
If a contour sample passes too close to the origin, the function returns a
nonzero value to force subdivision instead of accidentally discarding a pole.
"""
vals = np.asarray(values, dtype=complex)
vals = vals[np.isfinite(vals)]
if vals.size < 4:
return 0
if np.min(np.abs(vals)) < denominator_floor:
# The contour probably passes too close to a zero/pole. Returning a
# nonzero value forces subdivision, which is safer than discarding it.
return 1
phases = np.unwrap(np.angle(vals))
return int(np.rint((phases[-1] - phases[0]) / (2 * np.pi)))
def _near_branch_point(box: ComplexBox, branch_points: Sequence[complex], guard: float) -> bool:
if guard <= 0 or not branch_points:
return False
c = box.center
expanded_real = (box.real_min - guard, box.real_max + guard)
expanded_imag = (box.imag_min - guard, box.imag_max + guard)
for bp in branch_points:
if expanded_real[0] <= np.real(bp) <= expanded_real[1] and expanded_imag[0] <= np.imag(bp) <= expanded_imag[1]:
return True
if abs(c - bp) <= guard:
return True
return False
def _bisect_longest(box: ComplexBox) -> tuple[ComplexBox, ComplexBox]:
depth = box.depth + 1
if box.width >= box.height:
mid = 0.5 * (box.real_min + box.real_max)
return (
ComplexBox(box.real_min, mid, box.imag_min, box.imag_max, depth),
ComplexBox(mid, box.real_max, box.imag_min, box.imag_max, depth),
)
mid = 0.5 * (box.imag_min + box.imag_max)
return (
ComplexBox(box.real_min, box.real_max, box.imag_min, mid, depth),
ComplexBox(box.real_min, box.real_max, mid, box.imag_max, depth),
)
def _box_winding(denominator: Callable[[complex], complex], box: ComplexBox, config: PoleSearchConfig) -> int:
pts = contour_points(box, config.contour_points_per_side)
try:
vals = np.array([denominator(q) for q in pts], dtype=complex)
except Exception:
return 0
if not np.all(np.isfinite(vals)):
return 0
return winding_number(vals, config.denominator_floor)
def _refine_root(
denominator: Callable[[complex], complex],
start: complex,
tol: float,
) -> tuple[complex, float, bool]:
def f_xy(xy: np.ndarray) -> np.ndarray:
z = complex(float(xy[0]), float(xy[1]))
value = denominator(z)
return np.array([np.real(value), np.imag(value)], dtype=float)
sol = root(f_xy, np.array([np.real(start), np.imag(start)], dtype=float), tol=tol)
q = complex(float(sol.x[0]), float(sol.x[1]))
res = float(abs(denominator(q)))
return q, res, bool(sol.success)
[docs]
def complex_derivative(func: Callable[[complex], complex], z: complex, step: float) -> complex:
r"""Symmetric finite-difference derivative in the complex plane.
.. math::
f'(z) \approx \frac{f(z+h)-f(z-h)}{2h}.
"""
h = float(step)
if h <= 0:
h = max(1e-9, 1e-7 * max(1.0, abs(z)))
return (func(z + h) - func(z - h)) / (2.0 * h)
[docs]
def find_poles_by_winding(
denominator: Callable[[complex], complex],
polarization: str,
config: PoleSearchConfig,
) -> list[SommerfeldPole]:
"""Find denominator zeros in a rectangle using recursive winding numbers.
The algorithm implements the spirit of the pole-locating method described in
the FIPWA literature: use contour information to find all boxes that contain
singularities/zeros, recursively subdivide them, then refine the candidates
with a local nonlinear root solve.
"""
root_box = ComplexBox(config.real_min, config.real_max, config.imag_min, config.imag_max, 0)
queue: list[ComplexBox] = [root_box]
candidates: list[tuple[ComplexBox, int]] = []
boxes_tested = 0
while queue:
box = queue.pop()
boxes_tested += 1
if boxes_tested > config.max_boxes:
warnings.warn(
"Pole search reached max_boxes and is returning partial results. "
"Narrow the search region or increase max_boxes for a complete search.",
RuntimeWarning,
stacklevel=2,
)
break
if _near_branch_point(box, config.skip_branch_points, config.branch_point_guard):
# A coarse search box can contain both a guarded branch point and a
# distinct modal pole. Subdivide first so only terminal boxes local
# to the branch point are discarded.
if (
box.depth < config.max_depth
and max(box.width, box.height) > config.min_box_size
and max(box.width, box.height) > config.branch_point_guard
):
guarded_children = []
regular_children = []
for child in _bisect_longest(box):
if _near_branch_point(
child,
config.skip_branch_points,
config.branch_point_guard,
):
guarded_children.append(child)
else:
regular_children.append(child)
# The queue is LIFO, so regular regions are searched before
# further refinement of the excluded branch neighborhood.
queue.extend(guarded_children)
queue.extend(regular_children)
continue
w = _box_winding(denominator, box, config)
if abs(w) < config.winding_tol:
continue
if box.depth >= config.max_depth or max(box.width, box.height) <= config.min_box_size:
candidates.append((box, int(np.sign(w)) if w != 0 else 1))
continue
queue.extend(box.subdivide())
poles: list[SommerfeldPole] = []
for box, w in candidates:
q0 = box.center
q_refined, residual, success = _refine_root(denominator, q0, config.root_tol)
if not success and residual > config.residual_tol:
# Keep the box center only if it is already a plausible denominator minimum.
q_refined = q0
residual = float(abs(denominator(q_refined)))
if residual > config.residual_tol:
continue
if not (config.real_min <= np.real(q_refined) <= config.real_max):
continue
if not (config.imag_min <= np.imag(q_refined) <= config.imag_max):
continue
if _near_branch_point(
ComplexBox(np.real(q_refined), np.real(q_refined), np.imag(q_refined), np.imag(q_refined)),
config.skip_branch_points,
config.branch_point_guard,
):
continue
if any(abs(q_refined - pole.q) <= config.dedup_tol for pole in poles):
continue
deriv = complex_derivative(denominator, q_refined, max(config.min_box_size, 1e-7 * max(1.0, abs(q_refined))))
poles.append(
SommerfeldPole(
q=q_refined,
polarization=polarization,
residual=residual,
winding_number=w,
box=box,
derivative=deriv,
)
)
poles.sort(key=lambda p: (np.real(p.q), np.imag(p.q), p.polarization))
return poles
[docs]
def residue_vector_by_contour(
kernel: Callable[[complex], ArrayLikeComplex],
pole: SommerfeldPole,
radius: float,
points: int = 96,
) -> PoleResidue:
r"""Estimate vector residues by a circular contour integral.
.. math::
A_m = \frac{1}{2\pi i}\oint_{|q-q_p|=r} F_m(q)\,dq.
Args:
kernel: Callable returning one or more Bessel-free kernel values.
pole: Located pole around which the contour is drawn.
radius: Positive contour radius in SI units.
points: Number of endpoint-free trapezoid samples on the circle.
"""
if radius <= 0:
raise ValueError("radius must be positive.")
n = int(points)
if n < 16:
raise ValueError("points must be at least 16.")
theta = np.linspace(0.0, 2.0 * np.pi, n, endpoint=False)
q = pole.q + radius * np.exp(1j * theta)
dq_dtheta = 1j * radius * np.exp(1j * theta)
samples = np.array([_as_1d_complex(kernel(qq)) for qq in q], dtype=complex)
# Periodic trapezoid rule on endpoint-free samples. Using ``np.trapezoid``
# with endpoint-free theta would omit the closing segment.
dtheta = 2.0 * np.pi / n
integral = dtheta * np.sum(samples * dq_dtheta[:, None], axis=0)
residues = integral / (2.0j * np.pi)
return PoleResidue(pole=pole, residues=residues, contour_radius=radius, method="contour")
[docs]
def residue_vector_by_limit(
kernel: Callable[[complex], ArrayLikeComplex],
pole: SommerfeldPole,
step: float,
) -> PoleResidue:
r"""Estimate residues from the simple-pole limit.
.. math::
A_m = \lim_{q\to q_p} (q-q_p)F_m(q).
Four axial approach directions are averaged to reduce directional bias.
"""
if step <= 0:
raise ValueError("step must be positive.")
directions = np.array([1.0, -1.0, 1.0j, -1.0j], dtype=complex)
samples = []
for direction in directions:
q = pole.q + step * direction
samples.append(step * direction * _as_1d_complex(kernel(q)))
residues = np.mean(np.array(samples, dtype=complex), axis=0)
return PoleResidue(pole=pole, residues=residues, contour_radius=step, method="limit")
[docs]
def pole_integral_contribution(
residue: PoleResidue,
rho: float,
orders: Sequence[int],
prefactor: complex = 1j * pi,
use_hankel: bool = True,
) -> np.ndarray:
"""Convert kernel residues to scalar Sommerfeld-integral pole terms.
The default prefactor ``i*pi`` is the common half-line/Hankel contour factor
for outgoing waves, but sign conventions differ. Treat this routine as an
explicit diagnostic unless you have validated the prefactor against a known
analytic case for your exact Green-function normalization.
"""
q = residue.pole.q
values = []
for res, order in zip(residue.residues, orders):
special = hankel1(order, q * rho) if use_hankel else jv(order, q * rho)
values.append(prefactor * res * special)
return np.asarray(values, dtype=complex)
def _sqrt_physical(value: complex) -> complex:
root_value = np.lib.scimath.sqrt(value)
if np.imag(root_value) < 0 or (abs(np.imag(root_value)) < 1e-18 and np.real(root_value) < 0):
root_value = -root_value
return complex(root_value)
[docs]
def branch_cut_integral_kz(
kernel: Callable[[complex], ArrayLikeComplex],
branch_wavenumber: complex,
rho: float,
orders: Sequence[int],
config: BranchCutConfig,
) -> np.ndarray:
"""Diagnostic branch-cut integral parameterized by a layer vertical wavenumber.
This evaluates
integral [K(q_+) B(q_+ rho) dq_+/dt - K(q_-) B(q_- rho) dq_-/dt] dt
on ``t in [-t_limit, t_limit]`` if ``include_two_sides`` is true. With one
side it evaluates only the ``+`` path. The result is a vector with the same
length as ``orders`` and the kernel vector.
"""
orders = list(orders)
tmax = float(config.t_limit)
offset = float(config.side_offset)
if tmax <= 0:
raise ValueError("t_limit must be positive.")
def mapped_path(t: float, sign: float) -> tuple[complex, complex]:
kz = complex(float(t), sign * offset)
q = _sqrt_physical(branch_wavenumber**2 - kz**2)
if abs(q) < 1e-30:
dqdt = 0.0 + 0.0j
else:
dqdt = -kz / q
return q, dqdt
def side_value(t: float, sign: float) -> np.ndarray:
q, dqdt = mapped_path(t, sign)
kvec = _as_1d_complex(kernel(q))
if kvec.size != len(orders):
raise ValueError("kernel vector length must match the number of Bessel/Hankel orders.")
special = np.array(
[hankel1(order, q * rho) if config.use_hankel else jv(order, q * rho) for order in orders],
dtype=complex,
)
return kvec * special * dqdt
def integrand(t: float) -> np.ndarray:
plus = side_value(t, +1.0)
if not config.include_two_sides:
return plus
minus = side_value(t, -1.0)
return plus - minus
value, _ = quad_vec(
integrand,
-tmax,
tmax,
epsabs=config.epsabs,
epsrel=config.epsrel,
limit=config.limit,
)
return np.asarray(value, dtype=complex)
[docs]
def branch_cut_samples_kz(
kernel: Callable[[complex], ArrayLikeComplex],
branch_wavenumber: complex,
t_values: np.ndarray,
side_offset: float = 0.0,
side: int = +1,
) -> tuple[np.ndarray, np.ndarray]:
"""Sample a kernel along a kz-parameterized branch-cut side.
Returns ``(q_values, samples)``. This is intended for matrix-pencil fitting
experiments on a Sommerfeld branch cut.
"""
t_values = np.asarray(t_values, dtype=float)
q_values = np.empty_like(t_values, dtype=complex)
samples = []
for idx, t in enumerate(t_values):
kz = complex(float(t), float(side) * float(side_offset))
q = _sqrt_physical(branch_wavenumber**2 - kz**2)
q_values[idx] = q
samples.append(_as_1d_complex(kernel(q)))
return q_values, np.asarray(samples, dtype=complex)