#
#
# A medium in which there holds an equality
# $$\langle\langle I, II \rangle\rangle = \langle\langle II, I \rangle\rangle\rangle$$
# is called the Lorentz reciprocal medium. The Lorentz lemma refers to the case of electric currents only, and states that for any sufficiently good surface $S$ bounding some volume $V$ in a Lorentz reciprocal medium, there holds
# $$\oint_{S}\left(\boldsymbol{E}_{1}\times\boldsymbol{H}_{2} - \boldsymbol{E}_{2}\times\boldsymbol{H}_{1}\right) \cdot d\boldsymbol{S} = \int_{V} \left(\boldsymbol{E}_{2}\cdot\boldsymbol{J}_{1}- \boldsymbol{E}_{1}\cdot\boldsymbol{J}_{2}\right)dV$$
# ## Wave equation
#
# The mutual substitution of the Faraday's law and the circulation theorem for the magnetic field in the case of anisotropic media leads to the Helmholtz vector equations for the electric and magnetic fields. For the electric field, this equation writes as
# $$\nabla\times\hat{\mu}^{-1}\nabla\times\boldsymbol{E} - \omega^{2}\hat{\varepsilon}\boldsymbol{E} = i\omega\boldsymbol{J}_e - \nabla\times\hat{\mu}^{-1}\boldsymbol{J}_m$$
#
# In a homogeneous isotropic medium with constant dielectric constant $\varepsilon$ and magnetic susceptibility $\mu$ the equation is simplified
# $$\nabla\times\nabla\times\boldsymbol{E} - k^{2}\boldsymbol{E} = i\omega\boldsymbol{J}_e - \nabla\times\boldsymbol{J}_m$$
# where $k=\omega\sqrt{\varepsilon\mu}$ is the wavenumber. In the absence of free charges and their currents the equation reduces to the scalar Helmholtz equation
# $$\nabla\boldsymbol{E} = 0 \Rightarrow \Delta\boldsymbol{E} + k^2\boldsymbol{E} = 0 \Rightarrow (\Delta + k^2)\psi(\boldsymbol{r}) = 0$$
# with the Laplace operator $\Delta$.
# ## Problems in electrodynamics
#
# Electrodynamics problems can be divided into direct and inverse problems. In a direct problem, for given material tensors and given currents in space, the fields are calculated, i.e. the boundary value problem or the scattering problem is solved. For inverse problems for the given target parameters of the fields, the medium/structure is selected to best approximate these target parameters. As a rule, inverse problems are incorrect, i.e., they violate one or more of the properties of correctness: the condition of existence, uniqueness, and stability of solutions with respect to small perturbations.
# ## Eigen waves in a homogeneous isotropic medium
#
# For numerical solutions, analytical solutions of the eigenvalue problem for the Helmholtz equation in different coordinate systems are important. Let us write these solutions for a homogeneous isotropic medium with material constants $\varepsilon$ and $\mu$. In such a medium, the eigen solutions of the vector Helmholtz equation are transverse waves admitting a decomposition of the field into two linearly independent polarizations. In a general form, consider the vector and scalar equations:
# $$\nabla\times\nabla\times\boldsymbol{F} - k^{2}\boldsymbol{F} = 0, \;\; \left(\nabla^{2}+k^{2}\right)\psi\left(\boldsymbol{r}\right) = 0$$
# If the solutions of the scalar equation $\psi$ are known, the two sets of solutions to the vector equation are obtained from the relations
# \begin{align*}
# \boldsymbol{M} &= \nabla\times\boldsymbol{c}\psi(\boldsymbol{r}) \\
# \boldsymbol{N} &= \frac{1}{k} \nabla\times\boldsymbol{M}\left(\boldsymbol{r}\right)
# \end{align*}
# These vector functions are solenoidal and have zero divergence. The basis in three-dimensional space can be constructed by adding the irrotational field
# $$\boldsymbol{L} = \nabla\psi(\boldsymbol{r})$$
# This function turns out to be necessary to decompose the field in the source region.
# ### Cartesian coordinates and plane waves
#
# In Cartesian coordinates the solution of the scalar equation is the exponential function
# $$\psi(\boldsymbol{r}) = a \exp(i\boldsymbol{k}\cdot\boldsymbol{r})$$
# with wavevector $\boldsymbol{k}$, which components meet the dispersion equation
# $$k_{x}^{2}+k_{y}^{2}+k_{z}^{2} = k^{2}$$
# Often the propagation of plane waves is considered relative to some selected axis, for example $Z$, and the dispersion relation is written in the form of
# $$k_z = \sqrt{k^{2} - k_{x}^{2} - k_{y}^{2}}$$
# where the square root branch in the case of complex material constants is chosen from the physical requirement of exponential wave attenuation at infinity, $\Im m (k_z) \geq 0$. In this case, each projection of the wavevector on the plane $XY$ corresponds to two wavevectors describing the propagation of waves in the positive and negative directions relative to the axis $Z$:
# $$\boldsymbol{k}^{\pm} = (k_x,\;k_y\;,\pm k_z)^T$$
#
# For vector waves let us take $\boldsymbol{c} \equiv \hat{\boldsymbol{z}}$. Then,
# \begin{align*}
# \boldsymbol{M}\left(\boldsymbol{k},\boldsymbol{r}\right) &= i\boldsymbol{k}\times\hat{\boldsymbol{z}}e^{i\boldsymbol{k}\cdot\boldsymbol{r}} = i\kappa \hat{\boldsymbol{e}}^{s\pm}_{\boldsymbol{k}} e^{i\boldsymbol{k}\cdot\boldsymbol{r}} \\
# \boldsymbol{N}\left(\boldsymbol{k},\boldsymbol{r}\right) &= -\dfrac{\kappa}{k}\boldsymbol{k}\times\boldsymbol{k}\times\hat{\boldsymbol{z}}e^{i\boldsymbol{k}\cdot\boldsymbol{r}} = -\kappa \hat{\boldsymbol{e}}^{p\pm}_{\boldsymbol{k}} e^{i\boldsymbol{k}\cdot\boldsymbol{r}} \\
# \boldsymbol{L}\left(\boldsymbol{k},\boldsymbol{r}\right) &= i\boldsymbol{k}e^{i\boldsymbol{k}\cdot\boldsymbol{r}}
# \end{align*}
# Here the unit vectors for the TE ('s') and TM ('p') polarizations are:
# \begin{align*}
# \hat{\boldsymbol{e}}^{s\pm}_{\boldsymbol{k}} = \dfrac{\boldsymbol{k}^{\pm}\times\hat{\boldsymbol{e}}_{z}}{\left|\boldsymbol{k}^{\pm}\times\hat{\boldsymbol{e}}_{z}\right|} &= \dfrac{k_{y}}{\kappa}\hat{\boldsymbol{e}}_{x}-\dfrac{k_{x}}{\kappa}\hat{\boldsymbol{e}}_{y} \\
# \hat{\boldsymbol{e}}^{p\pm}_{\boldsymbol{k}} = \dfrac{\boldsymbol{k}^{\pm}\times\left(\boldsymbol{k}^{\pm}\times\hat{\boldsymbol{e}}_{z}\right)}{\left|\boldsymbol{k}^{\pm}\times\left(\boldsymbol{k}^{\pm}\times\hat{\boldsymbol{e}}_{z}\right)\right|} &= \pm\dfrac{k_{x}k_{z}}{\kappa k}\hat{\boldsymbol{e}}_{x}\pm\dfrac{k_{y}k_{z}}{\kappa k}\hat{\boldsymbol{e}}_{y}-\dfrac{\kappa}{k}\hat{\boldsymbol{e}}_{z}
# \end{align*}
# and the $XY$ projection length of the wavevector is: $\varkappa = \sqrt{k_x^2 + k_y^2}$.
#
# The fields are decomposed into plane waves using the orthogonality relations:
# \begin{equation*}
# \iiint_{-\infty}^{\infty}d^{3}\boldsymbol{r}\psi\left(\boldsymbol{k},\boldsymbol{r}\right)\psi\left(-\boldsymbol{k}',\boldsymbol{r}\right)=\left(2\pi\right)^{3}\delta\left(\boldsymbol{k}-\boldsymbol{k}'\right)
# \end{equation*}
#
# \begin{equation*}
# \iiint_{-\infty}^{\infty}d^{3}\boldsymbol{r}\boldsymbol{M}\left(\boldsymbol{k},\boldsymbol{r}\right)\cdot\boldsymbol{M}\left(-\boldsymbol{k}',\boldsymbol{r}\right) = \iiint_{-\infty}^{\infty}d^{3}\boldsymbol{k}\boldsymbol{N}\left(\boldsymbol{k},\boldsymbol{r}\right)\cdot\boldsymbol{N}\left(-\boldsymbol{k}',\boldsymbol{r}\right)=\left(2\pi\right)^{3}\varkappa^{2}\delta\left(\boldsymbol{k}-\boldsymbol{k}'\right)
# \end{equation*}
#
# \begin{equation*}
# \iiint_{-\infty}^{\infty}d^{3}\boldsymbol{r}\boldsymbol{L}\left(\boldsymbol{k},\boldsymbol{r}\right)\cdot\boldsymbol{L}\left(-\boldsymbol{k}',\boldsymbol{r}\right)=\left(2\pi\right)^{3}k^{2}\delta\left(\boldsymbol{k}-\boldsymbol{k}'\right)
# \end{equation*}
#
# ### Spherical coordinates and spherical waves
#
# In spherical coordinates $(r,\theta,\phi)$ the scalar Helmholtz equation has the following form
# \begin{equation*}
# \left(\dfrac{1}{r^{2}}\dfrac{\partial}{\partial r}r^{2}\dfrac{\partial}{\partial r}+\dfrac{1}{r^{2}\sin\theta}\dfrac{\partial}{\partial\theta}\sin\theta\dfrac{\partial}{\partial\theta}+\dfrac{1}{r^{2}\sin^{2}\theta}\dfrac{\partial^{2}}{\partial\phi^{2}}+k^{2}\right)\psi=0
# \end{equation*}
# Its solutions are written through the spherical Bessel functions $z_n(x)$ and the associated Legendre polynomials $P_n^m(\cos\theta)$:
# \begin{equation*}
# \psi_{mn} = \gamma_{mn} z_{n}\left(kr\right) P_{n}^{m}\left(\theta\right) e^{im\phi} = z_{n}\left(kr\right) Y_{mn} (\theta, \phi)
# \end{equation*}
# The indices are $0\leq n\leq\infty$, $-n\leq m\leq n$. All solutions are superpositions of two types of solutions that are usually chosen from four types of spherical functions - Bessel functions of the first kind $j_n(x)$, Neumann functions $y_n(x)$, and Hankel functions of the first and second kind $h_n^{(1,2)}(x)=j_n(x)\pm iy_n(x)$ (analogously to the two types of plane waves propagating in the positive and negative directions with respect to the $Z$ axis). Typically, one chooses a pair of the regular at the origin $j_n(x)$ and the $h_n^{(1)}(x)$, which is a divergent spherical wave for large arguments. For the function $z_n^{\sigma}$, we introduce an upper index $\sigma=1,3$ and assume $z^1_n\equiv j_n$, $z^3_n\equiv h^{(1)}_n$. The solutions are often normalized using the multiplier
# \begin{equation*}
# \gamma_{mn} = \sqrt{ \dfrac{\left(2n+1\right)}{4\pi} \dfrac{\left(n-m\right)!}{\left(n+m\right)!} }
# \end{equation*}
#
# To construct vector spherical waves one chooses $\boldsymbol{c}\equiv\boldsymbol{r}$, then
# \begin{align*}
# \boldsymbol{M}_{nm}(k\boldsymbol{r}) = \nabla\times\boldsymbol{r} z_{n}(kr) Y_{nm}(\theta,\phi) =& \frac{1}{\sqrt{2\pi n\left(n+1\right)}} z_{n}(kr) \left[ im\pi_{nm}(\theta)\hat{\boldsymbol{e}}_{\theta} - \tau_{nm}(\theta) \hat{\boldsymbol{e}}_{\phi} \right] e^{im\phi} \\
# \boldsymbol{N}_{nm}(k\boldsymbol{r}) = \dfrac{1}{k} \nabla\times\nabla\times\boldsymbol{r} z_{n}(kr) Y_{nm}(\theta,\phi) =& \dfrac{\sqrt{n\left(n+1\right)}}{\sqrt{2\pi}} \frac{z_{n}(kr)}{kr} P_{n}^{m}(\theta) e^{im\phi} \hat{\boldsymbol{e}}_{r} + \\
# &+ \frac{1}{\sqrt{2\pi n\left(n+1\right)}} \frac{\tilde{z}_{n}(kr)}{kr} \left[ \tau_{nm}(\theta) \hat{\boldsymbol{e}}_{\theta} + im\pi_{nm}(\theta) \hat{\boldsymbol{e}}_{\varphi} \right] e^{im\phi} \\
# \boldsymbol{L}_{nm}(k\boldsymbol{r}) = \dfrac{1}{k} \nabla j_{n}(kr) Y_{nm}(\theta,\phi) =& z_{n}'(kr) P_{n}^{m}(\theta) e^{im\phi} \hat{\boldsymbol{e}}_{r} + \\
# &+ \frac{z_{n}(kr)}{kr} \left[ \tau_{nm}(\theta) \hat{\boldsymbol{e}}_{\theta} + im\pi_{nm}(\theta)\hat{\boldsymbol{e}}_{\phi} \right] e^{im\phi}
# \end{align*}
# где $\tilde{z}_{n}\left(x\right)=d\left(xz_{n}\left(x\right)\right)/dx$, $\pi_{nm}(\theta) = P_{n}^{m}(\theta)/\sin\theta$, и $\tau_{nm}(\theta) = dP_{n}^{m}(\theta)/d\theta$.
#
# Orthogonality relations read
# \begin{equation*}
# \int_V d^3\boldsymbol{r} \psi^{\sigma}_{nm}(k\boldsymbol{r}) \psi^{\sigma}_{n',-m'}(k'\boldsymbol{r}) = \int_{V} d^3\boldsymbol{r} \boldsymbol{L}^{\sigma}_{nm}(k\boldsymbol{r}) \cdot \boldsymbol{L}^{\sigma}_{n',-m'}(k'\boldsymbol{r}) = \pi\delta_{mm'}\delta_{nn'}\dfrac{\delta\left(k-k'\right)}{2k^{2}}
# \end{equation*}
#
# \begin{equation*}
# \int_{V} d^3\boldsymbol{r} \boldsymbol{M}^{\sigma}_{nm}(k\boldsymbol{r}) \cdot \boldsymbol{M}^{\sigma}_{n',-m'}(k'\boldsymbol{r}) = \int_{V} d^3\boldsymbol{r} \boldsymbol{N}_{nm}(k\boldsymbol{r}) \cdot \boldsymbol{N}^{\sigma}_{n',-m'}(k'\boldsymbol{r}) = n(n+1)\pi\delta_{mm'}\delta_{nn'}\dfrac{\delta(k-k')}{2k^{2}}
# \end{equation*}
#
# Cylindrical waves can be derived analogously.
# References:
# 1. Orfanidis, J. Electromagnetic waves and antennas, Ch. 1-2, Rutgers University (2016)
# 2. W. C. Chew, Waves and Fields in Inhomogeneous Media, Ch 1, 7.2, IEEE Press (1995)
# 3. A. G. Sveshnikov, I. E. Mogilevskiy, Mathematical problems of the diffraction theory, Ch. 1, MSU (2012) (In Russian)
# 4. J. A. Kong, Theorems of bianisotropic media, Proceedings of the IEEE 60, 9, 1036-1046 (1972)
#
# Translated partially with DeepL.com (free version)
#