GF_Sommerfeld

class mqed.Dyadic_GF.GF_Sommerfeld.Greens_function_analytical(metal_epsi: complex, omega: float, eps_0: float = 1.0, qmax: float | None = None, epsabs: float = 1e-10, epsrel: float = 1e-10, limit: int = 400, split_propagating: bool = True)

Dyadic Green’s function for a planar interface with cylindrical symmetry.

Architecture Note (N-Layer Extensibility)

The reflection coefficients and z-wavevectors are computed via dedicated methods rather than inline lambdas. This design separates the integral machinery (Sommerfeld integrals, quadrature, tensor assembly — shared by all planar geometries) from the layer physics (Fresnel coefficients, phase factors — specific to the number of layers).

To extend to N layers, subclass this and override:

• _kz() — return kz in any layer

• _rs() / _rp() — return generalized reflection coefficients (e.g., recursive Tomas 1995 or transfer-matrix formalism)

The Sommerfeld integrals in compute_integrals() call these methods and remain unchanged.

s- and p-polarized Fresnel coefficients (two-layer):

\[r_s(q) = \frac{K_{z,0}(q) - K_{z,1}(q)}{K_{z,0}(q) + K_{z,1}(q)}\]

\[r_p(q) = \frac{\epsilon_1 K_{z,0}(q) - \epsilon_0 K_{z,1}(q)}{\epsilon_1 K_{z,0}(q) + \epsilon_0 K_{z,1}(q)}\]

where

\[K_{z,i}(q) = \sqrt{\epsilon_i k_0^2 - q^2}.\]

Total field is vacuum plus reflected part:

\[\overline{\overline{\mathbf{G}}}(\mathbf{r}_\alpha,\mathbf{r}_\beta,\omega) = \overline{\overline{\mathbf{G}}}_0(\mathbf{r}_\alpha,\mathbf{r}_\beta,\omega) + \overline{\overline{\mathbf{G}}}_{\text{refl}}^{(i)}(\mathbf{r}_\alpha,\mathbf{r}_\beta,\omega).\]

__init__(metal_epsi: complex, omega: float, eps_0: float = 1.0, qmax: float | None = None, epsabs: float = 1e-10, epsrel: float = 1e-10, limit: int = 400, split_propagating: bool = True)

Initializes the two-layer Green’s function calculator.

The two media are:

• Layer 0 (upper half-space, z > 0): permittivity eps_0 — typically vacuum (1.0).

• Layer 1 (lower half-space, z < 0): permittivity metal_epsi — typically a dispersive metal like Ag or Au.

Parameters:

• metal_epsi – Complex permittivity of the metal (layer 1).

• omega – Angular frequency [rad/s].

• eps_0 – Permittivity of the upper half-space (layer 0). Default 1.0 (vacuum).

• qmax – Upper integration limit for Sommerfeld integrals. None → integrate to infinity.

• epsabs – Absolute error tolerance for quad_vec.

• epsrel – Relative error tolerance for quad_vec.

• limit – Maximum number of adaptive subintervals.

• split_propagating – If True, split the integral at q = k0 (propagating/evanescent boundary) for better numerical accuracy.

static _beta_phys(eps: complex, k0: float, q)

Physical z-wavevector in a medium with permittivity eps.

\[K_{z}(q) = \sqrt{\epsilon\, k_0^2 - q^2}\]

The branch is chosen so that Im(K_z) >= 0 (evanescent waves decay away from the interface) and, when Im(K_z) ≈ 0, Re(K_z) >= 0 (propagating waves travel upward).

A tiny imaginary offset +i·10⁻¹² is added under the square root to lift the branch cut off the real axis — this is the standard i0⁺ prescription.

Parameters:

• eps – Complex permittivity of the medium.

• k0 – Free-space wavenumber ω/c [1/m].

• q – In-plane (transverse) wavevector component.

Returns:

Complex Kz value(s) with correct branch.

Note

Made a @staticmethod so it can be called without an instance (useful in parallel workers or standalone utilities). The k0 argument replaces the previous self.k0 access.

_kz(layer: int, q)

z-wavevector component in a given layer.

For the current two-layer geometry:

• layer 0 → upper half-space (permittivity self.eps_0)

• layer 1 → lower half-space (permittivity self.metal_epsi)

Override in an N-layer subclass to index into an array of permittivities: self.eps[layer].

Parameters:

• layer – Layer index (0 = upper, 1 = metal).

• q – In-plane wavevector.

Returns:

Complex Kz in the requested layer.

_rp(q)

p-polarized (TM) Fresnel reflection coefficient at the interface.

Two-layer formula:

\[r_p(q) = \frac{\epsilon_1 K_{z,0} - \epsilon_0 K_{z,1}} {\epsilon_1 K_{z,0} + \epsilon_0 K_{z,1}}\]

For N-layer extension, replace with the generalized reflection coefficient.

Parameters:

q – In-plane wavevector (scalar or array).

Returns:

Complex reflection coefficient.

_rs(q)

s-polarized (TE) Fresnel reflection coefficient at the interface.

Two-layer formula:

\[r_s(q) = \frac{K_{z,0} - K_{z,1}}{K_{z,0} + K_{z,1}}\]

For N-layer extension, replace with the generalized reflection coefficient computed via recursive Tomas (1995) relations or transfer-matrix method.

Parameters:

q – In-plane wavevector (scalar or array).

Returns:

Complex reflection coefficient.

calculate_total_Green_function(x: float, y: float, z1: float, z2: float)

Total Green’s function = vacuum + scattering.

\[\overline{\overline{\mathbf{G}}} = \overline{\overline{\mathbf{G}}}_0 + \overline{\overline{\mathbf{G}}}_{\text{SC}}^{(i)}.\]

Parameters:

• x – x-distance between the two points.

• y – y-distance between the two points.

• z1 – z-coordinate of the first point.

• z2 – z-coordinate of the second point.

Returns:

3x3 total Green’s tensor.

Return type:

np.ndarray

complex_quad(func: Callable[[float], complex | numpy.ndarray], a: float, qmax: float | None = None, epsabs: float | None = None, epsrel: float | None = None, limit: int | None = None, split_propagating: bool | None = None) → complex | numpy.ndarray

Integrates a complex-valued function using scipy’s quad_vec function. :param func: The complex-valued function to integrate. :type func: callable :param a: The lower limit of integration. :type a: float :param qmax: Maximum q value for integration. Defaults to np.inf. :type qmax: float :param epsabs: Absolute error tolerance. Defaults to 1e-10. :type epsabs: float, optional :param epsrel: Relative error tolerance. Defaults to 1e-10 :type epsrel: float, optional :param limit: Maximum number of subintervals. Defaults to 400. :type limit: int, optional :param returns: complex: The result of the integration.

compute_integrals(rho: float, z1: float, z2: float)

Compute the six Sommerfeld integrals in one vector-valued quadrature.

\[I_1 = \int_{0}^{\infty} dq\; R_s(q) \frac{q}{2K_{z,0}} J_0(q\rho) e^{iK_{z,0}(z_1+z_2)}\]

\[I_2 = \int_{0}^{\infty} dq\; R_s(q) \frac{q}{2K_{z,0}} J_2(q\rho) e^{iK_{z,0}(z_1+z_2)}\]

\[I_3 = \int_{0}^{\infty} dq\; R_p(q) \frac{qK_{z,0}}{2k_0^2} J_0(q\rho) e^{iK_{z,0}(z_1+z_2)}\]

\[I_4 = \int_{0}^{\infty} dq\; R_p(q) \frac{qK_{z,0}}{2k_0^2} J_2(q\rho) e^{iK_{z,0}(z_1+z_2)}\]

\[I_5 = \int_{0}^{\infty} dq\; R_p(q) \frac{i q^2}{k_0^2} J_1(q\rho) e^{iK_{z,0}(z_1+z_2)}\]

\[I_6 = \int_{0}^{\infty} dq\; R_p(q) \frac{q^3}{K_{z,0}k_0^2} J_0(q\rho) e^{iK_{z,0}(z_1+z_2)}\]

scatter_component(x: float, y: float, z1: float, z2: float)

Scattering part of the Green’s function.

\[\overline{\overline{\mathbf{G}}}_{\text{refl}}^{(i)}(\rho, \phi, z, z', \omega) = \int_{0}^{+\infty} \frac{i\,dk_{\rho}}{4\pi} \Big[ R_s\,\overline{\overline{\mathbf{M}}}_s + R_p\,\overline{\overline{\mathbf{M}}}_p \Big] e^{iK_{z,i}(k_{\rho}, \omega)(z+z')}.\]

Parameters:

• x – x-distance between the two points.

• y – y-distance between the two points.

• z1 – z-coordinate of the first point.

• z2 – z-coordinate of the second point.

scattering_p_component(x: float, y: float, z1: float, z2: float, integrals: numpy.ndarray | None = None)

p-polarized scattering tensor.

\[\begin{split}\overline{\overline{\mathbf{M}}}_p(k_{\rho}, \omega) = \frac{-k_{\rho}K_{z,i}}{2k_i^2} \begin{bmatrix} J_0 - \cos(2\phi)J_2 & -\sin(2\phi)J_2 & \frac{2ik_{\rho}}{K_{z,i}}\cos\phi\,J_1 \\ -\sin(2\phi)J_2 & J_0 + \cos(2\phi)J_2 & \frac{2ik_{\rho}}{K_{z,i}}\sin\phi\,J_1 \\ -\frac{2ik_{\rho}}{K_{z,i}}\cos\phi\,J_1 & -\frac{2ik_{\rho}}{K_{z,i}}\sin\phi\,J_1 & -\frac{2k_{\rho}^2}{K_{z,i}^2}J_0 \end{bmatrix}.\end{split}\]

Parameters:

• x – x-distance between the two points.

• y – y-distance between the two points.

• z1 – z-coordinate of the first point.

• z2 – z-coordinate of the second point.

Returns:

3x3 p-polarized scattering tensor.

Return type:

np.ndarray

scattering_s_component(x: float, y: float, z1: float, z2: float, integrals: numpy.ndarray | None = None)

s-polarized scattering tensor.

\[\begin{split}\overline{\overline{\mathbf{M}}}_s(k_{\rho}, \omega) = \frac{k_{\rho}}{2K_{z,i}} \begin{bmatrix} J_0 + \cos(2\phi)J_2 & \sin(2\phi)J_2 & 0 \\ \sin(2\phi)J_2 & J_0 - \cos(2\phi)J_2 & 0 \\ 0 & 0 & 0 \end{bmatrix}.\end{split}\]

Parameters:

• x – x-distance between the two points.

• y – y-distance between the two points.

• z1 – z-coordinate of the first point.

• z2 – z-coordinate of the second point.

Returns:

3x3 s-polarized scattering tensor.

Return type:

np.ndarray

vacuum_component(x: float, y: float, z1: float, z2: float)

Vacuum dyadic Green’s function.

\[\overline{\overline{\mathbf{G}}}_0(\mathbf{r}_\alpha,\mathbf{r}_\beta,\omega) = \frac{e^{ik_0 R_{\alpha\beta}}}{4\pi R_{\alpha\beta}}\Big[ \left(\overline{\overline{\mathbf{I}}}_3-\mathbf{e}_\mathrm{R}\mathbf{e}_\mathrm{R}\right) + \left(3\mathbf{e}_\mathrm{R}\mathbf{e}_\mathrm{R}-\overline{\overline{\mathbf{I}}}_3\right)\left(\frac{1}{(k_0 R_{\alpha\beta})^{2}}-\frac{i}{k_0 R_{\alpha\beta}}\right)\Big].\]

Parameters:

• x – x-distance between the two points.

• y – y-distance between the two points.

• z1 – z-coordinate of the first point.

• z2 – z-coordinate of the second point.

Returns:

3x3 vacuum Green’s tensor.

Return type:

np.ndarray