waveguides

meanas.fdfd.waveguide_2d

Operators and helper functions for waveguides with unchanging cross-section.

The propagation direction is chosen to be along the z axis, and all fields are given an implicit z-dependence of the form exp(-1 * wavenumber * z).

As the z-dependence is known, all the functions in this file assume a 2D grid (i.e. dxes = [[[dx_e[0], dx_e[1], ...], [dy_e[0], ...]], [[dx_h[0], ...], [dy_h[0], ...]]]).

===============

Consider Maxwell's equations in continuous space, in the frequency domain. Assuming a structure with some (x, y) cross-section extending uniformly into the z dimension, with a diagonal \(\epsilon\) tensor, we have

\[ \begin{aligned} \nabla \times \vec{E}(x, y, z) &= -\imath \omega \mu \vec{H} \\ \nabla \times \vec{H}(x, y, z) &= \imath \omega \epsilon \vec{E} \\ \vec{E}(x,y,z) &= (\vec{E}_t(x, y) + E_z(x, y)\vec{z}) e^{-\imath \beta z} \\ \vec{H}(x,y,z) &= (\vec{H}_t(x, y) + H_z(x, y)\vec{z}) e^{-\imath \beta z} \\ \end{aligned} \]

Expanding the first two equations into vector components, we get

\[ \begin{aligned} -\imath \omega \mu_{xx} H_x &= \partial_y E_z - \partial_z E_y \\ -\imath \omega \mu_{yy} H_y &= \partial_z E_x - \partial_x E_z \\ -\imath \omega \mu_{zz} H_z &= \partial_x E_y - \partial_y E_x \\ \imath \omega \epsilon_{xx} E_x &= \partial_y H_z - \partial_z H_y \\ \imath \omega \epsilon_{yy} E_y &= \partial_z H_x - \partial_x H_z \\ \imath \omega \epsilon_{zz} E_z &= \partial_x H_y - \partial_y H_x \\ \end{aligned} \]

Substituting in our expressions for \(\vec{E}\), \(\vec{H}\) and discretizing:

\[ \begin{aligned} -\imath \omega \mu_{xx} H_x &= \tilde{\partial}_y E_z + \imath \beta E_y \\ -\imath \omega \mu_{yy} H_y &= -\imath \beta E_x - \tilde{\partial}_x E_z \\ -\imath \omega \mu_{zz} H_z &= \tilde{\partial}_x E_y - \tilde{\partial}_y E_x \\ \imath \omega \epsilon_{xx} E_x &= \hat{\partial}_y H_z + \imath \beta H_y \\ \imath \omega \epsilon_{yy} E_y &= -\imath \beta H_x - \hat{\partial}_x H_z \\ \imath \omega \epsilon_{zz} E_z &= \hat{\partial}_x H_y - \hat{\partial}_y H_x \\ \end{aligned} \]

Rewrite the last three equations as

\[ \begin{aligned} \imath \beta H_y &= \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z \\ \imath \beta H_x &= -\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z \\ \imath \omega E_z &= \frac{1}{\epsilon_{zz}} \hat{\partial}_x H_y - \frac{1}{\epsilon_{zz}} \hat{\partial}_y H_x \\ \end{aligned} \]

Now apply \(\imath \beta \tilde{\partial}_x\) to the last equation, then substitute in for \(\imath \beta H_x\) and \(\imath \beta H_y\):

\[ \begin{aligned} \imath \beta \tilde{\partial}_x \imath \omega E_z &= \imath \beta \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x H_y - \imath \beta \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y H_x \\ &= \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x ( \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z) - \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y (-\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z) \\ &= \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x ( \imath \omega \epsilon_{xx} E_x) - \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y (-\imath \omega \epsilon_{yy} E_y) \\ \imath \beta \tilde{\partial}_x E_z &= \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) \\ \end{aligned} \]

With a similar approach (but using \(\imath \beta \tilde{\partial}_y\) instead), we can get

\[ \begin{aligned} \imath \beta \tilde{\partial}_y E_z &= \tilde{\partial}_y \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + \tilde{\partial}_y \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) \\ \end{aligned} \]

We can combine this equation for \(\imath \beta \tilde{\partial}_y E_z\) with the unused \(\imath \omega \mu_{xx} H_x\) and \(\imath \omega \mu_{yy} H_y\) equations to get

\[ \begin{aligned} -\imath \omega \mu_{xx} \imath \beta H_x &= -\beta^2 E_y + \imath \beta \tilde{\partial}_y E_z \\ -\imath \omega \mu_{xx} \imath \beta H_x &= -\beta^2 E_y + \tilde{\partial}_y ( \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) )\\ \end{aligned} \]

and

\[ \begin{aligned} -\imath \omega \mu_{yy} \imath \beta H_y &= \beta^2 E_x - \imath \beta \tilde{\partial}_x E_z \\ -\imath \omega \mu_{yy} \imath \beta H_y &= \beta^2 E_x - \tilde{\partial}_x ( \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) )\\ \end{aligned} \]

However, based on our rewritten equation for \(\imath \beta H_x\) and the so-far unused equation for \(\imath \omega \mu_{zz} H_z\) we can also write

\[ \begin{aligned} -\imath \omega \mu_{xx} (\imath \beta H_x) &= -\imath \omega \mu_{xx} (-\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z) \\ &= -\omega^2 \mu_{xx} \epsilon_{yy} E_y + \imath \omega \mu_{xx} \hat{\partial}_x ( \frac{1}{-\imath \omega \mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x)) \\ &= -\omega^2 \mu_{xx} \epsilon_{yy} E_y -\mu_{xx} \hat{\partial}_x \frac{1}{\mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x) \\ \end{aligned} \]

and, similarly,

\[ \begin{aligned} -\imath \omega \mu_{yy} (\imath \beta H_y) &= \omega^2 \mu_{yy} \epsilon_{xx} E_x +\mu_{yy} \hat{\partial}_y \frac{1}{\mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x) \\ \end{aligned} \]

By combining both pairs of expressions, we get

\[ \begin{aligned} \beta^2 E_x - \tilde{\partial}_x ( \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) ) &= \omega^2 \mu_{yy} \epsilon_{xx} E_x +\mu_{yy} \hat{\partial}_y \frac{1}{\mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x) \\ -\beta^2 E_y + \tilde{\partial}_y ( \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x) + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) ) &= -\omega^2 \mu_{xx} \epsilon_{yy} E_y -\mu_{xx} \hat{\partial}_x \frac{1}{\mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x) \\ \end{aligned} \]

Using these, we can construct the eigenvalue problem

\[ \beta^2 \begin{bmatrix} E_x \\ E_y \end{bmatrix} = (\omega^2 \begin{bmatrix} \mu_{yy} \epsilon_{xx} & 0 \\ 0 & \mu_{xx} \epsilon_{yy} \end{bmatrix} + \begin{bmatrix} -\mu_{yy} \hat{\partial}_y \\ \mu_{xx} \hat{\partial}_x \end{bmatrix} \mu_{zz}^{-1} \begin{bmatrix} -\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} + \begin{bmatrix} \tilde{\partial}_x \\ \tilde{\partial}_y \end{bmatrix} \epsilon_{zz}^{-1} \begin{bmatrix} \hat{\partial}_x \epsilon_{xx} & \hat{\partial}_y \epsilon_{yy} \end{bmatrix}) \begin{bmatrix} E_x \\ E_y \end{bmatrix} \]

In the literature, \(\beta\) is usually used to denote the lossless/real part of the propagation constant, but in meanas it is allowed to be complex.

An equivalent eigenvalue problem can be formed using the \(H_x\) and \(H_y\) fields, if those are more convenient.

Note that \(E_z\) was never discretized, so \(\beta\) will need adjustment to account for numerical dispersion if the result is introduced into a space with a discretized z-axis.

operator_e

operator_e(
    omega: complex,
    dxes: dx_lists2_t,
    epsilon: vfdslice,
    mu: vfdslice | None = None,
) -> sparse.sparray

Waveguide operator of the form

omega**2 * mu * epsilon +
mu * [[-Dy], [Dx]] / mu * [-Dy, Dx] +
[[Dx], [Dy]] / epsilon * [Dx, Dy] * epsilon

for use with a field vector of the form cat([E_x, E_y]).

More precisely, the operator is

\[ \omega^2 \begin{bmatrix} \mu_{yy} \epsilon_{xx} & 0 \\ 0 & \mu_{xx} \epsilon_{yy} \end{bmatrix} + \begin{bmatrix} -\mu_{yy} \hat{\partial}_y \\ \mu_{xx} \hat{\partial}_x \end{bmatrix} \mu_{zz}^{-1} \begin{bmatrix} -\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} + \begin{bmatrix} \tilde{\partial}_x \\ \tilde{\partial}_y \end{bmatrix} \epsilon_{zz}^{-1} \begin{bmatrix} \hat{\partial}_x \epsilon_{xx} & \hat{\partial}_y \epsilon_{yy} \end{bmatrix} \]

\(\tilde{\partial}_x\) and \(\hat{\partial}_x\) are the forward and backward derivatives along x, and each \(\epsilon_{xx}\), \(\mu_{yy}\), etc. is a diagonal matrix containing the vectorized material property distribution.

This operator can be used to form an eigenvalue problem of the form operator_e(...) @ [E_x, E_y] = wavenumber**2 * [E_x, E_y]

which can then be solved for the eigenmodes of the system (an exp(-i * wavenumber * z) z-dependence is assumed for the fields).

Parameters:

Name	Type	Description	Default
omega	complex	The angular frequency of the system.	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
epsilon	vfdslice	Vectorized dielectric constant grid	required
mu	vfdslice | None	Vectorized magnetic permeability grid (default 1 everywhere)	None

Returns:

Type	Description
sparray	Sparse matrix representation of the operator.

operator_h

operator_h(
    omega: complex,
    dxes: dx_lists2_t,
    epsilon: vfdslice,
    mu: vfdslice | None = None,
) -> sparse.sparray

Waveguide operator of the form

omega**2 * epsilon * mu +
epsilon * [[-Dy], [Dx]] / epsilon * [-Dy, Dx] +
[[Dx], [Dy]] / mu * [Dx, Dy] * mu

for use with a field vector of the form cat([H_x, H_y]).

More precisely, the operator is

\[ \omega^2 \begin{bmatrix} \epsilon_{yy} \mu_{xx} & 0 \\ 0 & \epsilon_{xx} \mu_{yy} \end{bmatrix} + \begin{bmatrix} -\epsilon_{yy} \tilde{\partial}_y \\ \epsilon_{xx} \tilde{\partial}_x \end{bmatrix} \epsilon_{zz}^{-1} \begin{bmatrix} -\hat{\partial}_y & \hat{\partial}_x \end{bmatrix} + \begin{bmatrix} \hat{\partial}_x \\ \hat{\partial}_y \end{bmatrix} \mu_{zz}^{-1} \begin{bmatrix} \tilde{\partial}_x \mu_{xx} & \tilde{\partial}_y \mu_{yy} \end{bmatrix} \]

\(\tilde{\partial}_x\) and \(\hat{\partial}_x\) are the forward and backward derivatives along x, and each \(\epsilon_{xx}\), \(\mu_{yy}\), etc. is a diagonal matrix containing the vectorized material property distribution.

This operator can be used to form an eigenvalue problem of the form operator_h(...) @ [H_x, H_y] = wavenumber**2 * [H_x, H_y]

which can then be solved for the eigenmodes of the system (an exp(-i * wavenumber * z) z-dependence is assumed for the fields).

Parameters:

Name	Type	Description	Default
omega	complex	The angular frequency of the system.	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
epsilon	vfdslice	Vectorized dielectric constant grid	required
mu	vfdslice | None	Vectorized magnetic permeability grid (default 1 everywhere)	None

Returns:

Type	Description
sparray	Sparse matrix representation of the operator.

normalized_fields_e

normalized_fields_e(
    e_xy: vcfdfield2,
    wavenumber: complex,
    omega: complex,
    dxes: dx_lists2_t,
    epsilon: vfdslice,
    mu: vfdslice | None = None,
    prop_phase: float = 0,
) -> tuple[vcfdslice_t, vcfdslice_t]

Given a vector e_xy containing the vectorized E_x and E_y fields, returns normalized, vectorized E and H fields for the system.

Parameters:

Name	Type	Description	Default
e_xy	vcfdfield2	Vector containing E_x and E_y fields	required
wavenumber	complex	Wavenumber assuming fields have z-dependence of exp(-i * wavenumber * z). It should satisfy operator_e() @ e_xy == wavenumber**2 * e_xy	required
omega	complex	The angular frequency of the system	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
epsilon	vfdslice	Vectorized dielectric constant grid	required
mu	vfdslice | None	Vectorized magnetic permeability grid (default 1 everywhere)	None
prop_phase	float	Phase shift (dz * corrected_wavenumber) over 1 cell in propagation direction. Default 0 (continuous propagation direction, i.e. dz->0).	0

Returns:

Type	Description
vcfdslice_t	(e, h), where each field is vectorized, normalized,
vcfdslice_t	and contains all three vector components.

Notes

e_xy is only the transverse electric eigenvector. This helper first reconstructs the full three-component E and H fields with exy2e(...) and exy2h(...), then normalizes them to unit forward power using _normalized_fields(...).

The normalization target is

\[ \Re\left[\mathrm{inner\_product}(e, h, \mathrm{conj\_h}=True)\right] = 1, \]

so the returned fields represent a unit-power forward mode under the discrete Yee-grid Poynting inner product.

normalized_fields_h

normalized_fields_h(
    h_xy: vcfdfield2,
    wavenumber: complex,
    omega: complex,
    dxes: dx_lists2_t,
    epsilon: vfdslice,
    mu: vfdslice | None = None,
    prop_phase: float = 0,
) -> tuple[vcfdslice_t, vcfdslice_t]

Given a vector h_xy containing the vectorized H_x and H_y fields, returns normalized, vectorized E and H fields for the system.

Parameters:

Name	Type	Description	Default
h_xy	vcfdfield2	Vector containing H_x and H_y fields	required
wavenumber	complex	Wavenumber assuming fields have z-dependence of exp(-i * wavenumber * z). It should satisfy operator_h() @ h_xy == wavenumber**2 * h_xy	required
omega	complex	The angular frequency of the system	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
epsilon	vfdslice	Vectorized dielectric constant grid	required
mu	vfdslice | None	Vectorized magnetic permeability grid (default 1 everywhere)	None
prop_phase	float	Phase shift (dz * corrected_wavenumber) over 1 cell in propagation direction. Default 0 (continuous propagation direction, i.e. dz->0).	0

Returns:

Type	Description
vcfdslice_t	(e, h), where each field is vectorized, normalized,
vcfdslice_t	and contains all three vector components.

Notes

This is the H_x/H_y analogue of normalized_fields_e(...). The final normalized mode should describe the same physical solution, but because the overall complex phase and sign are chosen heuristically, normalized_fields_e(...) and normalized_fields_h(...) need not return identical representatives for nearly symmetric modes.

exy2h

exy2h(
    wavenumber: complex,
    omega: complex,
    dxes: dx_lists2_t,
    epsilon: vfdslice,
    mu: vfdslice | None = None,
) -> sparse.sparray

Operator which transforms the vector e_xy containing the vectorized E_x and E_y fields, into a vectorized H containing all three H components

Parameters:

Name	Type	Description	Default
wavenumber	complex	Wavenumber assuming fields have z-dependence of exp(-i * wavenumber * z). It should satisfy operator_e() @ e_xy == wavenumber**2 * e_xy	required
omega	complex	The angular frequency of the system	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
epsilon	vfdslice	Vectorized dielectric constant grid	required
mu	vfdslice | None	Vectorized magnetic permeability grid (default 1 everywhere)	None

Returns:

Type	Description
sparray	Sparse matrix representing the operator.

hxy2e

hxy2e(
    wavenumber: complex,
    omega: complex,
    dxes: dx_lists2_t,
    epsilon: vfdslice,
    mu: vfdslice | None = None,
) -> sparse.sparray

Operator which transforms the vector h_xy containing the vectorized H_x and H_y fields, into a vectorized E containing all three E components

Parameters:

Name	Type	Description	Default
wavenumber	complex	Wavenumber assuming fields have z-dependence of exp(-i * wavenumber * z). It should satisfy operator_h() @ h_xy == wavenumber**2 * h_xy	required
omega	complex	The angular frequency of the system	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
epsilon	vfdslice	Vectorized dielectric constant grid	required
mu	vfdslice | None	Vectorized magnetic permeability grid (default 1 everywhere)	None

Returns:

Type	Description
sparray	Sparse matrix representing the operator.

hxy2h

hxy2h(
    wavenumber: complex,
    dxes: dx_lists2_t,
    mu: vfdslice | None = None,
) -> sparse.sparray

Operator which transforms the vector h_xy containing the vectorized H_x and H_y fields, into a vectorized H containing all three H components

Parameters:

Name	Type	Description	Default
wavenumber	complex	Wavenumber assuming fields have z-dependence of exp(-i * wavenumber * z). It should satisfy operator_h() @ h_xy == wavenumber**2 * h_xy	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
mu	vfdslice | None	Vectorized magnetic permeability grid (default 1 everywhere)	None

Returns:

Type	Description
sparray	Sparse matrix representing the operator.

exy2e

exy2e(
    wavenumber: complex,
    dxes: dx_lists2_t,
    epsilon: vfdslice,
) -> sparse.sparray

Operator which transforms the vector e_xy containing the vectorized E_x and E_y fields, into a vectorized E containing all three E components

From the operator derivation (see module docs), we have

\[ \imath \omega \epsilon_{zz} E_z = \hat{\partial}_x H_y - \hat{\partial}_y H_x \\ \]

as well as the intermediate equations

\[ \begin{aligned} \imath \beta H_y &= \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z \\ \imath \beta H_x &= -\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z \\ \end{aligned} \]

Combining these, we get

\[ \begin{aligned} E_z &= \frac{1}{- \omega \beta \epsilon_{zz}} (( \hat{\partial}_y \hat{\partial}_x H_z -\hat{\partial}_x \hat{\partial}_y H_z) + \imath \omega (\hat{\partial}_x \epsilon_{xx} E_x + \hat{\partial}_y \epsilon{yy} E_y)) &= \frac{1}{\imath \beta \epsilon_{zz}} (\hat{\partial}_x \epsilon_{xx} E_x + \hat{\partial}_y \epsilon{yy} E_y) \end{aligned} \]

Parameters:

Name	Type	Description	Default
wavenumber	complex	Wavenumber assuming fields have z-dependence of exp(-i * wavenumber * z) It should satisfy operator_e() @ e_xy == wavenumber**2 * e_xy	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
epsilon	vfdslice	Vectorized dielectric constant grid	required

Returns:

Type	Description
sparray	Sparse matrix representing the operator.

e2h

e2h(
    wavenumber: complex,
    omega: complex,
    dxes: dx_lists2_t,
    mu: vfdslice | None = None,
) -> sparse.sparray

Returns an operator which, when applied to a vectorized E eigenfield, produces the vectorized H eigenfield slice.

Parameters:

Name	Type	Description	Default
wavenumber	complex	Wavenumber assuming fields have z-dependence of exp(-i * wavenumber * z)	required
omega	complex	The angular frequency of the system	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
mu	vfdslice | None	Vectorized magnetic permeability grid (default 1 everywhere)	None

Returns:

Type	Description
sparray	Sparse matrix representation of the operator.

h2e

h2e(
    wavenumber: complex,
    omega: complex,
    dxes: dx_lists2_t,
    epsilon: vfdslice,
) -> sparse.sparray

Returns an operator which, when applied to a vectorized H eigenfield, produces the vectorized E eigenfield slice.

Parameters:

Name	Type	Description	Default
wavenumber	complex	Wavenumber assuming fields have z-dependence of exp(-i * wavenumber * z)	required
omega	complex	The angular frequency of the system	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
epsilon	vfdslice	Vectorized dielectric constant grid	required

Returns:

Type	Description
sparray	Sparse matrix representation of the operator.

curl_e

curl_e(
    wavenumber: complex, dxes: dx_lists2_t
) -> sparse.sparray

Discretized curl operator for use with the waveguide E field slice.

Parameters:

Name	Type	Description	Default
wavenumber	complex	Wavenumber assuming fields have z-dependence of exp(-i * wavenumber * z)	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required

Returns:

Type	Description
sparray	Sparse matrix representation of the operator.

curl_h

curl_h(
    wavenumber: complex, dxes: dx_lists2_t
) -> sparse.sparray

Discretized curl operator for use with the waveguide H field slice.

Parameters:

Name	Type	Description	Default
wavenumber	complex	Wavenumber assuming fields have z-dependence of exp(-i * wavenumber * z)	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required

Returns:

Type	Description
sparray	Sparse matrix representation of the operator.

h_err

h_err(
    h: vcfdslice,
    wavenumber: complex,
    omega: complex,
    dxes: dx_lists2_t,
    epsilon: vfdslice,
    mu: vfdslice | None = None,
) -> float

Calculates the relative error in the H field

Parameters:

Name	Type	Description	Default
h	vcfdslice	Vectorized H field	required
wavenumber	complex	Wavenumber assuming fields have z-dependence of exp(-i * wavenumber * z)	required
omega	complex	The angular frequency of the system	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
epsilon	vfdslice	Vectorized dielectric constant grid	required
mu	vfdslice | None	Vectorized magnetic permeability grid (default 1 everywhere)	None

Returns:

Type	Description
float	Relative error norm(A_h @ h) / norm(h).

e_err

e_err(
    e: vcfdslice,
    wavenumber: complex,
    omega: complex,
    dxes: dx_lists2_t,
    epsilon: vfdslice,
    mu: vfdslice | None = None,
) -> float

Calculates the relative error in the E field

Parameters:

Name	Type	Description	Default
e	vcfdslice	Vectorized E field	required
wavenumber	complex	Wavenumber assuming fields have z-dependence of exp(-i * wavenumber * z)	required
omega	complex	The angular frequency of the system	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
epsilon	vfdslice	Vectorized dielectric constant grid	required
mu	vfdslice | None	Vectorized magnetic permeability grid (default 1 everywhere)	None

Returns:

Type	Description
float	Relative error norm(A_e @ e) / norm(e).

sensitivity

sensitivity(
    e_norm: vcfdslice,
    h_norm: vcfdslice,
    wavenumber: complex,
    omega: complex,
    dxes: dx_lists2_t,
    epsilon: vfdslice,
    mu: vfdslice | None = None,
) -> vcfdslice_t

Given a waveguide structure (dxes, epsilon, mu) and mode fields (e_norm, h_norm, wavenumber, omega), calculates the sensitivity of the wavenumber \(\beta\) to changes in the dielectric structure \(\epsilon\).

The output is a vector of the same size as vec(epsilon), with each element specifying the sensitivity of wavenumber to changes in the corresponding element in vec(epsilon), i.e.

\[ sens_{i} = \frac{\partial\beta}{\partial\epsilon_i} \]

An adjoint approach is used to calculate the sensitivity; the derivation is provided here:

Starting with the eigenvalue equation

\[ \beta^2 E_{xy} = A_E E_{xy} \]

where \(A_E\) is the waveguide operator from operator_e(), and \(E_{xy} = \begin{bmatrix} E_x \\ E_y \end{bmatrix}\), we can differentiate with respect to one of the \(\epsilon\) elements (i.e. at one Yee grid point), \(\epsilon_i\):

\[ (2 \beta) \partial_{\epsilon_i}(\beta) E_{xy} + \beta^2 \partial_{\epsilon_i} E_{xy} = \partial_{\epsilon_i}(A_E) E_{xy} + A_E \partial_{\epsilon_i} E_{xy} \]

We then multiply by \(H_{yx}^\star = \begin{bmatrix}H_y^\star \\ -H_x^\star \end{bmatrix}\) from the left:

\[ (2 \beta) \partial_{\epsilon_i}(\beta) H_{yx}^\star E_{xy} + \beta^2 H_{yx}^\star \partial_{\epsilon_i} E_{xy} = H_{yx}^\star \partial_{\epsilon_i}(A_E) E_{xy} + H_{yx}^\star A_E \partial_{\epsilon_i} E_{xy} \]

However, \(H_{yx}^\star\) is actually a left-eigenvector of \(A_E\). This can be verified by inspecting the form of operator_h (\(A_H\)) and comparing its conjugate transpose to operator_e (\(A_E\)). Also, note \(H_{yx}^\star \cdot E_{xy} = H^\star \times E\) recalls the mode orthogonality relation. See doi:10.5194/ars-9-85-201 for a similar approach. Therefore,

\[ H_{yx}^\star A_E \partial_{\epsilon_i} E_{xy} = \beta^2 H_{yx}^\star \partial_{\epsilon_i} E_{xy} \]

and we can simplify to

\[ \partial_{\epsilon_i}(\beta) = \frac{1}{2 \beta} \frac{H_{yx}^\star \partial_{\epsilon_i}(A_E) E_{xy} }{H_{yx}^\star E_{xy}} \]

This expression can be quickly calculated for all \(i\) by writing out the various terms of \(\partial_{\epsilon_i} A_E\) and recognizing that the vector-matrix-vector products (i.e. scalars) \(sens_i = \vec{v}_{left} \partial_{\epsilon_i} (\epsilon_{xyz}) \vec{v}_{right}\), indexed by \(i\), can be expressed as elementwise multiplications \(\vec{sens} = \vec{v}_{left} \star \vec{v}_{right}\)

Parameters:

Name	Type	Description	Default
e_norm	vcfdslice	Normalized, vectorized E_xyz field for the mode. E.g. as returned by normalized_fields_e.	required
h_norm	vcfdslice	Normalized, vectorized H_xyz field for the mode. E.g. as returned by normalized_fields_e.	required
wavenumber	complex	Propagation constant for the mode. The z-axis is assumed to be continuous (i.e. without numerical dispersion).	required
omega	complex	The angular frequency of the system.	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
epsilon	vfdslice	Vectorized dielectric constant grid	required
mu	vfdslice | None	Vectorized magnetic permeability grid (default 1 everywhere)	None

Returns:

Type	Description
vcfdslice_t	Sparse matrix representation of the operator.

solve_modes

solve_modes(
    mode_numbers: Sequence[int],
    omega: complex,
    dxes: dx_lists2_t,
    epsilon: vfdslice,
    mu: vfdslice | None = None,
    mode_margin: int = 2,
) -> tuple[
    NDArray[numpy.complex128], NDArray[numpy.complex128]
]

Given a 2D region, attempts to solve for the eigenmode with the specified mode numbers.

Parameters:

Name	Type	Description	Default
mode_numbers	Sequence[int]	List of 0-indexed mode numbers to solve for	required
omega	complex	Angular frequency of the simulation	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types	required
epsilon	vfdslice	Dielectric constant	required
mu	vfdslice | None	Magnetic permeability (default 1 everywhere)	None
mode_margin	int	The eigensolver will actually solve for (max(mode_number) + mode_margin) modes, but only return the target mode. Increasing this value can improve the solver's ability to find the correct mode. Default 2.	2

Returns:

Name	Type	Description
e_xys	NDArray[complex128]	NDArray of vfdfield_t specifying fields. First dimension is mode number.
wavenumbers	NDArray[complex128]	list of wavenumbers

solve_mode

solve_mode(
    mode_number: int, *args: Any, **kwargs: Any
) -> tuple[vcfdfield2_t, complex]

Wrapper around solve_modes() that solves for a single mode.

Parameters:

Name	Type	Description	Default
mode_number	int	0-indexed mode number to solve for	required
*args	Any	passed to solve_modes()	()
**kwargs	Any	passed to solve_modes()	{}

Returns:

Type	Description
tuple[vcfdfield2_t, complex]	(e_xy, wavenumber)

inner_product

inner_product(
    e1: vcfdfield2,
    h2: vcfdfield2,
    dxes: dx_lists2_t,
    prop_phase: float = 0,
    conj_h: bool = False,
    trapezoid: bool = False,
) -> complex

Compute the discrete waveguide overlap / Poynting inner product.

This is the 2D transverse integral corresponding to the time-averaged longitudinal Poynting flux,

\[ \frac{1}{2}\int (E_x H_y - E_y H_x) \, dx \, dy \]

with the Yee-grid staggering and optional propagation-phase adjustment used by the waveguide helpers in this module.

Parameters:

Name	Type	Description	Default
e1	vcfdfield2	Vectorized electric field, typically from exy2e(...) or normalized_fields_e(...).	required
h2	vcfdfield2	Vectorized magnetic field, typically from hxy2h(...), exy2h(...), or one of the normalization helpers.	required
dxes	dx_lists2_t	Two-dimensional Yee-grid spacing lists [dx_e, dx_h].	required
prop_phase	float	Phase advance over one propagation cell. This is used to shift the H field into the same longitudinal reference plane as the E field.	0
conj_h	bool	Whether to conjugate h2 before forming the overlap. Use True for the usual time-averaged power normalization.	False
trapezoid	bool	Whether to use trapezoidal quadrature instead of the default rectangular Yee-cell sum.	False

Returns:

Type	Description
complex	Complex overlap / longitudinal power integral.

meanas.fdfd.waveguide_3d

Tools for working with waveguide modes in 3D domains.

This module relies heavily on waveguide_2d and mostly just transforms its parameters into 2D equivalents and expands the results back into 3D.

The intended workflow is:

1. Select a single-cell slice normal to the propagation axis.
2. Solve the corresponding 2D mode problem with solve_mode(...).
3. Turn that mode into a one-sided source with compute_source(...).
4. Build an overlap window with compute_overlap_e(...) for port readout.

polarity is part of the public convention throughout this module:

• +1 means forward propagation toward increasing index along axis
• -1 means backward propagation toward decreasing index along axis

That same convention controls which side of the selected slice is used for the overlap window and how the expanded field is phased.

solve_mode

solve_mode(
    mode_number: int,
    omega: complex,
    dxes: dx_lists_t,
    axis: int,
    polarity: int,
    slices: Sequence[slice],
    epsilon: fdfield,
    mu: fdfield | None = None,
) -> dict[str, complex | NDArray[complexfloating]]

Given a 3D grid, selects a slice from the grid and attempts to solve for an eigenmode propagating through that slice.

Parameters:

Name	Type	Description	Default
mode_number	int	Number of the mode, 0-indexed	required
omega	complex	Angular frequency of the simulation	required
dxes	dx_lists_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types	required
axis	int	Propagation axis (0=x, 1=y, 2=z)	required
polarity	int	Propagation direction (+1 for +ve, -1 for -ve)	required
slices	Sequence[slice]	epsilon[tuple(slices)] is used to select the portion of the grid to use as the waveguide cross-section. slices[axis] must select exactly one item.	required
epsilon	fdfield	Dielectric constant	required
mu	fdfield | None	Magnetic permeability (default 1 everywhere)	None

Returns:

Type	Description
dict[str, complex | NDArray[complexfloating]]	Dictionary containing:
dict[str, complex | NDArray[complexfloating]]	E: full-grid electric field for the solved mode
dict[str, complex | NDArray[complexfloating]]	H: full-grid magnetic field for the solved mode
dict[str, complex | NDArray[complexfloating]]	wavenumber: propagation constant corrected for the discretized propagation axis
dict[str, complex | NDArray[complexfloating]]	wavenumber_2d: propagation constant of the reduced 2D eigenproblem

Notes

The returned fields are normalized through the waveguide_2d normalization convention before being expanded back to 3D.

compute_source

compute_source(
    E: cfdfield,
    wavenumber: complex,
    omega: complex,
    dxes: dx_lists_t,
    axis: int,
    polarity: int,
    slices: Sequence[slice],
    epsilon: fdfield,
    mu: fdfield | None = None,
) -> cfdfield_t

Given an eigenmode obtained by solve_mode, returns the current source distribution necessary to position a unidirectional source at the slice location.

Parameters:

Name	Type	Description	Default
E	cfdfield	E-field of the mode	required
wavenumber	complex	Wavenumber of the mode	required
omega	complex	Angular frequency of the simulation	required
dxes	dx_lists_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types	required
axis	int	Propagation axis (0=x, 1=y, 2=z)	required
polarity	int	Propagation direction (+1 for +ve, -1 for -ve)	required
slices	Sequence[slice]	epsilon[tuple(slices)] is used to select the portion of the grid to use as the waveguide cross-section. slices[axis] should select only one item.	required
mu	fdfield | None	Magnetic permeability (default 1 everywhere)	None

Returns:

Type	Description
cfdfield_t	J distribution for a one-sided electric-current source.

Notes

The source is built from the expanded mode field and a boundary-source operator. The resulting current is intended to be injected with the same sign convention used elsewhere in the package:

E -= dt * J / epsilon

compute_overlap_e

compute_overlap_e(
    E: cfdfield_t,
    wavenumber: complex,
    dxes: dx_lists_t,
    axis: int,
    polarity: int,
    slices: Sequence[slice],
    omega: float,
) -> cfdfield_t

Build an overlap field for projecting another 3D electric field onto a mode.

The returned field is intended for the discrete overlap expression

\[ \sum \mathrm{overlap\_e} \; E_\mathrm{other}^* \]

where the sum is over the full Yee-grid field storage.

The construction uses a two-cell window immediately upstream of the selected slice:

• for polarity=+1, the two cells just before slices[axis].start
• for polarity=-1, the two cells just after slices[axis].stop

The window is clipped to the simulation domain if necessary. A clipped but non-empty window raises RuntimeWarning; an empty window raises ValueError.

The derivation below assumes reflection symmetry and the standard waveguide overlap relation involving

\[ \int ((E \times H_\mathrm{mode}) + (E_\mathrm{mode} \times H)) \cdot dn. \]

E x H_mode + E_mode x H -> Ex Hmy - EyHmx + Emx Hy - Emy Hx (Z-prop) Ex we/B Emx + Ex i/B dy Hmz - Ey (-we/B Emy) - Ey i/B dx Hmz we/B (Ex Emx + Ey Emy) + i/B (Ex dy Hmz - Ey dx Hmz) we/B (Ex Emx + Ey Emy) + i/B (Ex dy (dx Emy - dy Emx) - Ey dx (dx Emy - dy Emx)) we/B (Ex Emx + Ey Emy) + i/B (Ex dy dx Emy - Ex dy dy Emx - Ey dx dx Emy - Ey dx dy Emx)

Ex j/wu (-jB Emx - dx Emz) - Ey j/wu (dy Emz + jB Emy) B/wu (Ex Emx + Ey Emy) - j/wu (Ex dx Emz + Ey dy Emz)

Parameters:

Name	Type	Description	Default
E	cfdfield_t	E-field of the mode	required
wavenumber	complex	Wavenumber of the mode	required
dxes	dx_lists_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types	required
axis	int	Propagation axis (0=x, 1=y, 2=z)	required
polarity	int	Propagation direction (+1 for +ve, -1 for -ve)	required
slices	Sequence[slice]	epsilon[tuple(slices)] is used to select the portion of the grid to use as the waveguide cross-section. slices[axis] should select only one item.	required

Returns:

Type	Description
cfdfield_t	overlap_e normalized so that numpy.sum(overlap_e * E.conj()) == 1
cfdfield_t	over the retained overlap window.

expand_e

expand_e(
    E: cfdfield,
    wavenumber: complex,
    dxes: dx_lists_t,
    axis: int,
    polarity: int,
    slices: Sequence[slice],
) -> cfdfield_t

Given an eigenmode obtained by solve_mode, expands the E-field from the 2D slice where the mode was calculated to the entire domain (along the propagation axis). This assumes the epsilon cross-section remains constant throughout the entire domain; it is up to the caller to truncate the expansion to any regions where it is valid.

Parameters:

Name	Type	Description	Default
E	cfdfield	E-field of the mode	required
wavenumber	complex	Wavenumber of the mode	required
dxes	dx_lists_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types	required
axis	int	Propagation axis (0=x, 1=y, 2=z)	required
polarity	int	Propagation direction (+1 for +ve, -1 for -ve)	required
slices	Sequence[slice]	epsilon[tuple(slices)] is used to select the portion of the grid to use as the waveguide cross-section. slices[axis] should select only one item.	required

Returns:

Type	Description
cfdfield_t	E, with the original field expanded along the specified axis.

Notes

This helper assumes that the waveguide cross-section remains constant along the propagation axis and applies the phase factor

\[ e^{-i \, \mathrm{polarity} \, wavenumber \, \Delta z} \]

to each copied slice.

meanas.fdfd.waveguide_cyl

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 coordinate-transformed discrete curl equations used here are

\[ \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, \\ \imath \beta H_y &= -\imath \omega T_b \epsilon_{rr} E_r - T_b \hat{\partial}_y H_z, \\ \imath \beta H_r &= \imath \omega T_a \epsilon_{yy} E_y - T_b T_a^{-1} \hat{\partial}_r (T_b H_z), \\ \imath \omega E_z &= T_a \epsilon_{zz}^{-1} \left(\hat{\partial}_r H_y - \hat{\partial}_y H_r\right). \end{aligned} \]

The first three equations are the cylindrical analogue of the straight-guide relations for H_r, H_y, and H_z. The next two are the metric-weighted versions of the straight-guide identities for \imath \beta H_y and \imath \beta H_r, and the last equation plays the same role as the longitudinal E_z reconstruction in waveguide_2d.

Following the same elimination steps as in waveguide_2d, apply \imath \beta \tilde{\partial}_r and \imath \beta \tilde{\partial}_y to the equation for E_z, substitute for \imath \beta H_r and \imath \beta H_y, and then eliminate H_z with

\[ H_z = \frac{1}{-\imath \omega \mu_{zz}} \left(\tilde{\partial}_r E_y - \tilde{\partial}_y E_r\right). \]

This yields the transverse electric 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}. \]

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, ...]]].

cylindrical_operator

cylindrical_operator(
    omega: float,
    dxes: dx_lists2_t,
    epsilon: vfdslice,
    rmin: float,
) -> sparse.sparray

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)

Parameters:

Name	Type	Description	Default
omega	float	The angular frequency of the system	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
epsilon	vfdslice	Vectorized dielectric constant grid	required
rmin	float	Radius at the left edge of the simulation domain (at minimum 'x')	required

Returns:

Type	Description
sparray	Sparse matrix representation of the operator

solve_modes

solve_modes(
    mode_numbers: Sequence[int],
    omega: float,
    dxes: dx_lists2_t,
    epsilon: vfdslice,
    rmin: float,
    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.

Parameters:

Name	Type	Description	Default
mode_numbers	Sequence[int]	Mode numbers to solve, 0-indexed.	required
omega	float	Angular frequency of the simulation	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types. The first coordinate is assumed to be r, the second is y.	required
epsilon	vfdslice	Dielectric constant	required
rmin	float	Radius of curvature for the simulation. This should be the minimum value of r within the simulation domain.	required

Returns:

Name	Type	Description
e_xys	NDArray[complex128]	NDArray of vfdfield_t specifying fields. First dimension is mode number.
angular_wavenumbers	NDArray[complex128]	list of wavenumbers in 1/rad units.

solve_mode

solve_mode(
    mode_number: int, *args: Any, **kwargs: Any
) -> tuple[vcfdslice, complex]

Wrapper around solve_modes() that solves for a single mode.

Parameters:

Name	Type	Description	Default
mode_number	int	0-indexed mode number to solve for	required
*args	Any	passed to solve_modes()	()
**kwargs	Any	passed to solve_modes()	{}

Returns:

Type	Description
tuple[vcfdslice, complex]	(e_xy, angular_wavenumber)

linear_wavenumbers

linear_wavenumbers(
    e_xys: list[vcfdfield2_t],
    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.

Parameters:

Name	Type	Description	Default
e_xys	list[vcfdfield2_t]	Vectorized mode fields with shape (num_modes, 2 * x *y)	required
angular_wavenumbers	ArrayLike	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	required
epsilon	vfdslice	Vectorized dielectric constant grid with shape (3, x, y)	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
rmin	float	Radius at the left edge of the simulation domain (at minimum 'x')	required

Returns:

Type	Description
NDArray[complex128]	NDArray containing the calculated linear (1/distance) wavenumbers

exy2h

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

Parameters:

Name	Type	Description	Default
angular_wavenumber	complex	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	required
omega	float	The angular frequency of the system	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
rmin	float	Radius at the left edge of the simulation domain (at minimum 'x')	required
epsilon	vfdslice	Vectorized dielectric constant grid	required
mu	vfdslice | None	Vectorized magnetic permeability grid (default 1 everywhere)	None

Returns:

Type	Description
sparray	Sparse matrix representing the operator.

exy2e

exy2e(
    angular_wavenumber: complex,
    omega: float,
    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.

Parameters:

Name	Type	Description	Default
angular_wavenumber	complex	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	required
omega	float	The angular frequency of the system	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
rmin	float	Radius at the left edge of the simulation domain (at minimum 'x')	required
epsilon	vfdslice	Vectorized dielectric constant grid	required

Returns:

Type	Description
sparray	Sparse matrix representing the operator.

e2h

e2h(
    angular_wavenumber: complex,
    omega: float,
    dxes: dx_lists2_t,
    rmin: float,
    mu: vfdslice | None = None,
) -> sparse.sparray

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} \]

Parameters:

Name	Type	Description	Default
angular_wavenumber	complex	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	required
omega	float	The angular frequency of the system	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
rmin	float	Radius at the left edge of the simulation domain (at minimum 'x')	required
mu	vfdslice | None	Vectorized magnetic permeability grid (default 1 everywhere)	None

Returns:

Type	Description
sparray	Sparse matrix representation of the operator.

dxes2T

dxes2T(
    dxes: dx_lists2_t, rmin: float
) -> tuple[NDArray[numpy.float64], NDArray[numpy.float64]]

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.

Parameters:

Name	Type	Description	Default
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
rmin	float	Radius at the left edge of the simulation domain (at minimum 'x')	required

Returns:

Type	Description
tuple[NDArray[float64], NDArray[float64]]	Sparse diagonal matrices (T_a, T_b).

normalized_fields_e

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]

Given a vector e_xy containing the vectorized E_r and E_y fields, returns normalized, vectorized E and H fields for the system.

Parameters:

Name	Type	Description	Default
e_xy	vcfdfield2	Vector containing E_r and E_y fields	required
angular_wavenumber	complex	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	required
omega	float	The angular frequency of the system	required
dxes	dx_lists2_t	Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types (2D)	required
rmin	float	Radius at the left edge of the simulation domain (at minimum 'x')	required
epsilon	vfdslice	Vectorized dielectric constant grid	required
mu	vfdslice | None	Vectorized magnetic permeability grid (default 1 everywhere)	None
prop_phase	float	Phase shift (dz * corrected_wavenumber) over 1 cell in propagation direction. Default 0 (continuous propagation direction, i.e. dz->0).	0

Returns:

Type	Description
vcfdslice_t	(e, h), where each field is vectorized, normalized,
vcfdslice_t	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.