603 lines
22 KiB
Python
603 lines
22 KiB
Python
r"""
|
|
Operators and helper functions for cylindrical waveguides with unchanging cross-section.
|
|
|
|
Waveguide operator is derived according to 10.1364/OL.33.001848.
|
|
|
|
As in `waveguide_2d`, the propagation dependence is separated from the
|
|
transverse solve. Here the propagation coordinate is the bend angle `\theta`,
|
|
and the fields are assumed to have the form
|
|
|
|
$$
|
|
\vec{E}(r, y, \theta), \vec{H}(r, y, \theta) \propto e^{-\imath m \theta},
|
|
$$
|
|
|
|
where `m` is the angular wavenumber returned by `solve_mode(s)`. It is often
|
|
convenient to introduce the corresponding linear wavenumber
|
|
|
|
$$
|
|
\beta = \frac{m}{r_{\min}},
|
|
$$
|
|
|
|
so that the cylindrical problem resembles the straight-waveguide problem with
|
|
additional metric factors.
|
|
|
|
Those metric factors live on the staggered radial Yee grids. If the left edge of
|
|
the computational window is at `r = r_{\min}`, define the electric-grid and
|
|
magnetic-grid radial sample locations by
|
|
|
|
$$
|
|
\begin{aligned}
|
|
r_a(n) &= r_{\min} + \sum_{j \le n} \Delta r_{e, j}, \\
|
|
r_b\!\left(n + \tfrac{1}{2}\right) &= r_{\min} + \tfrac{1}{2}\Delta r_{e, n}
|
|
+ \sum_{j < n} \Delta r_{h, j},
|
|
\end{aligned}
|
|
$$
|
|
|
|
and from them the diagonal metric matrices
|
|
|
|
$$
|
|
\begin{aligned}
|
|
T_a &= \operatorname{diag}(r_a / r_{\min}), \\
|
|
T_b &= \operatorname{diag}(r_b / r_{\min}).
|
|
\end{aligned}
|
|
$$
|
|
|
|
With the same forward/backward derivative notation used in `waveguide_2d`, the
|
|
implementation treats the transverse electric eigenproblem as the canonical
|
|
cylindrical discretization. It reduces to `waveguide_2d.operator_e(...)` in the
|
|
large-radius limit `T_a, T_b \to I`. The eigenproblem implemented by
|
|
`cylindrical_operator(...)`:
|
|
|
|
$$
|
|
\beta^2
|
|
\begin{bmatrix} E_r \\ E_y \end{bmatrix}
|
|
=
|
|
\left(
|
|
\omega^2
|
|
\begin{bmatrix}
|
|
T_b^2 \mu_{yy} \epsilon_{xx} & 0 \\
|
|
0 & T_a^2 \mu_{xx} \epsilon_{yy}
|
|
\end{bmatrix}
|
|
+
|
|
\begin{bmatrix}
|
|
-T_b \mu_{yy} \hat{\partial}_y \\
|
|
T_a \mu_{xx} \hat{\partial}_x
|
|
\end{bmatrix}
|
|
T_b \mu_{zz}^{-1}
|
|
\begin{bmatrix}
|
|
-\tilde{\partial}_y & \tilde{\partial}_x
|
|
\end{bmatrix}
|
|
+
|
|
\begin{bmatrix}
|
|
\tilde{\partial}_x \\
|
|
\tilde{\partial}_y
|
|
\end{bmatrix}
|
|
T_a \epsilon_{zz}^{-1}
|
|
\begin{bmatrix}
|
|
\hat{\partial}_x T_b \epsilon_{xx} &
|
|
\hat{\partial}_y T_a \epsilon_{yy}
|
|
\end{bmatrix}
|
|
\right)
|
|
\begin{bmatrix} E_r \\ E_y \end{bmatrix}.
|
|
$$
|
|
|
|
The full fields reconstructed by `exy2e(...)` and `e2h(...)` use the matching
|
|
large-radius-compatible identities
|
|
|
|
$$
|
|
E_z =
|
|
\frac{1}{\imath \beta} T_a \epsilon_{zz}^{-1}
|
|
\begin{bmatrix}
|
|
\hat{\partial}_r T_b \epsilon_{rr} &
|
|
\hat{\partial}_y T_a \epsilon_{yy}
|
|
\end{bmatrix}
|
|
\begin{bmatrix} E_r \\ E_y \end{bmatrix},
|
|
$$
|
|
|
|
and
|
|
|
|
$$
|
|
\begin{bmatrix} H_r \\ H_y \\ H_z \end{bmatrix}
|
|
=
|
|
\frac{1}{-\imath \omega}\mu^{-1}
|
|
\begin{bmatrix}
|
|
0 & \imath\beta T_a^{-1} & \tilde{\partial}_y \\
|
|
-\imath\beta T_b^{-1} & 0 & -T_b^{-1}\tilde{\partial}_r T_a \\
|
|
-\tilde{\partial}_y & \tilde{\partial}_r & 0
|
|
\end{bmatrix}
|
|
\begin{bmatrix} E_r \\ E_y \\ E_z \end{bmatrix}.
|
|
$$
|
|
|
|
Since `\beta = m / r_{\min}`, the solver implemented in this file returns the
|
|
angular wavenumber `m`, while the operator itself is most naturally written in
|
|
terms of the linear quantity `\beta`. The helpers below reconstruct the full
|
|
field components from the solved transverse eigenvector and then normalize the
|
|
mode to unit forward power with the same discrete longitudinal Poynting inner
|
|
product used by `waveguide_2d`.
|
|
|
|
As in the straight-waveguide case, all functions here assume a 2D grid:
|
|
|
|
`dxes = [[[dr_e_0, dr_e_1, ...], [dy_e_0, ...]], [[dr_h_0, ...], [dy_h_0, ...]]]`.
|
|
"""
|
|
from typing import Any, cast
|
|
from collections.abc import Sequence
|
|
import logging
|
|
|
|
import numpy
|
|
from numpy.typing import NDArray, ArrayLike
|
|
from scipy import sparse
|
|
|
|
from ..fdmath import vec, unvec, dx_lists2_t, vcfdslice_t, vfdslice, vcfdslice, vcfdfield2
|
|
from ..fdmath.operators import deriv_forward, deriv_back
|
|
from ..eigensolvers import signed_eigensolve, rayleigh_quotient_iteration
|
|
from . import waveguide_2d
|
|
|
|
logger = logging.getLogger(__name__)
|
|
|
|
|
|
def cylindrical_operator(
|
|
omega: float,
|
|
dxes: dx_lists2_t,
|
|
epsilon: vfdslice,
|
|
rmin: float,
|
|
mu: vfdslice | None = None,
|
|
) -> sparse.sparray:
|
|
r"""
|
|
Cylindrical coordinate waveguide operator of the form
|
|
|
|
$$
|
|
(\omega^2 \begin{bmatrix} T_b T_b \mu_{yy} \epsilon_{xx} & 0 \\
|
|
0 & T_a T_a \mu_{xx} \epsilon_{yy} \end{bmatrix} +
|
|
\begin{bmatrix} -T_b \mu_{yy} \hat{\partial}_y \\
|
|
T_a \mu_{xx} \hat{\partial}_x \end{bmatrix} T_b \mu_{zz}^{-1}
|
|
\begin{bmatrix} -\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} +
|
|
\begin{bmatrix} \tilde{\partial}_x \\
|
|
\tilde{\partial}_y \end{bmatrix} T_a \epsilon_{zz}^{-1}
|
|
\begin{bmatrix} \hat{\partial}_x T_b \epsilon_{xx} & \hat{\partial}_y T_a \epsilon_{yy} \end{bmatrix})
|
|
\begin{bmatrix} E_r \\
|
|
E_y \end{bmatrix}
|
|
$$
|
|
|
|
for use with a field vector of the form `[E_r, E_y]`.
|
|
|
|
This operator can be used to form an eigenvalue problem of the form
|
|
A @ [E_r, E_y] = beta**2 * [E_r, E_y]
|
|
|
|
which can then be solved for the eigenmodes of the system
|
|
(an `exp(-i * angular_wavenumber * theta)` theta-dependence is assumed for
|
|
the fields, with `beta = angular_wavenumber / rmin`).
|
|
|
|
(NOTE: See module docs and 10.1364/OL.33.001848)
|
|
|
|
Args:
|
|
omega: The angular frequency of the system
|
|
dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
|
|
epsilon: Vectorized dielectric constant grid
|
|
rmin: Radius at the left edge of the simulation domain (at minimum 'x')
|
|
mu: Vectorized magnetic permeability grid (default 1 everywhere)
|
|
|
|
Returns:
|
|
Sparse matrix representation of the operator
|
|
"""
|
|
if mu is None:
|
|
mu = numpy.ones_like(epsilon)
|
|
|
|
Dfx, Dfy = deriv_forward(dxes[0])
|
|
Dbx, Dby = deriv_back(dxes[1])
|
|
|
|
Ta, Tb = dxes2T(dxes=dxes, rmin=rmin)
|
|
|
|
eps_parts = numpy.split(epsilon, 3)
|
|
eps_x = sparse.diags_array(eps_parts[0])
|
|
eps_y = sparse.diags_array(eps_parts[1])
|
|
eps_z_inv = sparse.diags_array(1 / eps_parts[2])
|
|
|
|
mu_parts = numpy.split(mu, 3)
|
|
mu_y = sparse.diags_array(mu_parts[1])
|
|
mu_x = sparse.diags_array(mu_parts[0])
|
|
mu_z_inv = sparse.diags_array(1 / mu_parts[2])
|
|
|
|
omega2 = omega * omega
|
|
diag = sparse.block_diag
|
|
|
|
sq0 = omega2 * diag((Tb @ Tb @ mu_y @ eps_x,
|
|
Ta @ Ta @ mu_x @ eps_y))
|
|
lin0 = sparse.vstack((-Tb @ mu_y @ Dby, Ta @ mu_x @ Dbx)) @ Tb @ mu_z_inv @ sparse.hstack((-Dfy, Dfx))
|
|
lin1 = sparse.vstack((Dfx, Dfy)) @ Ta @ eps_z_inv @ sparse.hstack((Dbx @ Tb @ eps_x,
|
|
Dby @ Ta @ eps_y))
|
|
op = sq0 + lin0 + lin1
|
|
return op
|
|
|
|
|
|
def solve_modes(
|
|
mode_numbers: Sequence[int],
|
|
omega: float,
|
|
dxes: dx_lists2_t,
|
|
epsilon: vfdslice,
|
|
rmin: float,
|
|
mu: vfdslice | None = None,
|
|
mode_margin: int = 2,
|
|
) -> tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]:
|
|
"""
|
|
Given a 2d (r, y) slice of epsilon, attempts to solve for the eigenmode
|
|
of the bent waveguide with the specified mode number.
|
|
|
|
Args:
|
|
mode_numbers: Mode numbers to solve, 0-indexed.
|
|
omega: Angular frequency of the simulation
|
|
dxes: Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types.
|
|
The first coordinate is assumed to be r, the second is y.
|
|
epsilon: Dielectric constant
|
|
rmin: Radius of curvature for the simulation. This should be the minimum value of
|
|
r within the simulation domain.
|
|
mu: Magnetic permeability (default 1 everywhere)
|
|
|
|
Returns:
|
|
e_xys: NDArray of vfdfield_t specifying fields. First dimension is mode number.
|
|
angular_wavenumbers: list of wavenumbers in 1/rad units.
|
|
"""
|
|
|
|
#
|
|
# Solve for the largest-magnitude eigenvalue of the real operator
|
|
#
|
|
dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes]
|
|
mu_real = None if mu is None else numpy.real(mu)
|
|
|
|
A_r = cylindrical_operator(numpy.real(omega), dxes_real, numpy.real(epsilon), rmin=rmin, mu=mu_real)
|
|
eigvals, eigvecs = signed_eigensolve(A_r, max(mode_numbers) + mode_margin)
|
|
keep_inds = -(numpy.array(mode_numbers) + 1)
|
|
e_xys = eigvecs[:, keep_inds].T
|
|
eigvals = eigvals[keep_inds]
|
|
|
|
#
|
|
# Now solve for the eigenvector of the full operator, using the real operator's
|
|
# eigenvector as an initial guess for Rayleigh quotient iteration.
|
|
#
|
|
A = cylindrical_operator(omega, dxes, epsilon, rmin=rmin, mu=mu)
|
|
for nn in range(len(mode_numbers)):
|
|
eigvals[nn], e_xys[nn, :] = rayleigh_quotient_iteration(A, e_xys[nn, :])
|
|
|
|
# Calculate the wave-vector (force the real part to be positive)
|
|
wavenumbers = numpy.sqrt(eigvals)
|
|
wavenumbers *= numpy.sign(numpy.real(wavenumbers))
|
|
|
|
# Wavenumbers assume the mode is at rmin, which is unlikely
|
|
# Instead, return the wavenumber in inverse radians
|
|
angular_wavenumbers = wavenumbers * cast('complex', rmin)
|
|
|
|
order = angular_wavenumbers.argsort()[::-1]
|
|
e_xys = e_xys[order]
|
|
angular_wavenumbers = angular_wavenumbers[order]
|
|
|
|
return e_xys, angular_wavenumbers
|
|
|
|
|
|
def solve_mode(
|
|
mode_number: int,
|
|
*args: Any,
|
|
**kwargs: Any,
|
|
) -> tuple[vcfdfield2, complex]:
|
|
"""
|
|
Wrapper around `solve_modes()` that solves for a single mode.
|
|
|
|
Args:
|
|
mode_number: 0-indexed mode number to solve for
|
|
*args: passed to `solve_modes()`
|
|
**kwargs: passed to `solve_modes()`
|
|
|
|
Returns:
|
|
(e_xy, angular_wavenumber)
|
|
"""
|
|
kwargs['mode_numbers'] = [mode_number]
|
|
e_xys, angular_wavenumbers = solve_modes(*args, **kwargs)
|
|
return e_xys[0], angular_wavenumbers[0]
|
|
|
|
|
|
def linear_wavenumbers(
|
|
e_xys: Sequence[vcfdfield2] | NDArray[numpy.complex128],
|
|
angular_wavenumbers: ArrayLike,
|
|
epsilon: vfdslice,
|
|
dxes: dx_lists2_t,
|
|
rmin: float,
|
|
) -> NDArray[numpy.complex128]:
|
|
"""
|
|
Calculate linear wavenumbers (1/distance) based on angular wavenumbers (1/rad)
|
|
and the mode's energy distribution.
|
|
|
|
Args:
|
|
e_xys: Vectorized mode fields with shape (num_modes, 2 * x *y)
|
|
angular_wavenumbers: Wavenumbers assuming fields have theta-dependence of
|
|
`exp(-i * angular_wavenumber * theta)`. They should satisfy
|
|
`operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy`
|
|
epsilon: Vectorized dielectric constant grid with shape (3, x, y)
|
|
dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
|
|
rmin: Radius at the left edge of the simulation domain (at minimum 'x')
|
|
|
|
Returns:
|
|
NDArray containing the calculated linear (1/distance) wavenumbers
|
|
"""
|
|
angular_wavenumbers = numpy.asarray(angular_wavenumbers)
|
|
mode_radii = numpy.empty_like(angular_wavenumbers, dtype=float)
|
|
|
|
shape2d = (len(dxes[0][0]), len(dxes[0][1]))
|
|
epsilon2d = unvec(epsilon, shape2d)[:2]
|
|
ra = rmin + numpy.cumsum(dxes[0][0])
|
|
rb = rmin + dxes[0][0] / 2.0 + numpy.concatenate((
|
|
numpy.zeros(1, dtype=dxes[1][0].dtype),
|
|
numpy.cumsum(dxes[1][0][:-1]),
|
|
))
|
|
for ii in range(angular_wavenumbers.size):
|
|
efield = unvec(e_xys[ii], shape2d, 2)
|
|
energy = numpy.real((efield * efield.conj()) * epsilon2d)
|
|
er_energy_vs_r = energy[0].sum(axis=1)
|
|
ey_energy_vs_r = energy[1].sum(axis=1)
|
|
energy_vs_r = er_energy_vs_r + ey_energy_vs_r
|
|
mode_radii[ii] = (
|
|
(rb * er_energy_vs_r).sum() + (ra * ey_energy_vs_r).sum()
|
|
) / energy_vs_r.sum()
|
|
|
|
logger.info(f'{mode_radii=}')
|
|
lin_wavenumbers = angular_wavenumbers / mode_radii
|
|
return lin_wavenumbers
|
|
|
|
|
|
def exy2h(
|
|
angular_wavenumber: complex,
|
|
omega: float,
|
|
dxes: dx_lists2_t,
|
|
rmin: float,
|
|
epsilon: vfdslice,
|
|
mu: vfdslice | None = None
|
|
) -> sparse.sparray:
|
|
"""
|
|
Operator which transforms the vector `e_xy` containing the vectorized E_r and E_y fields,
|
|
into a vectorized H containing all three H components
|
|
|
|
Args:
|
|
angular_wavenumber: Wavenumber assuming fields have theta-dependence of
|
|
`exp(-i * angular_wavenumber * theta)`. It should satisfy
|
|
`operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy`
|
|
omega: The angular frequency of the system
|
|
dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
|
|
rmin: Radius at the left edge of the simulation domain (at minimum 'x')
|
|
epsilon: Vectorized dielectric constant grid
|
|
mu: Vectorized magnetic permeability grid (default 1 everywhere)
|
|
|
|
Returns:
|
|
Sparse matrix representing the operator.
|
|
"""
|
|
e2hop = e2h(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, mu=mu)
|
|
return e2hop @ exy2e(angular_wavenumber=angular_wavenumber, dxes=dxes, rmin=rmin, epsilon=epsilon)
|
|
|
|
|
|
def exy2e(
|
|
angular_wavenumber: complex,
|
|
dxes: dx_lists2_t,
|
|
rmin: float,
|
|
epsilon: vfdslice,
|
|
) -> sparse.sparray:
|
|
"""
|
|
Operator which transforms the vector `e_xy` containing the vectorized E_r and E_y fields,
|
|
into a vectorized E containing all three E components
|
|
|
|
Unlike the straight waveguide case, the H_z components do not cancel and must be calculated
|
|
from E_r and E_y in order to then calculate E_z.
|
|
|
|
Args:
|
|
angular_wavenumber: Wavenumber assuming fields have theta-dependence of
|
|
`exp(-i * angular_wavenumber * theta)`. It should satisfy
|
|
`operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy`
|
|
dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
|
|
rmin: Radius at the left edge of the simulation domain (at minimum 'x')
|
|
epsilon: Vectorized dielectric constant grid
|
|
|
|
Returns:
|
|
Sparse matrix representing the operator.
|
|
"""
|
|
Dbx, Dby = deriv_back(dxes[1])
|
|
wavenumber = angular_wavenumber / rmin
|
|
|
|
Ta, Tb = dxes2T(dxes=dxes, rmin=rmin)
|
|
|
|
epsilon_parts = numpy.split(epsilon, 3)
|
|
epsilon_x, epsilon_y = (sparse.diags_array(epsi) for epsi in epsilon_parts[:2])
|
|
epsilon_z_inv = sparse.diags_array(1 / epsilon_parts[2])
|
|
|
|
n_pts = dxes[0][0].size * dxes[0][1].size
|
|
exy2ez = (
|
|
Ta @ epsilon_z_inv
|
|
@ sparse.hstack((Dbx @ Tb @ epsilon_x,
|
|
Dby @ Ta @ epsilon_y))
|
|
/ (1j * wavenumber)
|
|
)
|
|
|
|
op = sparse.vstack((sparse.eye_array(2 * n_pts),
|
|
exy2ez))
|
|
return op
|
|
|
|
|
|
def e2h(
|
|
angular_wavenumber: complex,
|
|
omega: float,
|
|
dxes: dx_lists2_t,
|
|
rmin: float,
|
|
mu: vfdslice | None = None
|
|
) -> sparse.sparray:
|
|
r"""
|
|
Returns an operator which, when applied to a vectorized E eigenfield, produces
|
|
the vectorized H eigenfield.
|
|
|
|
This operator is created directly from the initial coordinate-transformed equations:
|
|
$$
|
|
\begin{aligned}
|
|
-\imath \omega \mu_{rr} H_r &= \tilde{\partial}_y E_z + \imath \beta T_a^{-1} E_y, \\
|
|
-\imath \omega \mu_{yy} H_y &= -\imath \beta T_b^{-1} E_r
|
|
- T_b^{-1} \tilde{\partial}_r (T_a E_z), \\
|
|
-\imath \omega \mu_{zz} H_z &= \tilde{\partial}_r E_y - \tilde{\partial}_y E_r,
|
|
\end{aligned}
|
|
$$
|
|
|
|
Args:
|
|
angular_wavenumber: Wavenumber assuming fields have theta-dependence of
|
|
`exp(-i * angular_wavenumber * theta)`. It should satisfy
|
|
`operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy`
|
|
omega: The angular frequency of the system
|
|
dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
|
|
rmin: Radius at the left edge of the simulation domain (at minimum 'x')
|
|
mu: Vectorized magnetic permeability grid (default 1 everywhere)
|
|
|
|
Returns:
|
|
Sparse matrix representation of the operator.
|
|
"""
|
|
Dfx, Dfy = deriv_forward(dxes[0])
|
|
Ta, Tb = dxes2T(dxes=dxes, rmin=rmin)
|
|
Tai = sparse.diags_array(1 / Ta.diagonal())
|
|
Tbi = sparse.diags_array(1 / Tb.diagonal())
|
|
|
|
jB = 1j * angular_wavenumber / rmin
|
|
op = sparse.block_array([[ None, jB * Tai, Dfy],
|
|
[-jB * Tbi, None, -Tbi @ Dfx @ Ta],
|
|
[ -Dfy, Dfx, None]]) / (-1j * omega)
|
|
if mu is not None:
|
|
op = sparse.diags_array(1 / mu) @ op
|
|
return op
|
|
|
|
|
|
def dxes2T(
|
|
dxes: dx_lists2_t,
|
|
rmin: float,
|
|
) -> tuple[NDArray[numpy.float64], NDArray[numpy.float64]]:
|
|
r"""
|
|
Construct the cylindrical metric matrices $T_a$ and $T_b$.
|
|
|
|
`T_a` is sampled on the E-grid radial locations, while `T_b` is sampled on
|
|
the staggered H-grid radial locations. These are the diagonal matrices that
|
|
convert the straight-waveguide algebra into its cylindrical counterpart.
|
|
|
|
Args:
|
|
dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
|
|
rmin: Radius at the left edge of the simulation domain (at minimum 'x')
|
|
|
|
Returns:
|
|
Sparse diagonal matrices `(T_a, T_b)`.
|
|
"""
|
|
ra = rmin + numpy.cumsum(dxes[0][0]) # Radius at Ey points
|
|
rb = (
|
|
rmin
|
|
+ dxes[0][0] / 2.0
|
|
+ numpy.concatenate((
|
|
numpy.zeros(1, dtype=dxes[1][0].dtype),
|
|
numpy.cumsum(dxes[1][0][:-1]),
|
|
))
|
|
) # Radius at Er points
|
|
ta = ra / rmin
|
|
tb = rb / rmin
|
|
|
|
Ta = sparse.diags_array(vec(ta[:, None].repeat(dxes[0][1].size, axis=1)))
|
|
Tb = sparse.diags_array(vec(tb[:, None].repeat(dxes[1][1].size, axis=1)))
|
|
return Ta, Tb
|
|
|
|
|
|
def normalized_fields_e(
|
|
e_xy: vcfdfield2,
|
|
angular_wavenumber: complex,
|
|
omega: float,
|
|
dxes: dx_lists2_t,
|
|
rmin: float,
|
|
epsilon: vfdslice,
|
|
mu: vfdslice | None = None,
|
|
prop_phase: float = 0,
|
|
) -> tuple[vcfdslice_t, vcfdslice_t]:
|
|
r"""
|
|
Given a vector `e_xy` containing the vectorized E_r and E_y fields,
|
|
returns normalized, vectorized E and H fields for the system.
|
|
|
|
Args:
|
|
e_xy: Vector containing E_r and E_y fields
|
|
angular_wavenumber: Wavenumber assuming fields have theta-dependence of
|
|
`exp(-i * angular_wavenumber * theta)`. It should satisfy
|
|
`operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy`
|
|
omega: The angular frequency of the system
|
|
dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
|
|
rmin: Radius at the left edge of the simulation domain (at minimum 'x')
|
|
epsilon: Vectorized dielectric constant grid
|
|
mu: Vectorized magnetic permeability grid (default 1 everywhere)
|
|
prop_phase: Phase shift `(dz * corrected_wavenumber)` over 1 cell in propagation direction.
|
|
Default 0 (continuous propagation direction, i.e. dz->0).
|
|
|
|
Returns:
|
|
`(e, h)`, where each field is vectorized, normalized,
|
|
and contains all three vector components.
|
|
|
|
Notes:
|
|
The normalization step is delegated to `_normalized_fields(...)`, which
|
|
enforces unit forward power under the discrete inner product
|
|
|
|
$$
|
|
\frac{1}{2}\int (E_r H_y^* - E_y H_r^*) \, dr \, dy.
|
|
$$
|
|
|
|
The angular wavenumber `m` is first converted into the full three-component
|
|
fields, then the overall complex phase and sign are fixed so the result is
|
|
reproducible for symmetric modes.
|
|
"""
|
|
e = exy2e(angular_wavenumber=angular_wavenumber, dxes=dxes, rmin=rmin, epsilon=epsilon) @ e_xy
|
|
h = exy2h(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon, mu=mu) @ e_xy
|
|
e_norm, h_norm = _normalized_fields(
|
|
e=e, h=h, dxes=dxes, epsilon=epsilon, prop_phase=prop_phase,
|
|
)
|
|
return e_norm, h_norm
|
|
|
|
|
|
def _normalized_fields(
|
|
e: vcfdslice,
|
|
h: vcfdslice,
|
|
dxes: dx_lists2_t,
|
|
epsilon: vfdslice,
|
|
prop_phase: float = 0,
|
|
) -> tuple[vcfdslice_t, vcfdslice_t]:
|
|
r"""
|
|
Normalize a cylindrical waveguide mode to unit forward power.
|
|
|
|
The cylindrical helpers reuse the straight-waveguide inner product after the
|
|
field reconstruction step. The extra metric factors have already been folded
|
|
into the reconstructed `e`/`h` fields through `dxes2T(...)` and the
|
|
cylindrical `exy2e(...)` / `exy2h(...)` operators, so the same discrete
|
|
longitudinal Poynting integral can be used here.
|
|
|
|
The normalization procedure is:
|
|
|
|
1. Compute the time-averaged forward power with
|
|
`waveguide_2d.inner_product(..., conj_h=True)`.
|
|
2. Scale by `1 / sqrt(S_z)` so the resulting mode has unit forward power.
|
|
3. Remove the arbitrary complex phase and apply a quadrant-sum sign heuristic
|
|
so symmetric modes choose a stable representative.
|
|
|
|
`prop_phase` has the same meaning as in `waveguide_2d`: it compensates for
|
|
the half-cell longitudinal staggering between the E and H fields when the
|
|
propagation direction is itself discretized.
|
|
"""
|
|
shape = [s.size for s in dxes[0]]
|
|
|
|
# Find time-averaged Sz and normalize to it
|
|
# H phase is adjusted by a half-cell forward shift for Yee cell, and 1-cell reverse shift for Poynting
|
|
Sz_tavg = waveguide_2d.inner_product(e, h, dxes=dxes, prop_phase=prop_phase, conj_h=True).real # Note, using linear poynting vector
|
|
assert Sz_tavg > 0, f'Found a mode propagating in the wrong direction! {Sz_tavg=}'
|
|
|
|
energy = numpy.real(epsilon * e.conj() * e)
|
|
|
|
norm_amplitude = 1 / numpy.sqrt(Sz_tavg)
|
|
norm_angle = -numpy.angle(e[energy.argmax()]) # Will randomly add a negative sign when mode is symmetric
|
|
|
|
# Try to break symmetry to assign a consistent sign [experimental]
|
|
E_weighted = unvec(e * energy * numpy.exp(1j * norm_angle), shape)
|
|
sign = numpy.sign(E_weighted[:,
|
|
:max(shape[0] // 2, 1),
|
|
:max(shape[1] // 2, 1)].real.sum())
|
|
assert sign != 0
|
|
|
|
norm_factor = sign * norm_amplitude * numpy.exp(1j * norm_angle)
|
|
|
|
e *= norm_factor
|
|
h *= norm_factor
|
|
return vcfdslice_t(e), vcfdslice_t(h)
|