sommerfeld_singularities

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

• find_poles_by_winding() - locate zeros of an Airy/modal denominator by recursively applying the argument principle on rectangular boxes.

• residue_vector_by_contour() - estimate vector residues of a spectral kernel after a pole location has been found.

• 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.

class mqed.Dyadic_GF.sommerfeld_singularities.BranchCutConfig(branch_layer: int, t_limit: float, side_offset: float = 0.0, epsabs: float = 1e-08, epsrel: float = 1e-08, limit: int = 200, use_hankel: bool = True, include_two_sides: bool = True)

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

\[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.

Parameters:

• 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.

__delattr__(name)

Implement delattr(self, name).

__eq__(other)

Return self==value.

__init__(branch_layer: int, t_limit: float, side_offset: float = 0.0, epsabs: float = 1e-08, epsrel: float = 1e-08, limit: int = 200, use_hankel: bool = True, include_two_sides: bool = True) → None

__repr__()

Return repr(self).

__setattr__(name, value)

Implement setattr(self, name, value).

class mqed.Dyadic_GF.sommerfeld_singularities.ComplexBox(real_min: float, real_max: float, imag_min: float, imag_max: float, depth: int = 0)

Axis-aligned rectangular search box in the complex q plane.

Parameters:

• 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.

__delattr__(name)

Implement delattr(self, name).

__eq__(other)

Return self==value.

__init__(real_min: float, real_max: float, imag_min: float, imag_max: float, depth: int = 0) → None

__repr__()

Return repr(self).

__setattr__(name, value)

Implement setattr(self, name, value).

class mqed.Dyadic_GF.sommerfeld_singularities.PoleResidue(pole: SommerfeldPole, residues: numpy.ndarray, contour_radius: float, method: str = 'contour')

Residue of the seven Bessel-free kernels at one pole.

If a Bessel-free kernel has a simple pole at q_p,

\[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.

__delattr__(name)

Implement delattr(self, name).

__eq__(other)

Return self==value.

__init__(pole: SommerfeldPole, residues: numpy.ndarray, contour_radius: float, method: str = 'contour') → None

__repr__()

Return repr(self).

__setattr__(name, value)

Implement setattr(self, name, value).

class mqed.Dyadic_GF.sommerfeld_singularities.PoleSearchConfig(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-06, winding_tol: float = 0.35, denominator_floor: float = 1e-12, root_tol: float = 1e-10, residual_tol: float = 1e-06, dedup_tol: float = 1e-05, max_boxes: int = 20000, skip_branch_points: Sequence[complex] = <factory>, branch_point_guard: float = 0.0)

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

\[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 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 from_k0() to specify the region in dimensionless q/k0 units.

Parameters:

• 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.

__delattr__(name)

Implement delattr(self, name).

__eq__(other)

Return self==value.

__init__(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-06, winding_tol: float = 0.35, denominator_floor: float = 1e-12, root_tol: float = 1e-10, residual_tol: float = 1e-06, dedup_tol: float = 1e-05, max_boxes: int = 20000, skip_branch_points: Sequence[complex] = <factory>, branch_point_guard: float = 0.0) → None

__repr__()

Return repr(self).

__setattr__(name, value)

Implement setattr(self, name, value).

class mqed.Dyadic_GF.sommerfeld_singularities.SommerfeldPole(q: complex, polarization: str, residual: float, winding_number: int = 1, box: ComplexBox | None = None, derivative: complex | None = None)

One located Sommerfeld pole/root.

Parameters:

• q – Complex in-plane pole wave number in SI units.

• polarization – Polarization label, normally "s" or "p".

• residual – Absolute denominator 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 – Numerical derivative of the searched denominator at q.

__delattr__(name)

Implement delattr(self, name).

__eq__(other)

Return self==value.

__init__(q: complex, polarization: str, residual: float, winding_number: int = 1, box: ComplexBox | None = None, derivative: complex | None = None) → None

__repr__()

Return repr(self).

__setattr__(name, value)

Implement setattr(self, name, value).

mqed.Dyadic_GF.sommerfeld_singularities.branch_cut_integral_kz(kernel: Callable[[complex], complex | numpy.ndarray], branch_wavenumber: complex, rho: float, orders: Sequence[int], config: BranchCutConfig) → numpy.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.

mqed.Dyadic_GF.sommerfeld_singularities.branch_cut_samples_kz(kernel: Callable[[complex], complex | numpy.ndarray], branch_wavenumber: complex, t_values: numpy.ndarray, side_offset: float = 0.0, side: int = 1) → tuple[numpy.ndarray, numpy.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.

mqed.Dyadic_GF.sommerfeld_singularities.complex_derivative(func: Callable[[complex], complex], z: complex, step: float) → complex

Symmetric finite-difference derivative in the complex plane.

\[f'(z) \approx \frac{f(z+h)-f(z-h)}{2h}.\]

mqed.Dyadic_GF.sommerfeld_singularities.contour_points(box: ComplexBox, points_per_side: int) → numpy.ndarray

Counter-clockwise contour points around box without duplicate corners.

Parameters:

• 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.

mqed.Dyadic_GF.sommerfeld_singularities.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.

mqed.Dyadic_GF.sommerfeld_singularities.pole_integral_contribution(residue: PoleResidue, rho: float, orders: Sequence[int], prefactor: complex = 3.141592653589793j, use_hankel: bool = True) → numpy.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.

mqed.Dyadic_GF.sommerfeld_singularities.residue_vector_by_contour(kernel: Callable[[complex], complex | numpy.ndarray], pole: SommerfeldPole, radius: float, points: int = 96) → PoleResidue

Estimate vector residues by a circular contour integral.

\[A_m = \frac{1}{2\pi i}\oint_{|q-q_p|=r} F_m(q)\,dq.\]

Parameters:

• 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.

mqed.Dyadic_GF.sommerfeld_singularities.residue_vector_by_limit(kernel: Callable[[complex], complex | numpy.ndarray], pole: SommerfeldPole, step: float) → PoleResidue

Estimate residues from the simple-pole limit.

\[A_m = \lim_{q\to q_p} (q-q_p)F_m(q).\]

Four axial approach directions are averaged to reduce directional bias.

mqed.Dyadic_GF.sommerfeld_singularities.winding_number(values: numpy.ndarray, denominator_floor: float = 1e-12) → int

Integer winding number of complex samples around the origin.

The discrete estimate is

\[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.