Line data Source code
1 : # SPDX-FileCopyrightText: 2024 PairInteraction Developers
2 : # SPDX-License-Identifier: LGPL-3.0-or-later
3 :
4 : # ruff: noqa: N802, N806
5 :
6 1 : from __future__ import annotations
7 :
8 1 : import math
9 1 : from typing import TYPE_CHECKING, Literal, TypeVar
10 :
11 1 : import numpy as np
12 1 : import scipy.constants as const
13 1 : from numba import njit
14 1 : from scipy.integrate import quad
15 :
16 1 : from pairinteraction.green_tensor.bessel_function import (
17 : cached_bessel_function_0,
18 : cached_bessel_function_1,
19 : cached_bessel_function_2,
20 : )
21 :
22 : if TYPE_CHECKING:
23 : from collections.abc import Callable
24 :
25 : from pairinteraction.units import NDArray
26 :
27 : Entries = Literal["xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz"]
28 :
29 : from typing_extensions import ParamSpec
30 :
31 : P = ParamSpec("P")
32 : R = TypeVar("R")
33 :
34 : def njit(cache: bool) -> Callable[[Callable[P, R]], Callable[P, R]]: ... # type: ignore [no-redef]
35 :
36 :
37 1 : __all__ = ["dynamic_green_tensor_homogeneous", "dynamic_green_tensor_scattered"]
38 :
39 :
40 1 : def dynamic_green_tensor_homogeneous(
41 : pos1: NDArray, pos2: NDArray, omega: float, epsilon0: complex, *, only_real_part: bool = False
42 : ) -> NDArray:
43 : r"""Homogeneous Green Tensor for two atoms in cartesian coordinates in an infinite homogeneous medium.
44 :
45 : The function used is from equation 2 of the paper:
46 : "Dispersionless subradiant photon storage in one-dimensional emitter chains"
47 : https://doi.org/10.1103/PhysRevA.108.L051702
48 :
49 : We calculate the scaled Green Tensor, i.e.
50 : .. math::
51 : \frac{\omega^2}{\hbar \epsilon_0 c^2} G(r_\alpha, r_\beta, \omega)
52 :
53 : Args:
54 : pos1: Position vector of atom A in meters
55 : pos2: Position vector of atom B in meters
56 : omega: Angular frequency (i.e. 2*pi*f) in 1/s
57 : epsilon0: Electric permittivity of the medium (dimensionless, complex)
58 : only_real_part: If True, only the real part of the Green tensor is calculated (default: False)
59 :
60 : Returns: The 3x3 scaled homogeneous Green Tensor (general complex values) (m^(-3) [hbar]^(-1) [epsilon_0]^(-1))
61 : :math:`\omega^2 / (\hbar \epsilon_0 c^2) G(r_\alpha, r_\beta, \omega)`
62 :
63 : """
64 1 : k_vac = omega / const.c # magnitude of wave vector in vacuum
65 1 : k0 = k_vac * np.sqrt(epsilon0) # magnitude of wave vector in medium with permittivity epsilon0
66 1 : distance_vec = pos2 - pos1
67 1 : distance = np.linalg.norm(distance_vec)
68 :
69 : # this is missing a 1/k0**2 compared to the paper, we absorb this in the scaled prefactor, see line below
70 1 : prefactor = np.exp(1j * k0 * distance) / (4 * np.pi * distance**3)
71 1 : prefactor *= 1 / (epsilon0 * const.epsilon_0 * const.hbar) # epsilon0 from missing k0^2 compared to omega^2/c^2
72 1 : prefactor *= -1 # minus sign from H = - \hbar \sum V_{\alpha\beta} ...
73 :
74 1 : result: NDArray = prefactor * (
75 : (k0**2 * distance**2 + 1j * k0 * distance - 1) * np.eye(3)
76 : + (-(k0**2) * distance**2 - 3j * k0 * distance + 3) * np.outer(distance_vec, distance_vec) / distance**2
77 : )
78 :
79 1 : if only_real_part:
80 1 : return np.real(result)
81 0 : return result
82 :
83 :
84 1 : def dynamic_green_tensor_scattered(
85 : pos1: NDArray,
86 : pos2: NDArray,
87 : z1: float,
88 : z2: float,
89 : omega: float,
90 : epsilon0: complex,
91 : epsilon1: complex,
92 : epsilon2: complex,
93 : *,
94 : only_real_part: bool = False,
95 : ) -> NDArray:
96 : """Assemble the total scattering Green tensor.
97 :
98 : Args:
99 : pos1: Position vector of atom A (m)
100 : pos2: Position vector of atom B (m)
101 : z1: z-coordinate of the first surface (m)
102 : z2: z-coordinate of the second surface (m)
103 : omega: Angular frequency (i.e. 2*pi*f) in 1/s
104 : epsilon0: Electric permittivity of the medium between the two surfaces (dimensionless, complex)
105 : epsilon1: Electric permittivity of the upper medium (dimensionless, complex)
106 : epsilon2: Electric permittivity of the lower medium (dimensionless, complex)
107 : only_real_part: If True, only the real part of the Green tensor is calculated (default: False)
108 :
109 : Returns: The 3x3 Scattering Green Tensor (general complex values) (1/m)
110 :
111 : """
112 1 : if z1 > z2:
113 : # Ensure z1 is the lower surface and z2 is the upper surface
114 1 : z1, z2 = z2, z1
115 1 : epsilon1, epsilon2 = epsilon2, epsilon1
116 1 : if not (z1 < pos1[2] < z2 and z1 < pos2[2] < z2):
117 0 : raise ValueError("Both atoms must be located between the two surfaces (i.e. z1 < z_atom < z2).")
118 :
119 1 : distance = pos2 - pos1
120 1 : height = abs(z1 - z2)
121 :
122 1 : rho = np.sqrt(distance[0] ** 2 + distance[1] ** 2)
123 1 : phi = np.atan2(distance[1], distance[0]) if rho != 0 else 0
124 :
125 1 : z_ges = pos1[2] + pos2[2] - 2 * min(z1, z2)
126 1 : z_diff = pos1[2] - pos2[2]
127 :
128 1 : gt_scattered = np.zeros((3, 3), dtype=complex)
129 1 : for i, ix in enumerate(["x", "y", "z"]):
130 1 : for j, jx in enumerate(["x", "y", "z"]):
131 1 : entry: Entries = ix + jx # type: ignore [assignment]
132 1 : g_ij_elliptic = elliptic_integral(
133 : omega,
134 : height,
135 : rho,
136 : phi,
137 : epsilon0,
138 : epsilon1,
139 : epsilon2,
140 : z_ges,
141 : z_diff,
142 : entry,
143 : only_real_part=only_real_part,
144 : )
145 1 : g_ij_real = real_axis_integral(
146 : omega,
147 : height,
148 : rho,
149 : phi,
150 : epsilon0,
151 : epsilon1,
152 : epsilon2,
153 : z_ges,
154 : z_diff,
155 : entry,
156 : only_real_part=only_real_part,
157 : )
158 : # prefactor see comment in dynamic_green_tensor_homogeneous
159 1 : prefactor = -1 / (epsilon0 * const.epsilon_0 * const.hbar)
160 1 : value = prefactor * (g_ij_elliptic + g_ij_real)
161 1 : gt_scattered[i][j] = value
162 :
163 1 : return gt_scattered
164 :
165 :
166 1 : @njit(cache=True)
167 1 : def branch(epsilon: complex, k: float, k_rho: complex) -> complex:
168 : """Calculate the perpendicular wave vector component with positive imaginary part.
169 :
170 : Args:
171 : epsilon: Electric permittivity of the medium (dimensionless, complex)
172 : k: Magnitude of wave number in vacuum (1/m)
173 : k_rho: In-plane wave vector component (1/m)
174 :
175 : Returns: The perpendicular wave vector component (1/m)
176 :
177 : """
178 0 : return np.sqrt(epsilon * k**2 - k_rho**2 + 0j) # type: ignore [no-any-return]
179 :
180 :
181 1 : """The following functions are used from Appendix B of the paper:
182 : "Modified dipole-dipole interaction and dissipation in an atomic ensemble near surfaces"
183 : https://doi.org/10.1103/PhysRevA.97.053841
184 : and are needed to calculate the scattering Green Tensor for two atoms between two planar surfaces.
185 : We try to keep the same notation (in particular the same function names) as in the paper.
186 : """
187 :
188 :
189 1 : @njit(cache=True)
190 1 : def rs(kz: complex, k1z: complex) -> complex:
191 : """Calculate the Fresnel reflection coefficient for s-polarized light.
192 :
193 : Args:
194 : kz: Perpendicular wave vector component in the medium between the two surfaces (1/m)
195 : k1z: Perpendicular wave vector component in the upper or lower medium (1/m)
196 :
197 : Returns: The Fresnel reflection coefficient for s-polarized light (dimensionless, complex)
198 :
199 : """
200 0 : if kz == 0 and k1z == 0:
201 0 : return 0.0
202 0 : return (kz - k1z) / (kz + k1z)
203 :
204 :
205 1 : @njit(cache=True)
206 1 : def rp(kz: complex, k1z: complex, epsilon: complex) -> complex:
207 : """Calculate the Fresnel reflection coefficient for p-polarized light.
208 :
209 : Args:
210 : kz: Perpendicular wave vector component in the medium between the two surfaces (1/m)
211 : k1z: Perpendicular wave vector component in the upper or lower medium (1/m)
212 : epsilon: Electric permittivity of the upper or lower medium (dimensionless, complex)
213 :
214 : Returns: The Fresnel reflection coefficient for p-polarized light (dimensionless, complex)
215 :
216 : """
217 0 : if kz == 0 and k1z == 0:
218 0 : return 0.0
219 0 : return (epsilon * kz - k1z) / (epsilon * kz + k1z)
220 :
221 :
222 1 : @njit(cache=True)
223 1 : def D(r_plus: complex, r_minus: complex, kz: complex, h: float) -> complex:
224 : """Calculate the denominator term D used in the scattering Green Tensor matrix elements.
225 :
226 : Args:
227 : r_plus: Fresnel reflection coefficient for the upper surface (dimensionless, complex)
228 : r_minus: Fresnel reflection coefficient for the lower surface (dimensionless, complex)
229 : kz: Perpendicular wave vector component in the medium between the two surfaces (1/m)
230 : h: Distance between the two surfaces (m)
231 :
232 : Returns: The value of the denominator term D (dimensionless, complex)
233 :
234 : """
235 0 : return 1 - r_plus * r_minus * np.exp(2j * kz * h) # type: ignore [no-any-return]
236 :
237 :
238 1 : @njit(cache=True)
239 1 : def A_plus(r_plus: complex, r_minus: complex, kz: complex, h: float, z_ges: float, z_diff: float) -> complex:
240 : """Calculate the numerator term A_plus used in the scattering Green Tensor matrix elements.
241 :
242 : Args:
243 : r_plus: Fresnel reflection coefficient for the upper surface (dimensionless, complex)
244 : r_minus: Fresnel reflection coefficient for the lower surface (dimensionless, complex)
245 : kz: Perpendicular wave vector component in the medium between the two surfaces (1/m)
246 : h: Distance between the two surfaces (m)
247 : z_ges: Total z-coordinate (m)
248 : z_diff: Difference in z-coordinates (m)
249 :
250 : Returns: The value of the numerator term A_plus (dimensionless, complex)
251 :
252 : """
253 0 : return ( # type: ignore [no-any-return]
254 : r_minus * np.exp(1j * kz * (z_ges - h))
255 : + r_plus * np.exp(-1j * kz * (z_ges - h))
256 : + 2 * r_plus * r_minus * np.cos(kz * z_diff) * np.exp(1j * kz * h)
257 : ) / D(r_plus, r_minus, kz, h)
258 :
259 :
260 1 : @njit(cache=True)
261 1 : def A_minus(r_plus: complex, r_minus: complex, kz: complex, h: float, z_ges: float, z_diff: float) -> complex:
262 : """Calculate the numerator term A_minus used in the scattering Green Tensor matrix elements.
263 :
264 : Args:
265 : r_plus: Fresnel reflection coefficient for the upper surface (dimensionless, complex)
266 : r_minus: Fresnel reflection coefficient for the lower surface (dimensionless, complex)
267 : kz: Perpendicular wave vector component in the medium between the two surfaces (1/m)
268 : h: Distance between the two surfaces (m)
269 : z_ges: Total z-coordinate (m)
270 : z_diff: Difference in z-coordinates (m)
271 :
272 : Returns: The value of the numerator term A_minus (dimensionless, complex)
273 :
274 : """
275 0 : return ( # type: ignore [no-any-return]
276 : r_minus * np.exp(1j * kz * (z_ges - h))
277 : + r_plus * np.exp(-1j * kz * (z_ges - h))
278 : - 2 * r_plus * r_minus * np.cos(kz * z_diff) * np.exp(1j * kz * h)
279 : ) / D(r_plus, r_minus, kz, h)
280 :
281 :
282 1 : @njit(cache=True)
283 1 : def B_plus(r_plus: complex, r_minus: complex, kz: complex, h: float, z_ges: float, z_diff: float) -> complex:
284 : """Calculate the numerator term B_plus used in the scattering Green Tensor matrix elements.
285 :
286 : Args:
287 : r_plus: Fresnel reflection coefficient for the upper surface (dimensionless, complex)
288 : r_minus: Fresnel reflection coefficient for the lower surface (dimensionless, complex)
289 : kz: Perpendicular wave vector component in the medium between the two surfaces (1/m)
290 : h: Distance between the two surfaces (m)
291 : z_ges: Total z-coordinate (m)
292 : z_diff: Difference in z-coordinates (m)
293 :
294 : Returns: The value of the numerator term B_plus (dimensionless, complex)
295 :
296 : """
297 0 : return ( # type: ignore [no-any-return]
298 : r_minus * np.exp(1j * kz * (z_ges - h))
299 : - r_plus * np.exp(-1j * kz * (z_ges - h))
300 : + 2j * r_plus * r_minus * np.sin(kz * z_diff) * np.exp(1j * kz * h)
301 : ) / D(r_plus, r_minus, kz, h)
302 :
303 :
304 1 : @njit(cache=True)
305 1 : def B_minus(r_plus: complex, r_minus: complex, kz: complex, h: float, z_ges: float, z_diff: float) -> complex:
306 : """Calculate the numerator term B_minus used in the scattering Green Tensor matrix elements.
307 :
308 : Args:
309 : r_plus: Fresnel reflection coefficient for the upper surface (dimensionless, complex)
310 : r_minus: Fresnel reflection coefficient for the lower surface (dimensionless, complex)
311 : kz: Perpendicular wave vector component in the medium between the two surfaces (1/m)
312 : h: Distance between the two surfaces (m)
313 : z_ges: Total z-coordinate (m)
314 : z_diff: Difference in z-coordinates (m)
315 :
316 : Returns: The value of the numerator term B_minus (dimensionless, complex)
317 :
318 : """
319 0 : return ( # type: ignore [no-any-return]
320 : r_minus * np.exp(1j * kz * (z_ges - h))
321 : + r_plus * np.exp(-1j * kz * (z_ges - h))
322 : - 2j * r_plus * r_minus * np.sin(kz * z_diff) * np.exp(1j * kz * h)
323 : ) / D(r_plus, r_minus, kz, h)
324 :
325 :
326 1 : def Gs(
327 : kz: complex,
328 : h: float,
329 : k_rho: complex,
330 : rho: float,
331 : phi: float,
332 : rs_plus: complex,
333 : rs_minus: complex,
334 : z_ges: float,
335 : z_diff: float,
336 : entry: Entries,
337 : ) -> complex:
338 : """Calculate the Gs part of the scattering Green Tensor."""
339 1 : if entry in ["xz", "yz", "zx", "zy", "zz"]:
340 1 : return 0
341 1 : As_plus = A_plus(rs_plus, rs_minus, kz, h, z_ges, z_diff)
342 1 : J2 = cached_bessel_function_2(k_rho * rho)
343 1 : if entry in ["xy", "yx"]:
344 1 : return -As_plus / 2 * J2 * math.sin(2 * phi)
345 1 : J0 = cached_bessel_function_0(k_rho * rho)
346 1 : if entry == "xx":
347 1 : return As_plus / 2 * (J0 + J2 * math.cos(2 * phi))
348 1 : if entry == "yy":
349 1 : return As_plus / 2 * (J0 - J2 * math.cos(2 * phi))
350 :
351 0 : raise ValueError(f"Invalid entry '{entry}' for Gs function.")
352 :
353 :
354 1 : def Gp( # noqa: PLR0911
355 : kz: complex,
356 : h: float,
357 : k_rho: complex,
358 : rho: float,
359 : phi: float,
360 : rp_plus: complex,
361 : rp_minus: complex,
362 : z_ges: float,
363 : z_diff: float,
364 : entry: Entries,
365 : ) -> complex:
366 : """Calculate the Gp part of the scattering Green Tensor."""
367 1 : if entry in ["xz", "zx", "yz", "zy"]:
368 1 : J1 = cached_bessel_function_1(k_rho * rho)
369 1 : if entry == "zx":
370 1 : Bp_minus = B_minus(rp_plus, rp_minus, kz, h, z_ges, z_diff)
371 1 : return -1j * (k_rho / kz) * Bp_minus * J1 * math.cos(phi)
372 1 : if entry == "xz":
373 1 : Bp_plus = B_plus(rp_plus, rp_minus, kz, h, z_ges, z_diff)
374 1 : return 1j * (k_rho / kz) * Bp_plus * J1 * math.cos(phi)
375 1 : if entry == "yz":
376 1 : Bp_plus = B_plus(rp_plus, rp_minus, kz, h, z_ges, z_diff)
377 1 : return 1j * (k_rho / kz) * Bp_plus * J1 * math.sin(phi)
378 1 : if entry == "zy":
379 1 : Bp_minus = B_minus(rp_plus, rp_minus, kz, h, z_ges, z_diff)
380 1 : return -1j * (k_rho / kz) * Bp_minus * J1 * math.sin(phi)
381 :
382 1 : if entry == "zz":
383 1 : J0 = cached_bessel_function_0(k_rho * rho)
384 1 : Ap_plus = A_plus(rp_plus, rp_minus, kz, h, z_ges, z_diff)
385 1 : return -(k_rho**2 / kz**2) * Ap_plus * J0
386 :
387 1 : J2 = cached_bessel_function_2(k_rho * rho)
388 1 : Ap_minus = A_minus(rp_plus, rp_minus, kz, h, z_ges, z_diff)
389 1 : if entry in ["xy", "yx"]:
390 1 : return Ap_minus / 2 * J2 * math.sin(2 * phi)
391 :
392 1 : J0 = cached_bessel_function_0(k_rho * rho)
393 1 : if entry == "xx":
394 1 : return Ap_minus / 2 * (J0 - J2 * math.cos(2 * phi))
395 1 : if entry == "yy":
396 1 : return Ap_minus / 2 * (J0 + J2 * math.cos(2 * phi))
397 :
398 0 : raise ValueError(f"Invalid entry '{entry}' for Gp function.")
399 :
400 :
401 1 : """ The integrals for the scattering Green Tensor are evaluated in two parts:
402 : The first part is along an elliptical path from 0 to 2*k_maj in the complex plane to avoid singularities,
403 : and the second part is along the real axis from 2*k_maj to an upper limit.
404 :
405 : The methods for evaluating these integrals are explained in the papers:
406 : - "Accurate and efficient computation of the Green's tensor for stratified media"
407 : section III.A (https://doi.org/10.1103/PhysRevE.62.5797)
408 : - "Challenges in Computational Electromagnetics: Analysis and Optimization of Planar Multilayered Structures"
409 : section 2.4.1 (https://doi.org/10.5075/epfl-thesis-5122)
410 : """
411 :
412 :
413 1 : @njit(cache=True)
414 1 : def integrand_ellipse_partial(
415 : t: complex,
416 : k_maj: float,
417 : k_min: float,
418 : k0: float,
419 : epsilon0: float,
420 : epsilon1: complex,
421 : epsilon2: complex,
422 : h: float,
423 : ) -> tuple[complex, complex, complex, complex, complex, complex, complex]:
424 : # elliptical path, substitution
425 0 : k_rho = k_maj * (1 + np.cos(t)) - 1j * k_min * np.sin(t)
426 0 : dk_rho = -k_maj * np.sin(t) - 1j * k_min * np.cos(t)
427 :
428 : # Wave vector components
429 0 : kz = branch(epsilon0, k0, k_rho)
430 0 : k1z = branch(epsilon1, k0, k_rho)
431 0 : k2z = branch(epsilon2, k0, k_rho)
432 :
433 0 : rs_plus = rs(kz, k1z)
434 0 : rs_minus = rs(kz, k2z)
435 0 : rp_plus = rp(kz, k1z, epsilon1)
436 0 : rp_minus = rp(kz, k2z, epsilon2)
437 :
438 0 : if k_rho == 0 and dk_rho == 0: # noqa: SIM108
439 0 : prefactor = 0.0
440 : else:
441 0 : prefactor = 1j / (4 * np.pi) * (k_rho / kz) * np.exp(1j * kz * h) * dk_rho
442 0 : return k_rho, kz, rs_plus, rs_minus, rp_plus, rp_minus, prefactor
443 :
444 :
445 1 : def integrand_ellipse(
446 : t: complex,
447 : k_maj: float,
448 : k_min: float,
449 : k0: float,
450 : epsilon0: float,
451 : epsilon1: complex,
452 : epsilon2: complex,
453 : h: float,
454 : rho: float,
455 : phi: float,
456 : z_ges: float,
457 : z_diff: float,
458 : entry: Entries,
459 : real_or_imag: str,
460 : ) -> complex:
461 1 : k_rho, kz, rs_plus, rs_minus, rp_plus, rp_minus, prefactor = integrand_ellipse_partial(
462 : t, k_maj, k_min, k0, epsilon0, epsilon1, epsilon2, h
463 : )
464 1 : if k0 == 0 and kz == 0:
465 1 : return 0.0
466 :
467 0 : gs = Gs(kz, h, k_rho, rho, phi, rs_plus, rs_minus, z_ges, z_diff, entry)
468 0 : gp = Gp(kz, h, k_rho, rho, phi, rp_plus, rp_minus, z_ges, z_diff, entry)
469 :
470 0 : integrand = prefactor * (k0**2 * gs - kz**2 * gp)
471 0 : if real_or_imag == "real":
472 0 : return np.real(integrand)
473 0 : if real_or_imag == "imag":
474 0 : return np.imag(integrand)
475 0 : raise ValueError("real_or_imag must be 'real' or 'imag'")
476 :
477 :
478 1 : def elliptic_integral(
479 : omega: float,
480 : h: float,
481 : rho: float,
482 : phi: float,
483 : epsilon0: complex,
484 : epsilon1: complex,
485 : epsilon2: complex,
486 : z_ges: float,
487 : z_diff: float,
488 : entry: Entries,
489 : *,
490 : only_real_part: bool = False,
491 : ) -> complex:
492 : """Evaluate the elliptic part of the integral.
493 :
494 : Args:
495 : omega: Angular frequency (i.e. 2*pi*f) in 1/s
496 : h: Distance between the two surfaces in meters
497 : rho: In-plane distance between the two atoms in meters
498 : phi: Angle between the in-plane distance vector and the x-axis in radians
499 : epsilon0: Electric permittivity of the medium between the two surfaces (dimensionless, complex)
500 : epsilon1: Electric permittivity of the upper medium (dimensionless, complex)
501 : epsilon2: Electric permittivity of the lower medium (dimensionless, complex)
502 : z_ges: Sum of the z-positions of the two atoms in meters
503 : z_diff: Difference of the z-positions of the two atoms in meters
504 : entry: Entry of the Green tensor to calculate
505 : only_real_part: If True, only the real part of the integral is calculated (default: False)
506 :
507 : Returns: The value of the integral along the elliptical path as a complex number (1/m)
508 :
509 : """
510 1 : k_vac = omega / const.c # magnitude of wave vector in vacuum
511 1 : k0 = k_vac * np.sqrt(epsilon0)
512 :
513 : # Elliptical path in complex plane to avoid singularities (Integral from 0 to 2k_maj)
514 1 : k1 = k_vac * np.sqrt(epsilon1)
515 1 : k2 = k_vac * np.sqrt(epsilon2)
516 1 : kl_max = max(np.real(k0), np.real(k1), np.real(k2))
517 :
518 1 : k_maj = (kl_max + k_vac) / 2 # major axis of ellipse
519 1 : k_min = min(k_vac, 1 / rho) if rho != 0 else k_vac
520 :
521 1 : args = (k_maj, k_min, k0, epsilon0, epsilon1, epsilon2, h, rho, phi, z_ges, z_diff, entry)
522 :
523 1 : real_ellipse, _ = quad(integrand_ellipse, np.pi, 0, args=(*args, "real"), epsrel=1e-9, limit=1000) # type: ignore [arg-type]
524 1 : if only_real_part:
525 1 : return real_ellipse
526 0 : imag_ellipse, _ = quad(integrand_ellipse, np.pi, 0, args=(*args, "imag"), epsrel=1e-9, limit=1000) # type: ignore [arg-type]
527 0 : return real_ellipse + 1j * imag_ellipse
528 :
529 :
530 1 : def integrand_real(
531 : k_rho: complex,
532 : k0: float,
533 : epsilon0: float,
534 : epsilon1: complex,
535 : epsilon2: complex,
536 : h: float,
537 : rho: float,
538 : phi: float,
539 : z_ges: float,
540 : z_diff: float,
541 : entry: Entries,
542 : real_or_imag: str,
543 : ) -> complex:
544 1 : kz = branch(epsilon0, k0, k_rho)
545 1 : if kz == 0 and k0 == 0:
546 0 : return 0
547 :
548 1 : k1z = branch(epsilon1, k0, k_rho)
549 1 : k2z = branch(epsilon2, k0, k_rho)
550 :
551 1 : rs_plus = rs(kz, k1z)
552 1 : rs_minus = rs(kz, k2z)
553 1 : rp_plus = rp(kz, k1z, epsilon1)
554 1 : rp_minus = rp(kz, k2z, epsilon2)
555 :
556 : # Integrand
557 1 : integrand = (
558 : 1j
559 : / (4 * np.pi)
560 : * (
561 : k0**2 * Gs(kz, h, k_rho, rho, phi, rs_plus, rs_minus, z_ges, z_diff, entry)
562 : - kz**2 * Gp(kz, h, k_rho, rho, phi, rp_plus, rp_minus, z_ges, z_diff, entry)
563 : )
564 : * (k_rho / kz)
565 : * np.exp(1j * kz * h)
566 : )
567 1 : if real_or_imag == "real":
568 1 : return np.real(integrand) # type: ignore [no-any-return]
569 0 : if real_or_imag == "imag":
570 0 : return np.imag(integrand) # type: ignore [no-any-return]
571 0 : raise ValueError("real_or_imag must be 'real' or 'imag'")
572 :
573 :
574 1 : def real_axis_integral(
575 : omega: float,
576 : h: float,
577 : rho: float,
578 : phi: float,
579 : epsilon0: complex,
580 : epsilon1: complex,
581 : epsilon2: complex,
582 : z_ges: float,
583 : z_diff: float,
584 : entry: Entries,
585 : *,
586 : only_real_part: bool = False,
587 : ) -> complex:
588 : """Evaluate the real axis part of the integral.
589 :
590 : Args:
591 : omega: Angular frequency (i.e. 2*pi*f) in 1/s
592 : h: Distance between the two surfaces in meters
593 : rho: In-plane distance between the two atoms in meters
594 : phi: Angle between the in-plane distance vector and the x-axis in radians
595 : epsilon0: Electric permittivity of the medium between the two surfaces (dimensionless, complex)
596 : epsilon1: Electric permittivity of the upper medium (dimensionless, complex)
597 : epsilon2: Electric permittivity of the lower medium (dimensionless, complex)
598 : z_ges: Sum of the z-positions of the two atoms in meters
599 : z_diff: Difference of the z-positions of the two atoms in meters
600 : entry: Entry of the Green tensor to calculate
601 : upper_limit: Upper limit for the real axis integral (1/m)
602 : only_real_part: If True, only the real part of the integral is calculated (default: False)
603 :
604 : Returns: The value of the integral along the real axis as a complex number (1/m)
605 :
606 : """
607 1 : k_vac = omega / const.c # magnitude of wave vector in vacuum
608 1 : k0 = k_vac * np.sqrt(epsilon0)
609 1 : k1 = k_vac * np.sqrt(epsilon1)
610 1 : k2 = k_vac * np.sqrt(epsilon2)
611 :
612 1 : kl_max = max(np.real(k0), np.real(k1), np.real(k2))
613 1 : k_maj = (kl_max + k_vac) / 2
614 :
615 1 : args = (k0, epsilon0, epsilon1, epsilon2, h, rho, phi, z_ges, z_diff, entry)
616 :
617 : # Estimate the upper limit for the real axis integral
618 1 : upper_limit = np.sqrt((745 / h) ** 2 + 1)
619 :
620 1 : real_real, _ = quad(
621 : integrand_real, # type: ignore [arg-type]
622 : 2 * k_maj,
623 : upper_limit,
624 : args=(*args, "real"), # type: ignore [arg-type]
625 : limit=1000,
626 : epsrel=1e-9,
627 : )
628 1 : if only_real_part:
629 1 : return real_real
630 0 : imag_real, _ = quad(
631 : integrand_real, # type: ignore [arg-type]
632 : 2 * k_maj,
633 : upper_limit,
634 : args=(*args, "imag"), # type: ignore [arg-type]
635 : limit=1000,
636 : epsrel=1e-9,
637 : )
638 0 : return real_real + 1j * imag_real
|
