Skip to content

meanas¶

meanas is a Python package for finite-difference electromagnetic simulation. It combines:

  • meanas.fdfd for frequency-domain operators, sources, waveguide modes, and SCPML
  • meanas.fdtd for Yee-grid timestepping, CPML, energy/flux accounting, and phasor extraction
  • meanas.fdmath for the shared discrete operators and derivations underneath both solvers

This documentation is built directly from the package docstrings. The API pages are the source of truth for the mathematical derivations and calling conventions.

Examples and API Map¶

For most users, the tracked examples under examples/ are the right entry point. They show the intended combinations of tools for solving complete problems.

The API pages are better read as a toolbox map and derivation reference:

  • Use the FDTD API for time-domain stepping, CPML, and phasor extraction.
  • Use the FDFD API for driven frequency-domain solves and sparse operator algebra.
  • Use the Waveguide API for mode solving, port sources, and overlap windows.
  • Use the fdmath API for the lower-level finite-difference operators and the shared discrete derivations underneath both solvers.

Build outputs¶

The docs build generates two HTML views from the same source:

  • a normal multi-page site
  • a print-oriented combined page under site/print_page/

If htmlark is installed, ./make_docs.sh also writes a fully inlined site/standalone.html.

API

API Overview¶

The package is documented directly from its docstrings. The most useful entry points are:

  • meanas: top-level package overview
  • eigensolvers: generic eigenvalue utilities used by the mode solvers
  • fdfd: frequency-domain operators, sources, PML, solvers, and far-field transforms
  • waveguides: straight, cylindrical, and 3D waveguide mode helpers
  • fdtd: timestepping, CPML, energy/flux helpers, and phasor extraction
  • fdmath: shared discrete operators, vectorization helpers, and derivation background

The waveguide and FDTD pages are the best places to start if you want the mathematical derivations rather than just the callable reference.

meanas¶

meanas ¶

Electromagnetic simulation tools

See the tracked examples for end-to-end workflows, and help(meanas) for the toolbox overview and API derivations.

eigensolvers¶

meanas.eigensolvers ¶

Solvers for eigenvalue / eigenvector problems

power_iteration ¶

power_iteration(
    operator: spmatrix,
    guess_vector: NDArray[complex128] | None = None,
    iterations: int = 20,
) -> tuple[complex, NDArray[numpy.complex128]]

Use power iteration to estimate the dominant eigenvector of a matrix.

Parameters:

Name Type Description Default
operator spmatrix

Matrix to analyze.

 required
guess_vector NDArray[complex128] | None

Starting point for the eigenvector. Default is a randomly chosen vector.

 None
iterations int

Number of iterations to perform. Default 20.

 20

Returns:

Type Description
tuple[complex, NDArray[complex128]]

(Largest-magnitude eigenvalue, Corresponding eigenvector estimate)

rayleigh_quotient_iteration ¶

rayleigh_quotient_iteration(
    operator: spmatrix | LinearOperator,
    guess_vector: NDArray[complex128],
    iterations: int = 40,
    tolerance: float = 1e-13,
    solver: Callable[..., NDArray[complex128]]
    | None = None,
) -> tuple[complex, NDArray[numpy.complex128]]

Use Rayleigh quotient iteration to refine an eigenvector guess.

Parameters:

Name Type Description Default
operator spmatrix | LinearOperator

Matrix to analyze.

 required
guess_vector NDArray[complex128]

Eigenvector to refine.

 required
iterations int

Maximum number of iterations to perform. Default 40.

 40
tolerance float

Stop iteration if (A - I*eigenvalue) @ v < num_vectors * tolerance, Default 1e-13.

 1e-13
solver Callable[..., NDArray[complex128]] | None

Solver function of the form x = solver(A, b). By default, use scipy.sparse.spsolve for sparse matrices and scipy.sparse.bicgstab for general LinearOperator instances.

 None

Returns:

Type Description
tuple[complex, NDArray[complex128]]

(eigenvalues, eigenvectors)

signed_eigensolve ¶

signed_eigensolve(
    operator: spmatrix | LinearOperator,
    how_many: int,
    negative: bool = False,
) -> tuple[
    NDArray[numpy.complex128], NDArray[numpy.complex128]
]

Find the largest-magnitude positive-only (or negative-only) eigenvalues and eigenvectors of the provided matrix.

Parameters:

Name Type Description Default
operator spmatrix | LinearOperator

Matrix to analyze.

 required
how_many int

How many eigenvalues to find.

 required
negative bool

Whether to find negative-only eigenvalues. Default False (positive only).

 False

Returns:

Type Description
NDArray[complex128]

(sorted list of eigenvalues, 2D ndarray of corresponding eigenvectors)
NDArray[complex128]

eigenvectors[:, k] corresponds to the k-th eigenvalue

fdfd¶

meanas.fdfd ¶

Tools for finite difference frequency-domain (FDFD) simulations and calculations.

These mostly involve picking a single frequency, then setting up and solving a matrix equation (Ax=b) or eigenvalue problem.

Submodules:

  • operators, functional: General FDFD problem setup.
  • solvers: Solver interface and reference implementation.
  • scpml: Stretched-coordinate perfectly matched layer (SCPML) boundary conditions.
  • waveguide_2d: Operators and mode-solver for waveguides with constant cross-section.
  • waveguide_3d: Functions for transforming waveguide_2d results into 3D, including mode-source and overlap-window construction.
  • farfield, bloch, eme: specialized helper modules for near/far transforms, Bloch-periodic problems, and eigenmode expansion.

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

From the "Frequency domain" section of meanas.fdmath, we have

\[ \begin{aligned} \tilde{E}_{l, \vec{r}} &= \tilde{E}_{\vec{r}} e^{-\imath \omega l \Delta_t} \\ \tilde{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &= \tilde{H}_{\vec{r} + \frac{1}{2}} e^{-\imath \omega (l - \frac{1}{2}) \Delta_t} \\ \tilde{J}_{l, \vec{r}} &= \tilde{J}_{\vec{r}} e^{-\imath \omega (l - \frac{1}{2}) \Delta_t} \\ \tilde{M}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &= \tilde{M}_{\vec{r} + \frac{1}{2}} e^{-\imath \omega l \Delta_t} \\ \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{\vec{r}}) -\Omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} &= -\imath \Omega \tilde{J}_{\vec{r}} e^{\imath \omega \Delta_t / 2} \\ \Omega &= 2 \sin(\omega \Delta_t / 2) / \Delta_t \end{aligned} \]

resulting in

\[ \begin{aligned} \tilde{\partial}_t &\Rightarrow -\imath \Omega e^{-\imath \omega \Delta_t / 2}\\ \hat{\partial}_t &\Rightarrow -\imath \Omega e^{ \imath \omega \Delta_t / 2}\\ \end{aligned} \]

Maxwell's equations are then

\[ \begin{aligned} \tilde{\nabla} \times \tilde{E}_{\vec{r}} &= \imath \Omega e^{-\imath \omega \Delta_t / 2} \hat{B}_{\vec{r} + \frac{1}{2}} - \hat{M}_{\vec{r} + \frac{1}{2}} \\ \hat{\nabla} \times \hat{H}_{\vec{r} + \frac{1}{2}} &= -\imath \Omega e^{ \imath \omega \Delta_t / 2} \tilde{D}_{\vec{r}} + \tilde{J}_{\vec{r}} \\ \tilde{\nabla} \cdot \hat{B}_{\vec{r} + \frac{1}{2}} &= 0 \\ \hat{\nabla} \cdot \tilde{D}_{\vec{r}} &= \rho_{\vec{r}} \end{aligned} \]

With \(\Delta_t \to 0\), this simplifies to

\[ \begin{aligned} \tilde{E}_{l, \vec{r}} &\to \tilde{E}_{\vec{r}} \\ \tilde{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &\to \tilde{H}_{\vec{r} + \frac{1}{2}} \\ \tilde{J}_{l, \vec{r}} &\to \tilde{J}_{\vec{r}} \\ \tilde{M}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &\to \tilde{M}_{\vec{r} + \frac{1}{2}} \\ \Omega &\to \omega \\ \tilde{\partial}_t &\to -\imath \omega \\ \hat{\partial}_t &\to -\imath \omega \\ \end{aligned} \]

and then

\[ \begin{aligned} \tilde{\nabla} \times \tilde{E}_{\vec{r}} &= \imath \omega \hat{B}_{\vec{r} + \frac{1}{2}} - \hat{M}_{\vec{r} + \frac{1}{2}} \\ \hat{\nabla} \times \hat{H}_{\vec{r} + \frac{1}{2}} &= -\imath \omega \tilde{D}_{\vec{r}} + \tilde{J}_{\vec{r}} \\ \end{aligned} \]
\[ \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{\vec{r}}) -\omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} = -\imath \omega \tilde{J}_{\vec{r}} \\ \]

Core operator layers¶

meanas.fdfd.functional ¶

Functional versions of many FDFD operators. These can be useful for performing FDFD calculations without needing to construct large matrices in memory.

The functions generated here expect cfdfield_t inputs with shape (3, X, Y, Z), e.g. E = [E_x, E_y, E_z] where each (complex) component has shape (X, Y, Z)

e_full ¶

e_full(
    omega: complex,
    dxes: dx_lists_t,
    epsilon: fdfield,
    mu: fdfield | None = None,
) -> cfdfield_updater_t

Wave operator for use with E-field. See operators.e_full for details.

Parameters:

Name Type Description Default
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
epsilon fdfield

Dielectric constant

 required
mu fdfield | None

Magnetic permeability (default 1 everywhere)

 None

Returns:

Type Description
cfdfield_updater_t

Function f implementing the wave operator
cfdfield_updater_t

f(E) -> -i * omega * J

eh_full ¶

eh_full(
    omega: complex,
    dxes: dx_lists_t,
    epsilon: fdfield,
    mu: fdfield | None = None,
) -> Callable[
    [cfdfield, cfdfield], tuple[cfdfield_t, cfdfield_t]
]

Wave operator for full (both E and H) field representation. See operators.eh_full.

Parameters:

Name Type Description Default
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
epsilon fdfield

Dielectric constant

 required
mu fdfield | None

Magnetic permeability (default 1 everywhere)

 None

Returns:

Type Description
Callable[[cfdfield, cfdfield], tuple[cfdfield_t, cfdfield_t]]

Function f implementing the wave operator
Callable[[cfdfield, cfdfield], tuple[cfdfield_t, cfdfield_t]]

f(E, H) -> (J, -M)

e2h ¶

e2h(
    omega: complex,
    dxes: dx_lists_t,
    mu: fdfield | None = None,
) -> cfdfield_updater_t

Utility operator for converting the E field into the H field. For use with e_full -- assumes that there is no magnetic current M.

Parameters:

Name Type Description Default
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
mu fdfield | None

Magnetic permeability (default 1 everywhere)

 None

Returns:

Type Description
cfdfield_updater_t

Function f for converting E to H,
cfdfield_updater_t

f(E) -> H

m2j ¶

m2j(
    omega: complex,
    dxes: dx_lists_t,
    mu: fdfield | None = None,
) -> cfdfield_updater_t

Utility operator for converting magnetic current M distribution into equivalent electric current distribution J. For use with e.g. e_full.

Parameters:

Name Type Description Default
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
mu fdfield | None

Magnetic permeability (default 1 everywhere)

 None

Returns:

Type Description
cfdfield_updater_t

Function f for converting M to J,
cfdfield_updater_t

f(M) -> J

e_tfsf_source ¶

e_tfsf_source(
    TF_region: fdfield,
    omega: complex,
    dxes: dx_lists_t,
    epsilon: fdfield,
    mu: fdfield | None = None,
) -> cfdfield_updater_t

Operator that turns an E-field distribution into a total-field/scattered-field (TFSF) source.

If A is the full wave operator from e_full(...) and Q is the diagonal mask selecting the total-field region, then the TFSF source is the commutator

\[ \frac{A Q - Q A}{-i \omega} E. \]

This vanishes in the interior of the total-field and scattered-field regions and is supported only at their shared boundary, where the mask discontinuity makes A and Q fail to commute. The returned current is therefore the distributed source needed to inject the desired total field without also forcing the scattered-field region.

Parameters:

Name Type Description Default
TF_region fdfield

mask which is set to 1 in the total-field region, and 0 elsewhere (i.e. in the scattered-field region). Should have the same shape as the simulation grid, e.g. epsilon[0].shape.

 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
epsilon fdfield

Dielectric constant distribution

 required
mu fdfield | None

Magnetic permeability (default 1 everywhere)

 None

Returns:

Type Description
cfdfield_updater_t

Function f which takes an E field and returns a current distribution,
cfdfield_updater_t

f(E) -> J

poynting_e_cross_h ¶

poynting_e_cross_h(
    dxes: dx_lists_t,
) -> Callable[[cfdfield, cfdfield], cfdfield_t]

Generates a function that takes the single-frequency E and H fields and calculates the cross product E x H = \(E \times H\) as required for the Poynting vector, \(S = E \times H\).

On the Yee grid, the electric and magnetic components are not stored at the same locations. This helper therefore applies the same one-cell electric-field shifts used by the sparse operators.poynting_e_cross(...) construction so that the discrete cross product matches the face-centered energy flux used in meanas.fdtd.energy.poynting(...).

Note

This function also shifts the input E field by one cell as required for computing the Poynting cross product (see meanas.fdfd module docs).

Note

If E and H are peak amplitudes as assumed elsewhere in this code, the time-average of the poynting vector is <S> = Re(S)/2 = Re(E x H*) / 2. The factor of 1/2 can be omitted if root-mean-square quantities are used instead.

Parameters:

Name Type Description Default
dxes dx_lists_t

Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types

 required

Returns:

Type Description
Callable[[cfdfield, cfdfield], cfdfield_t]

Function f that returns the staggered-grid cross product E \times H.
Callable[[cfdfield, cfdfield], cfdfield_t]

For time-average power, call it as f(E, H.conj()) and take Re(...) / 2.

meanas.fdfd.operators ¶

Sparse matrix operators for use with electromagnetic wave equations.

These functions return sparse-matrix (scipy.sparse.sparray) representations of a variety of operators, intended for use with E and H fields vectorized using the meanas.fdmath.vectorization.vec() and meanas.fdmath.vectorization.unvec() functions.

E- and H-field values are defined on a Yee cell; epsilon values should be calculated for cells centered at each E component (mu at each H component).

Many of these functions require a dxes parameter, of type dx_lists_t; see the meanas.fdmath.types submodule for details.

The following operators are included:

  • E-only wave operator
  • H-only wave operator
  • EH wave operator
  • Curl for use with E, H fields
  • E to H conversion
  • M to J conversion
  • Poynting cross products
  • Circular shifts
  • Discrete derivatives
  • Averaging operators
  • Cross product matrices

e_full ¶

e_full(
    omega: complex,
    dxes: dx_lists_t,
    epsilon: vfdfield | vcfdfield,
    mu: vfdfield | None = None,
    pec: vfdfield | None = None,
    pmc: vfdfield | None = None,
) -> sparse.sparray

Wave operator $$ \nabla \times (\frac{1}{\mu} \nabla \times) - \Omega^2 \epsilon $$

del x (1/mu * del x) - omega**2 * epsilon

for use with the E-field, with wave equation $$ (\nabla \times (\frac{1}{\mu} \nabla \times) - \Omega^2 \epsilon) E = -\imath \omega J $$

(del x (1/mu * del x) - omega**2 * epsilon) E = -i * omega * J

To make this matrix symmetric, use the preconditioners from e_full_preconditioners().

Parameters:

Name Type Description Default
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
epsilon vfdfield | vcfdfield

Vectorized dielectric constant

 required
mu vfdfield | None

Vectorized magnetic permeability (default 1 everywhere).

 None
pec vfdfield | None

Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted as containing a perfect electrical conductor (PEC). The PEC is applied per-field-component (i.e. pec.size == epsilon.size)

 None
pmc vfdfield | None

Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted as containing a perfect magnetic conductor (PMC). The PMC is applied per-field-component (i.e. pmc.size == epsilon.size)

 None

Returns:

Type Description
sparray

Sparse matrix containing the wave operator.

e_full_preconditioners ¶

e_full_preconditioners(
    dxes: dx_lists_t,
) -> tuple[sparse.sparray, sparse.sparray]

Left and right preconditioners (Pl, Pr) for symmetrizing the e_full wave operator.

The preconditioned matrix A_symm = (Pl @ A @ Pr) is complex-symmetric (non-Hermitian unless there is no loss or PMLs).

The preconditioner matrices are diagonal and complex, with Pr = 1 / Pl

Parameters:

Name Type Description Default
dxes dx_lists_t

Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types

 required

Returns:

Type Description
tuple[sparray, sparray]

Preconditioner matrices (Pl, Pr).

h_full ¶

h_full(
    omega: complex,
    dxes: dx_lists_t,
    epsilon: vfdfield,
    mu: vfdfield | None = None,
    pec: vfdfield | None = None,
    pmc: vfdfield | None = None,
) -> sparse.sparray

Wave operator $$ \nabla \times (\frac{1}{\epsilon} \nabla \times) - \omega^2 \mu $$

del x (1/epsilon * del x) - omega**2 * mu

for use with the H-field, with wave equation $$ (\nabla \times (\frac{1}{\epsilon} \nabla \times) - \omega^2 \mu) E = \imath \omega M $$

(del x (1/epsilon * del x) - omega**2 * mu) E = i * omega * M

Parameters:

Name Type Description Default
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
epsilon vfdfield

Vectorized dielectric constant

 required
mu vfdfield | None

Vectorized magnetic permeability (default 1 everywhere)

 None
pec vfdfield | None

Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted as containing a perfect electrical conductor (PEC). The PEC is applied per-field-component (i.e. pec.size == epsilon.size)

 None
pmc vfdfield | None

Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted as containing a perfect magnetic conductor (PMC). The PMC is applied per-field-component (i.e. pmc.size == epsilon.size)

 None

Returns:

Type Description
sparray

Sparse matrix containing the wave operator.

eh_full ¶

eh_full(
    omega: complex,
    dxes: dx_lists_t,
    epsilon: vfdfield,
    mu: vfdfield | None = None,
    pec: vfdfield | None = None,
    pmc: vfdfield | None = None,
) -> sparse.sparray

Wave operator for [E, H] field representation. This operator implements Maxwell's equations without cancelling out either E or H. The operator is $$ \begin{bmatrix} -\imath \omega \epsilon & \nabla \times \ \nabla \times & \imath \omega \mu \end{bmatrix} $$

[[-i * omega * epsilon,  del x         ],
 [del x,                 i * omega * mu]]

for use with a field vector of the form cat(vec(E), vec(H)): $$ \begin{bmatrix} -\imath \omega \epsilon & \nabla \times \ \nabla \times & \imath \omega \mu \end{bmatrix} \begin{bmatrix} E \ H \end{bmatrix} = \begin{bmatrix} J \ -M \end{bmatrix} $$

Parameters:

Name Type Description Default
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
epsilon vfdfield

Vectorized dielectric constant

 required
mu vfdfield | None

Vectorized magnetic permeability (default 1 everywhere)

 None
pec vfdfield | None

Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted as containing a perfect electrical conductor (PEC). The PEC is applied per-field-component (i.e. pec.size == epsilon.size)

 None
pmc vfdfield | None

Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted as containing a perfect magnetic conductor (PMC). The PMC is applied per-field-component (i.e. pmc.size == epsilon.size)

 None

Returns:

Type Description
sparray

Sparse matrix containing the wave operator.

e2h ¶

e2h(
    omega: complex,
    dxes: dx_lists_t,
    mu: vfdfield | None = None,
    pmc: vfdfield | None = None,
) -> sparse.sparray

Utility operator for converting the E field into the H field. For use with e_full() -- assumes that there is no magnetic current M.

Parameters:

Name Type Description Default
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
mu vfdfield | None

Vectorized magnetic permeability (default 1 everywhere)

 None
pmc vfdfield | None

Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted as containing a perfect magnetic conductor (PMC). The PMC is applied per-field-component (i.e. pmc.size == epsilon.size)

 None

Returns:

Type Description
sparray

Sparse matrix for converting E to H.

m2j ¶

m2j(
    omega: complex,
    dxes: dx_lists_t,
    mu: vfdfield | None = None,
) -> sparse.sparray

Operator for converting a magnetic current M into an electric current J. For use with eg. e_full().

Parameters:

Name Type Description Default
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
mu vfdfield | None

Vectorized magnetic permeability (default 1 everywhere)

 None

Returns:

Type Description
sparray

Sparse matrix for converting M to J.

poynting_e_cross ¶

poynting_e_cross(
    e: vcfdfield, dxes: dx_lists_t
) -> sparse.sparray

Operator for computing the staggered-grid (E \times) part of the Poynting vector.

On the Yee grid the E and H components live on different edges, so the electric field must be shifted by one cell in the appropriate direction before forming the discrete cross product. This sparse operator encodes that shifted cross product directly and is the matrix equivalent of functional.poynting_e_cross_h(...).

Parameters:

Name Type Description Default
e vcfdfield

Vectorized E-field for the ExH cross product

 required
dxes dx_lists_t

Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types

 required

Returns:

Type Description
sparray

Sparse matrix containing the (E \times) part of the staggered Poynting
sparray

cross product.

poynting_h_cross ¶

poynting_h_cross(
    h: vcfdfield, dxes: dx_lists_t
) -> sparse.sparray

Operator for computing the staggered-grid (H \times) part of the Poynting vector.

Together with poynting_e_cross(...), this provides the matrix form of the Yee-grid cross product used in the flux helpers. The two are related by the usual antisymmetry of the cross product,

\[ H \times E = -(E \times H), \]

once the same staggered field placement is used on both sides.

Parameters:

Name Type Description Default
h vcfdfield

Vectorized H-field for the HxE cross product

 required
dxes dx_lists_t

Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types

 required

Returns:

Type Description
sparray

Sparse matrix containing the (H \times) part of the staggered Poynting
sparray

cross product.

e_tfsf_source ¶

e_tfsf_source(
    TF_region: vfdfield,
    omega: complex,
    dxes: dx_lists_t,
    epsilon: vfdfield,
    mu: vfdfield | None = None,
) -> sparse.sparray

Operator that turns a desired E-field distribution into a total-field/scattered-field (TFSF) source.

Let A be the full wave operator from e_full(...), and let Q = \mathrm{diag}(TF_region) be the projector onto the total-field region. Then the TFSF current operator is the commutator

\[ \frac{A Q - Q A}{-i \omega}. \]

Inside regions where Q is locally constant, A and Q commute and the source vanishes. Only cells at the TF/SF boundary contribute nonzero current, which is exactly the desired distributed source for injecting the chosen field into the total-field region without directly forcing the scattered-field region.

Parameters:

Name Type Description Default
TF_region vfdfield

Mask, which is set to 1 inside the total-field region and 0 in the scattered-field region

 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
epsilon vfdfield

Vectorized dielectric constant

 required
mu vfdfield | None

Vectorized magnetic permeability (default 1 everywhere).

 None

Returns:

Type Description
sparray

Sparse matrix that turns an E-field into a current (J) distribution.

e_boundary_source ¶

e_boundary_source(
    mask: vfdfield,
    omega: complex,
    dxes: dx_lists_t,
    epsilon: vfdfield,
    mu: vfdfield | None = None,
    periodic_mask_edges: bool = False,
) -> sparse.sparray

Operator that turns an E-field distrubtion into a current (J) distribution along the edges (external and internal) of the provided mask. This is just an e_tfsf_source() with an additional masking step.

Equivalently, this helper first constructs the TFSF commutator source for the full mask and then zeroes out all cells except the mask boundary. The boundary is defined as the set of cells whose mask value changes under a one-cell shift in any Cartesian direction. With periodic_mask_edges=False the shifts mirror at the domain boundary; with True they wrap periodically.

Parameters:

Name Type Description Default
mask vfdfield

The current distribution is generated at the edges of the mask, i.e. any points where shifting the mask by one cell in any direction would change its value.

 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
epsilon vfdfield

Vectorized dielectric constant

 required
mu vfdfield | None

Vectorized magnetic permeability (default 1 everywhere).

 None

Returns:

Type Description
sparray

Sparse matrix that turns an E-field into a current (J) distribution.

meanas.fdfd.solvers ¶

Solvers and solver interface for FDFD problems.

generic ¶

generic(
    omega: complex,
    dxes: dx_lists_t,
    J: vcfdfield,
    epsilon: vfdfield,
    mu: vfdfield | None = None,
    *,
    pec: vfdfield | None = None,
    pmc: vfdfield | None = None,
    adjoint: bool = False,
    matrix_solver: Callable[..., ArrayLike] = _scipy_qmr,
    matrix_solver_opts: dict[str, Any] | None = None,
    E_guess: vcfdfield | None = None,
) -> vcfdfield_t

Conjugate gradient FDFD solver using CSR sparse matrices.

All ndarray arguments should be 1D arrays, as returned by meanas.fdmath.vectorization.vec().

Parameters:

Name Type Description Default
omega complex

Complex frequency to solve at.

 required
dxes dx_lists_t

[[dx_e, dy_e, dz_e], [dx_h, dy_h, dz_h]] (complex cell sizes) as discussed in meanas.fdmath.types

 required
J vcfdfield

Electric current distribution (at E-field locations)

 required
epsilon vfdfield

Dielectric constant distribution (at E-field locations)

 required
mu vfdfield | None

Magnetic permeability distribution (at H-field locations)

 None
pec vfdfield | None

Perfect electric conductor distribution (at E-field locations; non-zero value indicates PEC is present)

 None
pmc vfdfield | None

Perfect magnetic conductor distribution (at H-field locations; non-zero value indicates PMC is present)

 None
adjoint bool

If true, solves the adjoint problem.

 False
matrix_solver Callable[..., ArrayLike]

Called as matrix_solver(A, b, **matrix_solver_opts) -> x, where A: scipy.sparse.csr_array; b: ArrayLike; x: ArrayLike; Default is a wrapped version of scipy.sparse.linalg.qmr() which doesn't return convergence info and logs the residual every 100 iterations.

 _scipy_qmr
matrix_solver_opts dict[str, Any] | None

Passed as kwargs to matrix_solver(...)

 None
E_guess vcfdfield | None

Guess at the solution E-field. matrix_solver must accept an x0 argument with the same purpose.

 None

Returns:

Type Description
vcfdfield_t

E-field which solves the system.

meanas.fdfd.scpml ¶

Functions for creating stretched coordinate perfectly matched layer (PML) absorbers.

s_function_t module-attribute ¶

s_function_t = Callable[
    [NDArray[float64]], NDArray[float64]
]

Typedef for s-functions, see prepare_s_function()

prepare_s_function ¶

prepare_s_function(
    ln_R: float = -16, m: float = 4
) -> s_function_t

Create an s_function to pass to the SCPML functions. This is used when you would like to customize the PML parameters.

Parameters:

Name Type Description Default
ln_R float

Natural logarithm of the desired reflectance

 -16
m float

Polynomial order for the PML (imaginary part increases as distance ** m)

 4

Returns:

Type Description
s_function_t

An s_function, which takes an ndarray (distances) and returns an ndarray (complex part
s_function_t

of the cell width; needs to be divided by sqrt(epilon_effective) * real(omega))
s_function_t

before use.

uniform_grid_scpml ¶

uniform_grid_scpml(
    shape: Sequence[int],
    thicknesses: Sequence[int],
    omega: float,
    epsilon_effective: float = 1.0,
    s_function: s_function_t | None = None,
) -> list[list[NDArray[numpy.float64]]]

Create dx arrays for a uniform grid with a cell width of 1 and a pml.

If you want something more fine-grained, check out stretch_with_scpml(...).

Parameters:

Name Type Description Default
shape Sequence[int]

Shape of the grid, including the PMLs (which are 2*thicknesses thick)

 required
thicknesses Sequence[int]

[th_x, th_y, th_z] Thickness of the PML in each direction. Both polarities are added. Each th_ of pml is applied twice, once on each edge of the grid along the given axis. th_* may be zero, in which case no pml is added.

 required
omega float

Angular frequency for the simulation

 required
epsilon_effective float

Effective epsilon of the PML. Match this to the material at the edge of your grid. Default 1.

 1.0
s_function s_function_t | None

created by prepare_s_function(...), allowing customization of pml parameters. Default uses prepare_s_function() with no parameters.

 None

Returns:

Type Description
list[list[NDArray[float64]]]

Complex cell widths (dx_lists_mut) as discussed in meanas.fdmath.types.

stretch_with_scpml ¶

stretch_with_scpml(
    dxes: list[list[NDArray[float64]]],
    axis: int,
    polarity: int,
    omega: float,
    epsilon_effective: float = 1.0,
    thickness: int = 10,
    s_function: s_function_t | None = None,
) -> list[list[NDArray[numpy.float64]]]

Stretch dxes to contain a stretched-coordinate PML (SCPML) in one direction along one axis.

Parameters:

Name Type Description Default
dxes list[list[NDArray[float64]]]

Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types

 required
axis int

axis to stretch (0=x, 1=y, 2=z)

 required
polarity int

direction to stretch (-1 for -ve, +1 for +ve)

 required
omega float

Angular frequency for the simulation

 required
epsilon_effective float

Effective epsilon of the PML. Match this to the material at the edge of your grid. Default 1.

 1.0
thickness int

number of cells to use for pml (default 10)

 10
s_function s_function_t | None

Created by prepare_s_function(...), allowing customization of pml parameters. Default uses prepare_s_function() with no parameters.

 None

Returns:

Type Description
list[list[NDArray[float64]]]

Complex cell widths (dx_lists_mut) as discussed in meanas.fdmath.types.
list[list[NDArray[float64]]]

Multiple calls to this function may be necessary if multiple absorpbing boundaries are needed.

meanas.fdfd.farfield ¶

Functions for performing near-to-farfield transformation (and the reverse).

near_to_farfield ¶

near_to_farfield(
    E_near: cfdfield_t,
    H_near: cfdfield_t,
    dx: float,
    dy: float,
    padded_size: list[int] | int | None = None,
) -> dict[str, Any]

Compute the farfield, i.e. the distribution of the fields after propagation through several wavelengths of uniform medium.

The input fields should be complex phasors.

Parameters:

Name Type Description Default
E_near cfdfield_t

List of 2 ndarrays containing the 2D phasor field slices for the transverse E fields (e.g. [Ex, Ey] for calculating the farfield toward the z-direction).

 required
H_near cfdfield_t

List of 2 ndarrays containing the 2D phasor field slices for the transverse H fields (e.g. [Hx, hy] for calculating the farfield towrad the z-direction).

 required
dx float

Cell size along x-dimension, in units of wavelength.

 required
dy float

Cell size along y-dimension, in units of wavelength.

 required
padded_size list[int] | int | None

Shape of the output. A single integer n will be expanded to (n, n). Powers of 2 are most efficient for FFT computation. Default is the smallest power of 2 larger than the input, for each axis.

 None

Returns:

Type Description
dict[str, Any]

Dict with keys
dict[str, Any]
  • E_far: Normalized E-field farfield; multiply by (i k exp(-i k r) / (4 pi r)) to get the actual field value.
dict[str, Any]
  • H_far: Normalized H-field farfield; multiply by (i k exp(-i k r) / (4 pi r)) to get the actual field value.
dict[str, Any]
  • kx, ky: Wavevector values corresponding to the x- and y- axes in E_far and H_far, normalized to wavelength (dimensionless).
dict[str, Any]
  • dkx, dky: step size for kx and ky, normalized to wavelength.
dict[str, Any]
  • theta: arctan2(ky, kx) corresponding to each (kx, ky). This is the angle in the x-y plane, counterclockwise from above, starting from +x.
dict[str, Any]
  • phi: arccos(kz / k) corresponding to each (kx, ky). This is the angle away from +z.

far_to_nearfield ¶

far_to_nearfield(
    E_far: cfdfield_t,
    H_far: cfdfield_t,
    dkx: float,
    dky: float,
    padded_size: list[int] | int | None = None,
) -> dict[str, Any]

Compute the farfield, i.e. the distribution of the fields after propagation through several wavelengths of uniform medium.

The input fields should be complex phasors.

Parameters:

Name Type Description Default
E_far cfdfield_t

List of 2 ndarrays containing the 2D phasor field slices for the transverse E fields (e.g. [Ex, Ey] for calculating the nearfield toward the z-direction). Fields should be normalized so that E_far = E_far_actual / (i k exp(-i k r) / (4 pi r))

 required
H_far cfdfield_t

List of 2 ndarrays containing the 2D phasor field slices for the transverse H fields (e.g. [Hx, hy] for calculating the nearfield toward the z-direction). Fields should be normalized so that H_far = H_far_actual / (i k exp(-i k r) / (4 pi r))

 required
dkx float

kx discretization, in units of wavelength.

 required
dky float

ky discretization, in units of wavelength.

 required
padded_size list[int] | int | None

Shape of the output. A single integer n will be expanded to (n, n). Powers of 2 are most efficient for FFT computation. Default is the smallest power of 2 larger than the input, for each axis.

 None

Returns:

Type Description
dict[str, Any]

Dict with keys
dict[str, Any]
  • E: E-field nearfield
dict[str, Any]
  • H: H-field nearfield
dict[str, Any]
  • dx, dy: spatial discretization, normalized to wavelength (dimensionless)

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_number

Number of the mode, 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 E_r - T_b^{-1} \tilde{\partial}} H_y &= -\imath \beta T_b^{-1r (T_a E_z), \ -\imath \omega \mu_y E_r, \end{aligned} $$} H_z &= \tilde{\partial}_r E_y - \tilde{\partial

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.

fdtd¶

meanas.fdtd ¶

Utilities for running finite-difference time-domain (FDTD) simulations

See the discussion of Maxwell's Equations in meanas.fdmath for basic mathematical background.

Timestep¶

From the discussion of "Plane waves and the Dispersion relation" in meanas.fdmath, we have

\[ c^2 \Delta_t^2 = \frac{\Delta_t^2}{\mu \epsilon} < 1/(\frac{1}{\Delta_x^2} + \frac{1}{\Delta_y^2} + \frac{1}{\Delta_z^2}) \]

or, if \(\Delta_x = \Delta_y = \Delta_z\), then \(c \Delta_t < \frac{\Delta_x}{\sqrt{3}}\).

Based on this, we can set

dt = sqrt(mu.min() * epsilon.min()) / sqrt(1/dx_min**2 + 1/dy_min**2 + 1/dz_min**2)

The dx_min, dy_min, dz_min should be the minimum value across both the base and derived grids.

Poynting Vector and Energy Conservation¶

Let

\[ \begin{aligned} \tilde{S}_{l, l', \vec{r}} &=& &\tilde{E}_{l, \vec{r}} \otimes \hat{H}_{l', \vec{r} + \frac{1}{2}} \\ &=& &\vec{x} (\tilde{E}^y_{l,m+1,n,p} \hat{H}^z_{l',\vec{r} + \frac{1}{2}} - \tilde{E}^z_{l,m+1,n,p} \hat{H}^y_{l', \vec{r} + \frac{1}{2}}) \\ & &+ &\vec{y} (\tilde{E}^z_{l,m,n+1,p} \hat{H}^x_{l',\vec{r} + \frac{1}{2}} - \tilde{E}^x_{l,m,n+1,p} \hat{H}^z_{l', \vec{r} + \frac{1}{2}}) \\ & &+ &\vec{z} (\tilde{E}^x_{l,m,n,p+1} \hat{H}^y_{l',\vec{r} + \frac{1}{2}} - \tilde{E}^y_{l,m,n,p+1} \hat{H}^z_{l', \vec{r} + \frac{1}{2}}) \end{aligned} \]

where \(\vec{r} = (m, n, p)\) and \(\otimes\) is a modified cross product in which the \(\tilde{E}\) terms are shifted as indicated.

By taking the divergence and rearranging terms, we can show that

\[ \begin{aligned} \hat{\nabla} \cdot \tilde{S}_{l, l', \vec{r}} &= \hat{\nabla} \cdot (\tilde{E}_{l, \vec{r}} \otimes \hat{H}_{l', \vec{r} + \frac{1}{2}}) \\ &= \hat{H}_{l', \vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l, \vec{r}} - \tilde{E}_{l, \vec{r}} \cdot \hat{\nabla} \times \hat{H}_{l', \vec{r} + \frac{1}{2}} \\ &= \hat{H}_{l', \vec{r} + \frac{1}{2}} \cdot (-\tilde{\partial}_t \mu_{\vec{r} + \frac{1}{2}} \hat{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} - \hat{M}_{l, \vec{r} + \frac{1}{2}}) - \tilde{E}_{l, \vec{r}} \cdot (\hat{\partial}_t \tilde{\epsilon}_{\vec{r}} \tilde{E}_{l'+\frac{1}{2}, \vec{r}} + \tilde{J}_{l', \vec{r}}) \\ &= \hat{H}_{l'} \cdot (-\mu / \Delta_t)(\hat{H}_{l + \frac{1}{2}} - \hat{H}_{l - \frac{1}{2}}) - \tilde{E}_l \cdot (\epsilon / \Delta_t )(\tilde{E}_{l'+\frac{1}{2}} - \tilde{E}_{l'-\frac{1}{2}}) - \hat{H}_{l'} \cdot \hat{M}_{l} - \tilde{E}_l \cdot \tilde{J}_{l'} \\ \end{aligned} \]

where in the last line the spatial subscripts have been dropped to emphasize the time subscripts \(l, l'\), i.e.

\[ \begin{aligned} \tilde{E}_l &= \tilde{E}_{l, \vec{r}} \\ \hat{H}_l &= \tilde{H}_{l, \vec{r} + \frac{1}{2}} \\ \tilde{\epsilon} &= \tilde{\epsilon}_{\vec{r}} \\ \end{aligned} \]

etc. For \(l' = l + \frac{1}{2}\) we get

\[ \begin{aligned} \hat{\nabla} \cdot \tilde{S}_{l, l + \frac{1}{2}} &= \hat{H}_{l + \frac{1}{2}} \cdot (-\mu / \Delta_t)(\hat{H}_{l + \frac{1}{2}} - \hat{H}_{l - \frac{1}{2}}) - \tilde{E}_l \cdot (\epsilon / \Delta_t)(\tilde{E}_{l+1} - \tilde{E}_l) - \hat{H}_{l'} \cdot \hat{M}_l - \tilde{E}_l \cdot \tilde{J}_{l + \frac{1}{2}} \\ &= (-\mu / \Delta_t)(\hat{H}^2_{l + \frac{1}{2}} - \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}}) - (\epsilon / \Delta_t)(\tilde{E}_{l+1} \cdot \tilde{E}_l - \tilde{E}^2_l) - \hat{H}_{l'} \cdot \hat{M}_l - \tilde{E}_l \cdot \tilde{J}_{l + \frac{1}{2}} \\ &= -(\mu \hat{H}^2_{l + \frac{1}{2}} +\epsilon \tilde{E}_{l+1} \cdot \tilde{E}_l) / \Delta_t \\ +(\mu \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}} +\epsilon \tilde{E}^2_l) / \Delta_t \\ - \hat{H}_{l+\frac{1}{2}} \cdot \hat{M}_l \\ - \tilde{E}_l \cdot \tilde{J}_{l+\frac{1}{2}} \\ \end{aligned} \]

and for \(l' = l - \frac{1}{2}\),

\[ \begin{aligned} \hat{\nabla} \cdot \tilde{S}_{l, l - \frac{1}{2}} &= (\mu \hat{H}^2_{l - \frac{1}{2}} +\epsilon \tilde{E}_{l-1} \cdot \tilde{E}_l) / \Delta_t \\ -(\mu \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}} +\epsilon \tilde{E}^2_l) / \Delta_t \\ - \hat{H}_{l-\frac{1}{2}} \cdot \hat{M}_l \\ - \tilde{E}_l \cdot \tilde{J}_{l-\frac{1}{2}} \\ \end{aligned} \]

These two results form the discrete time-domain analogue to Poynting's theorem. They hint at the expressions for the energy, which can be calculated at the same time-index as either the E or H field:

\[ \begin{aligned} U_l &= \epsilon \tilde{E}^2_l + \mu \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}} \\ U_{l + \frac{1}{2}} &= \epsilon \tilde{E}_l \cdot \tilde{E}_{l + 1} + \mu \hat{H}^2_{l + \frac{1}{2}} \\ \end{aligned} \]

Rewriting the Poynting theorem in terms of the energy expressions,

\[ \begin{aligned} (U_{l+\frac{1}{2}} - U_l) / \Delta_t &= -\hat{\nabla} \cdot \tilde{S}_{l, l + \frac{1}{2}} \\ - \hat{H}_{l+\frac{1}{2}} \cdot \hat{M}_l \\ - \tilde{E}_l \cdot \tilde{J}_{l+\frac{1}{2}} \\ (U_l - U_{l-\frac{1}{2}}) / \Delta_t &= -\hat{\nabla} \cdot \tilde{S}_{l, l - \frac{1}{2}} \\ - \hat{H}_{l-\frac{1}{2}} \cdot \hat{M}_l \\ - \tilde{E}_l \cdot \tilde{J}_{l-\frac{1}{2}} \\ \end{aligned} \]

This result is exact and should practically hold to within numerical precision. No time- or spatial-averaging is necessary.

Note that each value of \(J\) contributes to the energy twice (i.e. once per field update) despite only causing the value of \(E\) to change once (same for \(M\) and \(H\)).

Sources¶

It is often useful to excite the simulation with an arbitrary broadband pulse and then extract the frequency-domain response by performing an on-the-fly Fourier transform of the time-domain fields.

accumulate_phasor in meanas.fdtd.phasor performs the phase accumulation for one or more target frequencies, while leaving source normalization and simulation-loop policy to the caller. Convenience wrappers accumulate_phasor_e, accumulate_phasor_h, and accumulate_phasor_j apply the usual Yee time offsets. The helpers accumulate

\[ \Delta_t \sum_l w_l e^{-i \omega t_l} f_l \]

with caller-provided sample times and weights. In this codebase the matching electric-current convention is typically E -= dt * J / epsilon in FDTD and -i \omega J on the right-hand side of the FDFD wave equation.

The Ricker wavelet (normalized second derivative of a Gaussian) is commonly used for the pulse shape. It can be written

\[ f_r(t) = (1 - \frac{1}{2} (\omega (t - \tau))^2) e^{-(\frac{\omega (t - \tau)}{2})^2} \]

with \(\tau > \frac{2 * \pi}{\omega}\) as a minimum delay to avoid a discontinuity at t=0 (assuming the source is off for t<0 this gives \(\sim 10^{-3}\) error at t=0).

Boundary conditions¶

meanas.fdtd exposes two boundary-related building blocks:

  • conducting_boundary(...) for simple perfect-electric-conductor style field clamping at one face of the domain.
  • cpml_params(...) and updates_with_cpml(...) for convolutional perfectly matched layers (CPMLs) attached to one or more faces of the Yee grid.

updates_with_cpml(...) accepts a three-by-two table of CPML parameter blocks:

cpml_params[axis][polarity_index]

where axis is 0, 1, or 2 and polarity_index corresponds to (-1, +1). Passing None for one entry disables CPML on that face while leaving the other directions unchanged. This is how mixed boundary setups such as "absorbing in x, periodic in y/z" are expressed.

When comparing an FDTD run against an FDFD solve, use the same stretched coordinate system in both places:

  1. Build the FDTD update with the desired CPML parameters.
  2. Stretch the FDFD dxes with the matching SCPML transform.
  3. Compare the extracted phasor against the FDFD field or residual on those stretched dxes.

The electric-current sign convention used throughout the examples and tests is

\[ E \leftarrow E - \Delta_t J / \epsilon \]

which matches the FDFD right-hand side

\[ -i \omega J. \]

Core update and analysis helpers¶

meanas.fdtd.base ¶

Basic FDTD field updates

maxwell_e ¶

maxwell_e(
    dt: float, dxes: dx_lists_t | None = None
) -> fdfield_updater_t

Build a function which performs a portion the time-domain E-field update,

E += curl_back(H[t]) / epsilon

The full update should be

E += (curl_back(H[t]) + J) / epsilon

which requires an additional step of E += J / epsilon which is not performed by the generated function.

See meanas.fdmath for descriptions of

  • This update step: "Maxwell's equations" section
  • dxes: "Datastructure: dx_lists_t" section
  • epsilon: "Permittivity and Permeability" section

Also see the "Timestep" section of meanas.fdtd for a discussion of the dt parameter.

Parameters:

Name Type Description Default
dt float

Timestep. See meanas.fdtd for details.

 required
dxes dx_lists_t | None

Grid description; see meanas.fdmath.

 None

Returns:

Type Description
fdfield_updater_t

Function f(E_old, H_old, epsilon) -> E_new.

maxwell_h ¶

maxwell_h(
    dt: float, dxes: dx_lists_t | None = None
) -> fdfield_updater_t

Build a function which performs part of the time-domain H-field update,

H -= curl_forward(E[t]) / mu

The full update should be

H -= (curl_forward(E[t]) + M) / mu

which requires an additional step of H -= M / mu which is not performed by the generated function; this step can be omitted if there is no magnetic current M.

See meanas.fdmath for descriptions of

  • This update step: "Maxwell's equations" section
  • dxes: "Datastructure: dx_lists_t" section
  • mu: "Permittivity and Permeability" section

Also see the "Timestep" section of meanas.fdtd for a discussion of the dt parameter.

Parameters:

Name Type Description Default
dt float

Timestep. See meanas.fdtd for details.

 required
dxes dx_lists_t | None

Grid description; see meanas.fdmath.

 None

Returns:

Type Description
fdfield_updater_t

Function f(E_old, H_old, epsilon) -> E_new.

meanas.fdtd.pml ¶

Convolutional perfectly matched layer (CPML) support for FDTD updates.

The helpers in this module construct per-face CPML parameters and then wrap the standard Yee updates with the additional auxiliary psi fields needed by the CPML recurrence.

The intended call pattern is:

  1. Build a cpml_params[axis][polarity_index] table with cpml_params(...).
  2. Pass that table into updates_with_cpml(...) together with dt, dxes, and epsilon.
  3. Advance the returned update_E / update_H closures in the simulation loop.

Each face can be enabled or disabled independently by replacing one table entry with None.

cpml_params ¶

cpml_params(
    axis: int,
    polarity: int,
    dt: float,
    thickness: int = 8,
    ln_R_per_layer: float = -1.6,
    epsilon_eff: float = 1,
    mu_eff: float = 1,
    m: float = 3.5,
    ma: float = 1,
    cfs_alpha: float = 0,
) -> dict[str, Any]

Construct the parameter block for one CPML face.

Parameters:

Name Type Description Default
axis int

Which Cartesian axis the CPML is normal to (0, 1, or 2).

 required
polarity int

Which face along that axis (-1 for the low-index face, +1 for the high-index face).

 required
dt float

Timestep used by the Yee update.

 required
thickness int

Number of Yee cells occupied by the CPML region.

 8
ln_R_per_layer float

Logarithmic attenuation target per layer.

 -1.6
epsilon_eff float

Effective permittivity used when choosing the CPML scaling.

 1
mu_eff float

Effective permeability used when choosing the CPML scaling.

 1
m float

Polynomial grading exponent for sigma and kappa.

 3.5
ma float

Polynomial grading exponent for the complex-frequency shift alpha.

 1
cfs_alpha float

Maximum complex-frequency shift parameter.

 0

Returns:

Type Description
dict[str, Any]

Dictionary with:
dict[str, Any]
  • param_e: (p0, p1, p2) arrays for the E update
dict[str, Any]
  • param_h: (p0, p1, p2) arrays for the H update
dict[str, Any]
  • region: slice tuple selecting the CPML cells on that face

updates_with_cpml ¶

updates_with_cpml(
    cpml_params: Sequence[Sequence[dict[str, Any] | None]],
    dt: float,
    dxes: dx_lists_t,
    epsilon: fdfield,
    *,
    dtype: DTypeLike = numpy.float32,
) -> tuple[
    Callable[[fdfield_t, fdfield_t, fdfield_t], None],
    Callable[[fdfield_t, fdfield_t, fdfield_t], None],
]

Build Yee-step update closures augmented with CPML terms.

Parameters:

Name Type Description Default
cpml_params Sequence[Sequence[dict[str, Any] | None]]

Three-by-two sequence indexed as [axis][polarity_index]. Entries are the dictionaries returned by cpml_params(...); use None to disable CPML on one face.

 required
dt float

Timestep.

 required
dxes dx_lists_t

Yee-grid spacing lists [dx_e, dx_h].

 required
epsilon fdfield

Electric material distribution used by the E update.

 required
dtype DTypeLike

Storage dtype for the auxiliary CPML state arrays.

 float32

Returns:

Type Description
Callable[[fdfield_t, fdfield_t, fdfield_t], None]

(update_E, update_H) closures with the same call shape as the basic
Callable[[fdfield_t, fdfield_t, fdfield_t], None]

Yee updates:
tuple[Callable[[fdfield_t, fdfield_t, fdfield_t], None], Callable[[fdfield_t, fdfield_t, fdfield_t], None]]
  • update_E(e, h, epsilon)
tuple[Callable[[fdfield_t, fdfield_t, fdfield_t], None], Callable[[fdfield_t, fdfield_t, fdfield_t], None]]
  • update_H(e, h, mu)
tuple[Callable[[fdfield_t, fdfield_t, fdfield_t], None], Callable[[fdfield_t, fdfield_t, fdfield_t], None]]

The closures retain the CPML auxiliary state internally.

meanas.fdtd.boundaries ¶

Boundary conditions

TODO conducting boundary documentation¶

meanas.fdtd.energy ¶

poynting ¶

poynting(
    e: fdfield, h: fdfield, dxes: dx_lists_t | None = None
) -> fdfield_t

Calculate the poynting vector S (\(S\)).

This is the energy transfer rate (amount of energy U per dt transferred between adjacent cells) in each direction that happens during the half-step bounded by the two provided fields.

The returned vector field S is the energy flow across +x, +y, and +z boundaries of the corresponding U cell. For example,

    mx = numpy.roll(mask, -1, axis=0)
    my = numpy.roll(mask, -1, axis=1)
    mz = numpy.roll(mask, -1, axis=2)

    u_hstep = fdtd.energy_hstep(e0=es[ii - 1], h1=hs[ii], e2=es[ii],     **args)
    u_estep = fdtd.energy_estep(h0=hs[ii],     e1=es[ii], h2=hs[ii + 1], **args)
    delta_j_B = fdtd.delta_energy_j(j0=js[ii], e1=es[ii], dxes=dxes)
    du_half_h2e = u_estep - u_hstep - delta_j_B

    s_h2e = -fdtd.poynting(e=es[ii], h=hs[ii], dxes=dxes) * dt
    planes = [s_h2e[0, mask].sum(), -s_h2e[0, mx].sum(),
              s_h2e[1, mask].sum(), -s_h2e[1, my].sum(),
              s_h2e[2, mask].sum(), -s_h2e[2, mz].sum()]

    assert_close(sum(planes), du_half_h2e[mask])

(see meanas.tests.test_fdtd.test_poynting_planes)

The full relationship is $$ \begin{aligned} (U_{l+\frac{1}{2}} - U_l) / \Delta_t &= -\hat{\nabla} \cdot \tilde{S}{l, l + \frac{1}{2}} \ - \hat{H}}{2}} \cdot \hat{Ml \ - \tilde{E}_l \cdot \tilde{J} \ (U_l - U_{l-\frac{1}{2}}) / \Delta_t &= -\hat{\nabla} \cdot \tilde{S}}{2}{l, l - \frac{1}{2}} \ - \hat{H}}{2}} \cdot \hat{Ml \ - \tilde{E}_l \cdot \tilde{J} \ \end{aligned} $$}{2}

These equalities are exact and should practically hold to within numerical precision. No time- or spatial-averaging is necessary. (See meanas.fdtd docs for derivation.)

Parameters:

Name Type Description Default
e fdfield

E-field

 required
h fdfield

H-field (one half-timestep before or after e)

 required
dxes dx_lists_t | None

Grid description; see meanas.fdmath.

 None

Returns:

Name Type Description
s fdfield_t

Vector field. Components indicate the energy transfer rate from the corresponding energy cell into its +x, +y, and +z neighbors during the half-step from the time of the earlier input field until the time of later input field.

poynting_divergence ¶

poynting_divergence(
    s: fdfield | None = None,
    *,
    e: fdfield | None = None,
    h: fdfield | None = None,
    dxes: dx_lists_t | None = None,
) -> fdfield_t

Calculate the divergence of the poynting vector.

This is the net energy flow for each cell, i.e. the change in energy U per dt caused by transfer of energy to nearby cells (rather than absorption/emission by currents J or M).

See poynting and meanas.fdtd for more details. Args: s: Poynting vector, as calculated with poynting. Optional; caller can provide e and h instead. e: E-field (optional; need either s or both e and h) h: H-field (optional; need either s or both e and h) dxes: Grid description; see meanas.fdmath.

Returns:

Name Type Description
ds fdfield_t

Divergence of the poynting vector. Entries indicate the net energy flow out of the corresponding energy cell.

energy_hstep ¶

energy_hstep(
    e0: fdfield,
    h1: fdfield,
    e2: fdfield,
    epsilon: fdfield | None = None,
    mu: fdfield | None = None,
    dxes: dx_lists_t | None = None,
) -> fdfield_t

Calculate energy U at the time of the provided H-field h1.

TODO: Figure out what this means spatially.

Parameters:

Name Type Description Default
e0 fdfield

E-field one half-timestep before the energy.

 required
h1 fdfield

H-field (at the same timestep as the energy).

 required
e2 fdfield

E-field one half-timestep after the energy.

 required
epsilon fdfield | None

Dielectric constant distribution.

 None
mu fdfield | None

Magnetic permeability distribution.

 None
dxes dx_lists_t | None

Grid description; see meanas.fdmath.

 None

Returns:

Type Description
fdfield_t

Energy, at the time of the H-field h1.

energy_estep ¶

energy_estep(
    h0: fdfield,
    e1: fdfield,
    h2: fdfield,
    epsilon: fdfield | None = None,
    mu: fdfield | None = None,
    dxes: dx_lists_t | None = None,
) -> fdfield_t

Calculate energy U at the time of the provided E-field e1.

TODO: Figure out what this means spatially.

Parameters:

Name Type Description Default
h0 fdfield

H-field one half-timestep before the energy.

 required
e1 fdfield

E-field (at the same timestep as the energy).

 required
h2 fdfield

H-field one half-timestep after the energy.

 required
epsilon fdfield | None

Dielectric constant distribution.

 None
mu fdfield | None

Magnetic permeability distribution.

 None
dxes dx_lists_t | None

Grid description; see meanas.fdmath.

 None

Returns:

Type Description
fdfield_t

Energy, at the time of the E-field e1.

delta_energy_h2e ¶

delta_energy_h2e(
    dt: float,
    e0: fdfield,
    h1: fdfield,
    e2: fdfield,
    h3: fdfield,
    epsilon: fdfield | None = None,
    mu: fdfield | None = None,
    dxes: dx_lists_t | None = None,
) -> fdfield_t

Change in energy during the half-step from h1 to e2.

This is just from (e2 * e2 + h3 * h1) - (h1 * h1 + e0 * e2)

Parameters:

Name Type Description Default
e0 fdfield

E-field one half-timestep before the start of the energy delta.

 required
h1 fdfield

H-field at the start of the energy delta.

 required
e2 fdfield

E-field at the end of the energy delta (one half-timestep after h1).

 required
h3 fdfield

H-field one half-timestep after the end of the energy delta.

 required
epsilon fdfield | None

Dielectric constant distribution.

 None
mu fdfield | None

Magnetic permeability distribution.

 None
dxes dx_lists_t | None

Grid description; see meanas.fdmath.

 None

Returns:

Type Description
fdfield_t

Change in energy from the time of h1 to the time of e2.

delta_energy_e2h ¶

delta_energy_e2h(
    dt: float,
    h0: fdfield,
    e1: fdfield,
    h2: fdfield,
    e3: fdfield,
    epsilon: fdfield | None = None,
    mu: fdfield | None = None,
    dxes: dx_lists_t | None = None,
) -> fdfield_t

Change in energy during the half-step from e1 to h2.

This is just from (h2 * h2 + e3 * e1) - (e1 * e1 + h0 * h2)

Parameters:

Name Type Description Default
h0 fdfield

E-field one half-timestep before the start of the energy delta.

 required
e1 fdfield

H-field at the start of the energy delta.

 required
h2 fdfield

E-field at the end of the energy delta (one half-timestep after e1).

 required
e3 fdfield

H-field one half-timestep after the end of the energy delta.

 required
epsilon fdfield | None

Dielectric constant distribution.

 None
mu fdfield | None

Magnetic permeability distribution.

 None
dxes dx_lists_t | None

Grid description; see meanas.fdmath.

 None

Returns:

Type Description
fdfield_t

Change in energy from the time of e1 to the time of h2.

delta_energy_j ¶

delta_energy_j(
    j0: fdfield, e1: fdfield, dxes: dx_lists_t | None = None
) -> fdfield_t

Calculate the electric-current work term \(J \cdot E\) on the Yee grid.

This is the source contribution that appears beside the flux divergence in the discrete Poynting identities documented in meanas.fdtd.

Note that each value of J contributes twice in a full Yee cycle (once per half-step energy balance) even though it directly changes E only once.

Parameters:

Name Type Description Default
j0 fdfield

Electric-current density sampled at the same half-step as the current work term.

 required
e1 fdfield

Electric field sampled at the matching integer timestep.

 required
dxes dx_lists_t | None

Grid description; defaults to unit spacing.

 None

Returns:

Type Description
fdfield_t

Per-cell source-work contribution with the scalar field shape.

dxmul ¶

dxmul(
    ee: fdfield,
    hh: fdfield,
    epsilon: fdfield | float | None = None,
    mu: fdfield | float | None = None,
    dxes: dx_lists_t | None = None,
) -> fdfield_t

Multiply E- and H-like field products by material weights and cell volumes.

Parameters:

Name Type Description Default
ee fdfield

Three-component electric-field product, such as e0 * e2.

 required
hh fdfield

Three-component magnetic-field product, such as h1 * h1.

 required
epsilon fdfield | float | None

Electric material weight; defaults to 1.

 None
mu fdfield | float | None

Magnetic material weight; defaults to 1.

 None
dxes dx_lists_t | None

Grid description; defaults to unit spacing.

 None

Returns:

Type Description
fdfield_t

Scalar field containing the weighted electric plus magnetic contribution
fdfield_t

for each Yee cell.

meanas.fdtd.phasor ¶

Helpers for extracting single- or multi-frequency phasors from FDTD samples.

These helpers are intentionally low-level: callers own the accumulator arrays and the sampling policy. The accumulated quantity is

dt * sum(weight * exp(-1j * omega * t_step) * sample_step)

where t_step = (step + offset_steps) * dt.

The usual Yee offsets are:

  • accumulate_phasor_e(..., step=l) for E_l
  • accumulate_phasor_h(..., step=l) for H_{l + 1/2}
  • accumulate_phasor_j(..., step=l) for J_{l + 1/2}

These helpers do not choose warmup/accumulation windows or pulse-envelope normalization. They also do not impose a current sign convention. In this codebase, electric-current injection normally appears as E -= dt * J / epsilon, which matches the FDFD right-hand side -1j * omega * J.

accumulate_phasor ¶

accumulate_phasor(
    accumulator: NDArray[complexfloating],
    omegas: float
    | complex
    | Sequence[float | complex]
    | NDArray,
    dt: float,
    sample: ArrayLike,
    step: int,
    *,
    offset_steps: float = 0.0,
    weight: ArrayLike = 1.0,
) -> NDArray[numpy.complexfloating]

Add one time-domain sample into a phasor accumulator.

The added quantity is

dt * weight * exp(-1j * omega * t_step) * sample

where t_step = (step + offset_steps) * dt.

Note

This helper already multiplies by dt. If the caller's normalization factor was derived from a discrete sum that already includes dt, pass weight / dt here.

accumulate_phasor_e ¶

accumulate_phasor_e(
    accumulator: NDArray[complexfloating],
    omegas: float
    | complex
    | Sequence[float | complex]
    | NDArray,
    dt: float,
    sample: ArrayLike,
    step: int,
    *,
    weight: ArrayLike = 1.0,
) -> NDArray[numpy.complexfloating]

Accumulate an E-field sample taken at integer timestep step.

accumulate_phasor_h ¶

accumulate_phasor_h(
    accumulator: NDArray[complexfloating],
    omegas: float
    | complex
    | Sequence[float | complex]
    | NDArray,
    dt: float,
    sample: ArrayLike,
    step: int,
    *,
    weight: ArrayLike = 1.0,
) -> NDArray[numpy.complexfloating]

Accumulate an H-field sample corresponding to H_{step + 1/2}.

accumulate_phasor_j ¶

accumulate_phasor_j(
    accumulator: NDArray[complexfloating],
    omegas: float
    | complex
    | Sequence[float | complex]
    | NDArray,
    dt: float,
    sample: ArrayLike,
    step: int,
    *,
    weight: ArrayLike = 1.0,
) -> NDArray[numpy.complexfloating]

Accumulate a current sample corresponding to J_{step + 1/2}.

fdmath¶

meanas.fdmath ¶

Basic discrete calculus for finite difference (fd) simulations.

Fields, Functions, and Operators¶

Discrete fields are stored in one of two forms:

  • The fdfield_t form is a multidimensional numpy.NDArray
    • For a scalar field, this is just U[m, n, p], where m, n, and p are discrete indices referring to positions on the x, y, and z axes respectively.
    • For a vector field, the first index specifies which vector component is accessed: E[:, m, n, p] = [Ex[m, n, p], Ey[m, n, p], Ez[m, n, p]].
  • The vfdfield_t form is simply a vectorzied (i.e. 1D) version of the fdfield_t, as obtained by meanas.fdmath.vectorization.vec (effectively just numpy.ravel)
Operators which act on fields also come in two forms
  • Python functions, created by the functions in meanas.fdmath.functional. The generated functions act on fields in the fdfield_t form.
  • Linear operators, usually 2D sparse matrices using scipy.sparse, created by meanas.fdmath.operators. These operators act on vectorized fields in the vfdfield_t form.

The operations performed should be equivalent: functional.op(*args)(E) should be equivalent to unvec(operators.op(*args) @ vec(E), E.shape[1:]).

Generally speaking the field_t form is easier to work with, but can be harder or less efficient to compose (e.g. it is easy to generate a single matrix by multiplying a series of other matrices).

Discrete calculus¶

This documentation and approach is roughly based on W.C. Chew's excellent "Electromagnetic Theory on a Lattice" (doi:10.1063/1.355770), which covers a superset of this material with similar notation and more detail.

Scalar derivatives and cell shifts¶

Define the discrete forward derivative as $$ [\tilde{\partial}x f] - f_m) $$ where }{2}} = \frac{1}{\Delta_{x, m}} (f_{m + 1\(f\) is a function defined at discrete locations on the x-axis (labeled using \(m\)). The value at \(m\) occupies a length \(\Delta_{x, m}\) along the x-axis. Note that \(m\) is an index along the x-axis, not necessarily an x-coordinate, since each length \(\Delta_{x, m}, \Delta_{x, m+1}, ...\) is independently chosen.

If we treat f as a 1D array of values, with the i-th value f[i] taking up a length dx[i] along the x-axis, the forward derivative is

deriv_forward(f)[i] = (f[i + 1] - f[i]) / dx[i]

Likewise, discrete reverse derivative is $$ [\hat{\partial}x f ]) $$ or}{2}} = \frac{1}{\Delta_{x, m}} (f_{m} - f_{m - 1

deriv_back(f)[i] = (f[i] - f[i - 1]) / dx[i]

The derivatives' values are shifted by a half-cell relative to the original function, and will have different cell widths if all the dx[i] ( \(\Delta_{x, m}\) ) are not identical:

[figure: derivatives and cell sizes]
    dx0   dx1      dx2      dx3      cell sizes for function
   ----- ----- ----------- -----
   ______________________________
        |     |           |     |
     f0 |  f1 |     f2    |  f3 |    function
   _____|_____|___________|_____|
     |     |        |        |
     | Df0 |   Df1  |   Df2  | Df3   forward derivative (periodic boundary)
   __|_____|________|________|___

 dx'3] dx'0   dx'1     dx'2  [dx'3   cell sizes for forward derivative
   -- ----- -------- -------- ---
 dx'0] dx'1   dx'2     dx'3  [dx'0   cell sizes for reverse derivative
   ______________________________
     |     |        |        |
     | df1 |  df2   |   df3  | df0   reverse derivative (periodic boundary)
   __|_____|________|________|___

Periodic boundaries are used here and elsewhere unless otherwise noted.

In the above figure, f0 = \(f_0\), f1 = \(f_1\) Df0 = \([\tilde{\partial}f]_{0 + \frac{1}{2}}\) Df1 = \([\tilde{\partial}f]_{1 + \frac{1}{2}}\) df0 = \([\hat{\partial}f]_{0 - \frac{1}{2}}\) etc.

The fractional subscript \(m + \frac{1}{2}\) is used to indicate values defined at shifted locations relative to the original \(m\), with corresponding lengths $$ \Delta_{x, m + \frac{1}{2}} = \frac{1}{2} * (\Delta_{x, m} + \Delta_{x, m + 1}) $$

Just as \(m\) is not itself an x-coordinate, neither is \(m + \frac{1}{2}\); carefully note the positions of the various cells in the above figure vs their labels. If the positions labeled with \(m\) are considered the "base" or "original" grid, the positions labeled with \(m + \frac{1}{2}\) are said to lie on a "dual" or "derived" grid.

For the remainder of the Discrete calculus section, all figures will show constant-length cells in order to focus on the vector derivatives themselves. See the Grid description section below for additional information on this topic and generalization to three dimensions.

Gradients and fore-vectors¶

Expanding to three dimensions, we can define two gradients $$ [\tilde{\nabla} f]{m,n,p} = \vec{x} [\tilde{\partial}_x f] + \vec{y} [\tilde{\partial}}{2},n,py f] + \vec{z} [\tilde{\partial}}{2},pz f] $$ $$ [\hat{\nabla} f]}{2}{m,n,p} = \vec{x} [\hat{\partial}_x f] + \vec{y} [\hat{\partial}}{2},n,py f] + \vec{z} [\hat{\partial}}{2},pz f] $$}{2}

or

[code: gradients]
grad_forward(f)[i,j,k] = [Dx_forward(f)[i, j, k],
                          Dy_forward(f)[i, j, k],
                          Dz_forward(f)[i, j, k]]
                       = [(f[i + 1, j, k] - f[i, j, k]) / dx[i],
                          (f[i, j + 1, k] - f[i, j, k]) / dy[i],
                          (f[i, j, k + 1] - f[i, j, k]) / dz[i]]

grad_back(f)[i,j,k] = [Dx_back(f)[i, j, k],
                       Dy_back(f)[i, j, k],
                       Dz_back(f)[i, j, k]]
                    = [(f[i, j, k] - f[i - 1, j, k]) / dx[i],
                       (f[i, j, k] - f[i, j - 1, k]) / dy[i],
                       (f[i, j, k] - f[i, j, k - 1]) / dz[i]]

The three derivatives in the gradient cause shifts in different directions, so the x/y/z components of the resulting "vector" are defined at different points: the x-component is shifted in the x-direction, y in y, and z in z.

We call the resulting object a "fore-vector" or "back-vector", depending on the direction of the shift. We write it as $$ \tilde{g}{m,n,p} = \vec{x} g^x + \vec{y} g^y_{m,n + \frac{1}{2},p} + \vec{z} g^z_{m,n,p + \frac{1}{2}} $$ $$ \hat{g}}{2},n,p{m,n,p} = \vec{x} g^x + \vec{y} g^y_{m,n - \frac{1}{2},p} + \vec{z} g^z_{m,n,p - \frac{1}{2}} $$}{2},n,p

[figure: gradient / fore-vector]
   (m, n+1, p+1) ______________ (m+1, n+1, p+1)
                /:            /|
               / :           / |
              /  :          /  |
  (m, n, p+1)/_____________/   |     The forward derivatives are defined
             |   :         |   |     at the Dx, Dy, Dz points,
             |   :.........|...|     but the forward-gradient fore-vector
 z y        Dz  /          |  /      is the set of all three
 |/_x        | Dy          | /       and is said to be "located" at (m,n,p)
             |/            |/
    (m, n, p)|_____Dx______| (m+1, n, p)

Divergences¶

There are also two divergences,

$$ d_{n,m,p} = [\tilde{\nabla} \cdot \hat{g}]{n,m,p} = [\tilde{\partial}_x g^x] + [\tilde{\partial}y g^y] + [\tilde{\partial}z g^z] $$

$$ d_{n,m,p} = [\hat{\nabla} \cdot \tilde{g}]{n,m,p} = [\hat{\partial}_x g^x] + [\hat{\partial}y g^y] + [\hat{\partial}z g^z] $$

or

[code: divergences]
div_forward(g)[i,j,k] = Dx_forward(gx)[i, j, k] +
                        Dy_forward(gy)[i, j, k] +
                        Dz_forward(gz)[i, j, k]
                      = (gx[i + 1, j, k] - gx[i, j, k]) / dx[i] +
                        (gy[i, j + 1, k] - gy[i, j, k]) / dy[i] +
                        (gz[i, j, k + 1] - gz[i, j, k]) / dz[i]

div_back(g)[i,j,k] = Dx_back(gx)[i, j, k] +
                     Dy_back(gy)[i, j, k] +
                     Dz_back(gz)[i, j, k]
                   = (gx[i, j, k] - gx[i - 1, j, k]) / dx[i] +
                     (gy[i, j, k] - gy[i, j - 1, k]) / dy[i] +
                     (gz[i, j, k] - gz[i, j, k - 1]) / dz[i]

where g = [gx, gy, gz] is a fore- or back-vector field.

Since we applied the forward divergence to the back-vector (and vice-versa), the resulting scalar value is defined at the back-vector's (fore-vector's) location \((m,n,p)\) and not at the locations of its components \((m \pm \frac{1}{2},n,p)\) etc.

[figure: divergence]
                                ^^
     (m-1/2, n+1/2, p+1/2) _____||_______ (m+1/2, n+1/2, p+1/2)
                          /:    ||  ,,  /|
                         / :    || //  / |      The divergence at (m, n, p) (the center
                        /  :      //  /  |      of this cube) of a fore-vector field
  (m-1/2, n-1/2, p+1/2)/_____________/   |      is the sum of the outward-pointing
                       |   :         |   |      fore-vector components, which are
     z y            <==|== :.........|.====>    located at the face centers.
     |/_x              |  /          |  /
                       | /    //     | /       Note that in a nonuniform grid, each
                       |/    // ||   |/        dimension is normalized by the cell width.
  (m-1/2, n-1/2, p-1/2)|____//_______| (m+1/2, n-1/2, p-1/2)
                           ''   ||
                                VV

Curls¶

The two curls are then

$$ \begin{aligned} \hat{h}{m + \frac{1}{2}, n + \frac{1}{2}, p + \frac{1}{2}} &= \ [\tilde{\nabla} \times \tilde{g}] &= \vec{x} (\tilde{\partial}}{2}, n + \frac{1}{2}, p + \frac{1}{2}y g^z}{2}} - \tilde{\partialz g^y) \ &+ \vec{y} (\tilde{\partial}}{2},pz g^x}{2},n,p} - \tilde{\partialx g^z) \ &+ \vec{z} (\tilde{\partial}}{2}x g^y}{2},p} - \tilde{\partialy g^z) \end{aligned} $$}{2},n,p

and

$$ \tilde{h}{m - \frac{1}{2}, n - \frac{1}{2}, p - \frac{1}{2}} = [\hat{\nabla} \times \hat{g}] $$}{2}, n - \frac{1}{2}, p - \frac{1}{2}

where \(\hat{g}\) and \(\tilde{g}\) are located at \((m,n,p)\) with components at \((m \pm \frac{1}{2},n,p)\) etc., while \(\hat{h}\) and \(\tilde{h}\) are located at \((m \pm \frac{1}{2}, n \pm \frac{1}{2}, p \pm \frac{1}{2})\) with components at \((m, n \pm \frac{1}{2}, p \pm \frac{1}{2})\) etc.

[code: curls]
curl_forward(g)[i,j,k] = [Dy_forward(gz)[i, j, k] - Dz_forward(gy)[i, j, k],
                          Dz_forward(gx)[i, j, k] - Dx_forward(gz)[i, j, k],
                          Dx_forward(gy)[i, j, k] - Dy_forward(gx)[i, j, k]]

curl_back(g)[i,j,k] = [Dy_back(gz)[i, j, k] - Dz_back(gy)[i, j, k],
                       Dz_back(gx)[i, j, k] - Dx_back(gz)[i, j, k],
                       Dx_back(gy)[i, j, k] - Dy_back(gx)[i, j, k]]

For example, consider the forward curl, at (m, n, p), of a back-vector field g, defined on a grid containing (m + 1/2, n + 1/2, p + 1/2). The curl will be a fore-vector, so its z-component will be defined at (m, n, p + 1/2). Take the nearest x- and y-components of g in the xy plane where the curl's z-component is located; these are

[curl components]
(m,       n + 1/2, p + 1/2) : x-component of back-vector at (m + 1/2, n + 1/2, p + 1/2)
(m + 1,   n + 1/2, p + 1/2) : x-component of back-vector at (m + 3/2, n + 1/2, p + 1/2)
(m + 1/2, n      , p + 1/2) : y-component of back-vector at (m + 1/2, n + 1/2, p + 1/2)
(m + 1/2, n + 1  , p + 1/2) : y-component of back-vector at (m + 1/2, n + 3/2, p + 1/2)

These four xy-components can be used to form a loop around the curl's z-component; its magnitude and sign is set by their loop-oriented sum (i.e. two have their signs flipped to complete the loop).

[figure: z-component of curl]
                          :             |
    z y                   :    ^^       |
    |/_x                  :....||.<.....|  (m+1, n+1, p+1/2)
                          /    ||      /
                       | v     ||   | ^
                       |/           |/
         (m, n, p+1/2) |_____>______|  (m+1, n, p+1/2)

Maxwell's Equations¶

If we discretize both space (m,n,p) and time (l), Maxwell's equations become

$$ \begin{aligned} \tilde{\nabla} \times \tilde{E}{l,\vec{r}} &= -\tilde{\partial}_t \hat{B} - \hat{M}}{2}, \vec{r} + \frac{1}{2}{l, \vec{r} + \frac{1}{2}} \ \hat{\nabla} \times \hat{H}}{2},\vec{r} + \frac{1}{2}} &= \hat{\partialt \tilde{D} + \tilde{J}}{l-\frac{1}{2},\vec{r}} \ \tilde{\nabla} \cdot \hat{B} &= 0 \ \hat{\nabla} \cdot \tilde{D}}{2}, \vec{r} + \frac{1}{2}{l,\vec{r}} &= \rho \end{aligned} $$}

with

$$ \begin{aligned} \hat{B}{\vec{r}} &= \mu} + \frac{1}{2}} \cdot \hat{H{\vec{r} + \frac{1}{2}} \ \tilde{D} \end{aligned} $$}} &= \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}

where the spatial subscripts are abbreviated as \(\vec{r} = (m, n, p)\) and \(\vec{r} + \frac{1}{2} = (m + \frac{1}{2}, n + \frac{1}{2}, p + \frac{1}{2})\), \(\tilde{E}\) and \(\hat{H}\) are the electric and magnetic fields, \(\tilde{J}\) and \(\hat{M}\) are the electric and magnetic current distributions, and \(\epsilon\) and \(\mu\) are the dielectric permittivity and magnetic permeability.

The above is Yee's algorithm, written in a form analogous to Maxwell's equations. The time derivatives can be expanded to form the update equations:

[code: Maxwell's equations updates]
H[i, j, k] -= dt * (curl_forward(E)[i, j, k] + M[t, i, j, k]) /      mu[i, j, k]
E[i, j, k] += dt * (curl_back(   H)[i, j, k] + J[t, i, j, k]) / epsilon[i, j, k]

Note that the E-field fore-vector and H-field back-vector are offset by a half-cell, resulting in distinct locations for all six E- and H-field components:

[figure: Field components]

        (m - 1/2,=> ____________Hx__________[H] <= r + 1/2 = (m + 1/2,
         n + 1/2,  /:           /:          /|                n + 1/2,
   z y   p + 1/2) / :          / :         / |                p + 1/2)
   |/_x          /  :         /  :        /  |
                /   :       Ez__________Hy   |      Locations of the E- and
               /    :        :   :      /|   |      H-field components for the
 (m - 1/2,    /     :        :  Ey...../.|..Hz      [E] fore-vector at r = (m,n,p)
  n - 1/2, =>/________________________/  |  /|      (the large cube's center)
  p + 1/2)   |      :        : /      |  | / |      and [H] back-vector at r + 1/2
             |      :        :/       |  |/  |      (the top right corner)
             |      :       [E].......|.Ex   |
             |      :.................|......| <= (m + 1/2, n + 1/2, p + 1/2)
             |     /                  |     /
             |    /                   |    /
             |   /                    |   /         This is the Yee discretization
             |  /                     |  /          scheme ("Yee cell").
r - 1/2 =    | /                      | /
 (m - 1/2,   |/                       |/
  n - 1/2,=> |________________________| <= (m + 1/2, n - 1/2, p - 1/2)
  p - 1/2)

Each component forms its own grid, offset from the others:

[figure: E-fields for adjacent cells]

              H1__________Hx0_________H0
  z y        /:                       /|
  |/_x      / :                      / |    This figure shows H back-vector locations
           /  :                     /  |    H0, H1, etc. and their associated components
         Hy1  :                   Hy0  |    H0 = (Hx0, Hy0, Hz0) etc.
         /    :                   /    |
        /    Hz1                 /     Hz0
       H2___________Hx3_________H3     |    The equivalent drawing for E would have
       |      :                 |      |    fore-vectors located at the cube's
       |      :                 |      |    center (and the centers of adjacent cubes),
       |      :                 |      |    with components on the cube's faces.
       |      H5..........Hx4...|......H4
       |     /                  |     /
      Hz2   /                  Hz2   /
       |   /                    |   /
       | Hy6                    | Hy4
       | /                      | /
       |/                       |/
       H6__________Hx7__________H7

The divergence equations can be derived by taking the divergence of the curl equations and combining them with charge continuity, $$ \hat{\nabla} \cdot \tilde{J} + \hat{\partial}_t \rho = 0 $$ implying that the discrete Maxwell's equations do not produce spurious charges.

Wave equation¶

Taking the backward curl of the \(\tilde{\nabla} \times \tilde{E}\) equation and replacing the resulting \(\hat{\nabla} \times \hat{H}\) term using its respective equation, and setting \(\hat{M}\) to zero, we can form the discrete wave equation:

\[ \begin{aligned} \tilde{\nabla} \times \tilde{E}_{l,\vec{r}} &= -\tilde{\partial}_t \hat{B}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}} - \hat{M}_{l-1, \vec{r} + \frac{1}{2}} \\ \mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}} &= -\tilde{\partial}_t \hat{H}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}} \\ \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}}) &= \hat{\nabla} \times (-\tilde{\partial}_t \hat{H}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}}) \\ \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}}) &= -\tilde{\partial}_t \hat{\nabla} \times \hat{H}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}} \\ \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}}) &= -\tilde{\partial}_t \hat{\partial}_t \epsilon_{\vec{r}} \tilde{E}_{l, \vec{r}} + \hat{\partial}_t \tilde{J}_{l-\frac{1}{2},\vec{r}} \\ \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}}) + \tilde{\partial}_t \hat{\partial}_t \epsilon_{\vec{r}} \cdot \tilde{E}_{l, \vec{r}} &= \tilde{\partial}_t \tilde{J}_{l - \frac{1}{2}, \vec{r}} \end{aligned} \]

Frequency domain¶

We can substitute in a time-harmonic fields

\[ \begin{aligned} \tilde{E}_{l, \vec{r}} &= \tilde{E}_{\vec{r}} e^{-\imath \omega l \Delta_t} \\ \tilde{J}_{l, \vec{r}} &= \tilde{J}_{\vec{r}} e^{-\imath \omega (l - \frac{1}{2}) \Delta_t} \end{aligned} \]

resulting in

\[ \begin{aligned} \tilde{\partial}_t &\Rightarrow (e^{ \imath \omega \Delta_t} - 1) / \Delta_t = \frac{-2 \imath}{\Delta_t} \sin(\omega \Delta_t / 2) e^{-\imath \omega \Delta_t / 2} = -\imath \Omega e^{-\imath \omega \Delta_t / 2}\\ \hat{\partial}_t &\Rightarrow (1 - e^{-\imath \omega \Delta_t}) / \Delta_t = \frac{-2 \imath}{\Delta_t} \sin(\omega \Delta_t / 2) e^{ \imath \omega \Delta_t / 2} = -\imath \Omega e^{ \imath \omega \Delta_t / 2}\\ \Omega &= 2 \sin(\omega \Delta_t / 2) / \Delta_t \end{aligned} \]

This gives the frequency-domain wave equation,

\[ \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{\vec{r}}) -\Omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} = -\imath \Omega \tilde{J}_{\vec{r}} e^{\imath \omega \Delta_t / 2} \\ \]

Plane waves and Dispersion relation¶

With uniform material distribution and no sources

\[ \begin{aligned} \mu_{\vec{r} + \frac{1}{2}} &= \mu \\ \epsilon_{\vec{r}} &= \epsilon \\ \tilde{J}_{\vec{r}} &= 0 \\ \end{aligned} \]

the frequency domain wave equation simplifies to

\[ \hat{\nabla} \times \tilde{\nabla} \times \tilde{E}_{\vec{r}} - \Omega^2 \epsilon \mu \tilde{E}_{\vec{r}} = 0 \]

Since \(\hat{\nabla} \cdot \tilde{E}_{\vec{r}} = 0\), we can simplify

\[ \begin{aligned} \hat{\nabla} \times \tilde{\nabla} \times \tilde{E}_{\vec{r}} &= \tilde{\nabla}(\hat{\nabla} \cdot \tilde{E}_{\vec{r}}) - \hat{\nabla} \cdot \tilde{\nabla} \tilde{E}_{\vec{r}} \\ &= - \hat{\nabla} \cdot \tilde{\nabla} \tilde{E}_{\vec{r}} \\ &= - \tilde{\nabla}^2 \tilde{E}_{\vec{r}} \end{aligned} \]

and we get

\[ \tilde{\nabla}^2 \tilde{E}_{\vec{r}} + \Omega^2 \epsilon \mu \tilde{E}_{\vec{r}} = 0 \]

We can convert this to three scalar-wave equations of the form

\[ (\tilde{\nabla}^2 + K^2) \phi_{\vec{r}} = 0 \]

with \(K^2 = \Omega^2 \mu \epsilon\). Now we let

\[ \phi_{\vec{r}} = A e^{\imath (k_x m \Delta_x + k_y n \Delta_y + k_z p \Delta_z)} \]

resulting in

\[ \begin{aligned} \tilde{\partial}_x &\Rightarrow (e^{ \imath k_x \Delta_x} - 1) / \Delta_t = \frac{-2 \imath}{\Delta_x} \sin(k_x \Delta_x / 2) e^{ \imath k_x \Delta_x / 2} = \imath K_x e^{ \imath k_x \Delta_x / 2}\\ \hat{\partial}_x &\Rightarrow (1 - e^{-\imath k_x \Delta_x}) / \Delta_t = \frac{-2 \imath}{\Delta_x} \sin(k_x \Delta_x / 2) e^{-\imath k_x \Delta_x / 2} = \imath K_x e^{-\imath k_x \Delta_x / 2}\\ K_x &= 2 \sin(k_x \Delta_x / 2) / \Delta_x \\ \end{aligned} \]

with similar expressions for the y and z dimnsions (and \(K_y, K_z\)).

This implies

\[ \tilde{\nabla}^2 = -(K_x^2 + K_y^2 + K_z^2) \phi_{\vec{r}} \\ K_x^2 + K_y^2 + K_z^2 = \Omega^2 \mu \epsilon = \Omega^2 / c^2 \]

where \(c = \sqrt{\mu \epsilon}\).

Assuming real \((k_x, k_y, k_z), \omega\) will be real only if

\[ c^2 \Delta_t^2 = \frac{\Delta_t^2}{\mu \epsilon} < 1/(\frac{1}{\Delta_x^2} + \frac{1}{\Delta_y^2} + \frac{1}{\Delta_z^2}) \]

If \(\Delta_x = \Delta_y = \Delta_z\), this simplifies to \(c \Delta_t < \Delta_x / \sqrt{3}\). This last form can be interpreted as enforcing causality; the distance that light travels in one timestep (i.e., \(c \Delta_t\)) must be less than the diagonal of the smallest cell ( \(\Delta_x / \sqrt{3}\) when on a uniform cubic grid).

Grid description¶

As described in the section on scalar discrete derivatives above, cell widths (dx[i], dy[j], dz[k]) along each axis can be arbitrary and independently defined. Moreover, all field components are actually defined at "derived" or "dual" positions, in-between the "base" grid points on one or more axes.

To get a better sense of how this works, let's start by drawing a grid with uniform dy and dz and nonuniform dx. We will only draw one cell in the y and z dimensions to make the illustration simpler; we need at least two cells in the x dimension to demonstrate how nonuniform dx affects the various components.

Place the E fore-vectors at integer indices \(r = (m, n, p)\) and the H back-vectors at fractional indices \(r + \frac{1}{2} = (m + \frac{1}{2}, n + \frac{1}{2}, p + \frac{1}{2})\). Remember that these are indices and not coordinates; they can correspond to arbitrary (monotonically increasing) coordinates depending on the cell widths.

Draw lines to denote the planes on which the H components and back-vectors are defined. For simplicity, don't draw the equivalent planes for the E components and fore-vectors, except as necessary to show their locations -- it's easiest to just connect them to their associated H-equivalents.

The result looks something like this:

[figure: Component centers]
                                                             p=
          [H]__________Hx___________[H]_____Hx______[H]   __ +1/2
  z y     /:           /:           /:      /:      /|     |      |
  |/_x   / :          / :          / :     / :     / |     |      |
        /  :         /  :         /  :    /  :    /  |     |      |
      Hy   :       Ez...........Hy   :  Ez......Hy   |     |      |
      /:   :        :   :       /:   :   :   :  /|   |     |      |
     / :  Hz        :  Ey....../.:..Hz   :  Ey./.|..Hz    __ 0    | dz[0]
    /  :  /:        :  /      /  :  /:   :  / /  |  /|     |      |
   /_________________________/_______________/   | / |     |      |
   |   :/  :        :/       |   :/  :   :/  |   |/  |     |      |
   |  Ex   :       [E].......|..Ex   :  [E]..|..Ex   |     |      |
   |       :                 |       :       |       |     |      |
   |      [H]..........Hx....|......[H].....H|x.....[H]   __ --------- (n=+1/2, p=-1/2)
   |      /                  |      /        |      /     /       /
  Hz     /                  Hz     /        Hz     /     /       /
   |    /                    |    /          |    /     /       /
   |  Hy                     |  Hy           |  Hy    __ 0     / dy[0]
   |  /                      |  /            |  /     /       /
   | /                       | /             | /     /       /
   |/                        |/              |/     /       /
  [H]__________Hx___________[H]_____Hx______[H]   __ -1/2  /
                                                       =n
   |------------|------------|-------|-------|
 -1/2           0          +1/2     +1     +3/2 = m

    ------------------------- ----------------
              dx[0]                  dx[1]

  Part of a nonuniform "base grid", with labels specifying
  positions of the various field components. [E] fore-vectors
  are at the cell centers, and [H] back-vectors are at the
  vertices. H components along the near (-y) top (+z) edge
  have been omitted to make the insides of the cubes easier
  to visualize.

The above figure shows where all the components are located; however, it is also useful to show what volumes those components correspond to. Consider the Ex component at m = +1/2: it is shifted in the x-direction by a half-cell from the E fore-vector at m = 0 (labeled [E] in the figure). It corresponds to a volume between m = 0 and m = +1 (the other dimensions are not shifted, i.e. they are still bounded by n, p = +-1/2). (See figure below). Since m is an index and not an x-coordinate, the Ex component is not necessarily at the center of the volume it represents, and the x-length of its volume is the derived quantity dx'[0] = (dx[0] + dx[1]) / 2 rather than the base dx. (See also Scalar derivatives and cell shifts).

[figure: Ex volumes]
                                                             p=
           <_________________________________________>   __ +1/2
  z y     <<           /:           /       /:      >>    |      |
  |/_x   < <          / :          /       / :     > >    |      |
        <  <         /  :         /       /  :    >  >    |      |
       <   <        /   :        /       /   :   >   >    |      |
      <:   <       /    :        :      /    :  >:   >    |      |
     < :   <      /     :        :     /     : > :   >   __ 0    | dz[0]
    <  :   <     /      :        :    /      :>  :   >    |      |
   <____________/____________________/_______>   :   >    |      |
   <   :   <    |       :        :   |       >   :   >    |      |
   <  Ex   <    |       :       Ex   |       >  Ex   >    |      |
   <   :   <    |       :        :   |       >   :   >    |      |
   <   :   <....|.......:........:...|.......>...:...>   __ --------- (n=+1/2, p=-1/2)
   <   :  <     |      /         :  /|      />   :  >    /       /
   <   : <      |     /          : / |     / >   : >    /       /
   <   :<       |    /           :/  |    /  >   :>    /       /
   <   <        |   /            :   |   /   >   >    _ 0     / dy[0]
   <  <         |  /                 |  /    >  >    /       /
   < <          | /                  | /     > >    /       /
   <<           |/                   |/      >>    /       /
   <____________|____________________|_______>   __ -1/2  /
                                                     =n
   |------------|------------|-------|-------|
 -1/2           0          +1/2      +1    +3/2 = m

   ~------------ -------------------- -------~
     dx'[-1]          dx'[0]           dx'[1]

 The Ex values are positioned on the x-faces of the base
 grid. They represent the Ex field in volumes shifted by
 a half-cell in the x-dimension, as shown here. Only the
 center cell (with width dx'[0]) is fully shown; the
 other two are truncated (shown using >< markers).

 Note that the Ex positions are the in the same positions
 as the previous figure; only the cell boundaries have moved.
 Also note that the points at which Ex is defined are not
 necessarily centered in the volumes they represent; non-
 uniform cell sizes result in off-center volumes like the
 center cell here.

The next figure shows the volumes corresponding to the Hy components, which are shifted in two dimensions (x and z) compared to the base grid.

[figure: Hy volumes]
                                                             p=
  z y     mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm   __ +1/2  s
  |/_x   <<           m:                    m:      >>    |       |
        < <          m :                   m :     > >    |       | dz'[1]
       <  <         m  :                  m  :    >  >    |       |
     Hy........... m........Hy...........m......Hy   >    |       |
     <    <       m    :                m    :  >    >    |       |
    <     ______ m_____:_______________m_____:_>______   __ 0
   <      <     m     /:              m     / >      >    |       |
  mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm       >    |       |
  <       <    |    /  :             |    /  >       >    |       | dz'[0]
  <       <    |   /   :             |   /   >       >    |       |
  <       <    |  /    :             |  /    >       >    |       |
  <       wwwww|w/wwwwwwwwwwwwwwwwwww|w/wwwww>wwwwwwww   __       s
  <      <     |/     w              |/     w>      >    /         /
  _____________|_____________________|________     >    /         /
  <    <       |    w                |    w  >    >    /         /
  <  Hy........|...w........Hy.......|...w...>..Hy    _ 0       / dy[0]
  < <          |  w                  |  w    >  >    /         /
  <<           | w                   | w     > >    /         /
  <            |w                    |w      >>    /         /
  wwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww   __ -1/2    /

  |------------|------------|--------|-------|
-1/2           0          +1/2      +1     +3/2 = m

  ~------------ --------------------- -------~
     dx'[-1]            dx'[0]         dx'[1]

 The Hy values are positioned on the y-edges of the base
 grid. Again here, the 'Hy' labels represent the same points
 as in the basic grid figure above; the edges have shifted
 by a half-cell along the x- and z-axes.

 The grid lines _|:/ are edges of the area represented by
 each Hy value, and the lines drawn using <m>.w represent
 edges where a cell's faces extend beyond the drawn area
 (i.e. where the drawing is truncated in the x- or z-
 directions).

Datastructure: dx_lists_t¶

In this documentation, the E fore-vectors are placed on the base grid. An equivalent formulation could place the H back-vectors on the base grid instead. However, in the case of a non-uniform grid, the operation to get from the "base" cell widths to "derived" ones is not its own inverse.

The base grid's cell sizes could be fully described by a list of three 1D arrays, specifying the cell widths along all three axes:

[dx, dy, dz] = [[dx[0], dx[1], ...], [dy[0], ...], [dz[0], ...]]

Note that this is a list-of-arrays rather than a 2D array, as the simulation domain may have a different number of cells along each axis.

Knowing the base grid's cell widths and the boundary conditions (periodic unless otherwise noted) is enough information to calculate the cell widths dx', dy', and dz' for the derived grids.

However, since most operations are trivially generalized to allow either E or H to be defined on the base grid, they are written to take the a full set of base and derived cell widths, distinguished by which field they apply to rather than their "base" or "derived" status. This removes the need for each function to generate the derived widths, and makes the "base" vs "derived" distinction unnecessary in the code.

The resulting data structure containing all the cell widths takes the form of a list-of-lists-of-arrays. The first list-of-arrays provides the cell widths for the E-field fore-vectors, while the second list-of-arrays does the same for the H-field back-vectors:

 [[[dx_e[0], dx_e[1], ...], [dy_e[0], ...], [dz_e[0], ...]],
  [[dx_h[0], dx_h[1], ...], [dy_h[0], ...], [dz_h[0], ...]]]

where dx_e[0] is the x-width of the m=0 cells, as used when calculating dE/dx, and dy_h[0] is the y-width of the n=0 cells, as used when calculating dH/dy, etc.

Permittivity and Permeability¶

Since each vector component of E and H is defined in a different location and represents a different volume, the value of the spatially-discrete epsilon and mu can also be different for all three field components, even when representing a simple planar interface between two isotropic materials.

As a result, epsilon and mu are taken to have the same dimensions as the field, and composed of the three diagonal tensor components:

[equations: epsilon_and_mu]
epsilon = [epsilon_xx, epsilon_yy, epsilon_zz]
mu = [mu_xx, mu_yy, mu_zz]

or

\[ \epsilon = \begin{bmatrix} \epsilon_{xx} & 0 & 0 \\ 0 & \epsilon_{yy} & 0 \\ 0 & 0 & \epsilon_{zz} \end{bmatrix} $$ $$ \mu = \begin{bmatrix} \mu_{xx} & 0 & 0 \\ 0 & \mu_{yy} & 0 \\ 0 & 0 & \mu_{zz} \end{bmatrix} \]

where the off-diagonal terms (e.g. epsilon_xy) are assumed to be zero.

High-accuracy volumetric integration of shapes on multiple grids can be performed by the gridlock module.

The values of the vacuum permittivity and permability effectively become scaling factors that appear in several locations (e.g. between the E and H fields). In order to limit floating-point inaccuracy and simplify calculations, they are often set to 1 and relative permittivities and permeabilities are used in their places; the true values can be multiplied back in after the simulation is complete if non- normalized results are needed.

Functional and sparse operators¶

meanas.fdmath.functional ¶

Math functions for finite difference simulations

Basic discrete calculus etc.

deriv_forward ¶

deriv_forward(
    dx_e: Sequence[NDArray[floating | complexfloating]]
    | None = None,
) -> tuple[
    fdfield_updater_t, fdfield_updater_t, fdfield_updater_t
]

Utility operators for taking discretized derivatives (backward variant).

Parameters:

Name Type Description Default
dx_e Sequence[NDArray[floating | complexfloating]] | None

Lists of cell sizes for all axes [[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...].

 None

Returns:

Type Description
tuple[fdfield_updater_t, fdfield_updater_t, fdfield_updater_t]

List of functions for taking forward derivatives along each axis.

deriv_back ¶

deriv_back(
    dx_h: Sequence[NDArray[floating | complexfloating]]
    | None = None,
) -> tuple[
    fdfield_updater_t, fdfield_updater_t, fdfield_updater_t
]

Utility operators for taking discretized derivatives (forward variant).

Parameters:

Name Type Description Default
dx_h Sequence[NDArray[floating | complexfloating]] | None

Lists of cell sizes for all axes [[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...].

 None

Returns:

Type Description
tuple[fdfield_updater_t, fdfield_updater_t, fdfield_updater_t]

List of functions for taking forward derivatives along each axis.

curl_forward ¶

curl_forward(
    dx_e: Sequence[NDArray[floating | complexfloating]]
    | None = None,
) -> Callable[[TT], TT]

Curl operator for use with the E field.

Parameters:

Name Type Description Default
dx_e Sequence[NDArray[floating | complexfloating]] | None

Lists of cell sizes for all axes [[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...].

 None

Returns:

Type Description
Callable[[TT], TT]

Function f for taking the discrete forward curl of a field,
Callable[[TT], TT]

f(E) -> curlE \(= \nabla_f \times E\)

curl_back ¶

curl_back(
    dx_h: Sequence[NDArray[floating | complexfloating]]
    | None = None,
) -> Callable[[TT], TT]

Create a function which takes the backward curl of a field.

Parameters:

Name Type Description Default
dx_h Sequence[NDArray[floating | complexfloating]] | None

Lists of cell sizes for all axes [[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...].

 None

Returns:

Type Description
Callable[[TT], TT]

Function f for taking the discrete backward curl of a field,
Callable[[TT], TT]

f(H) -> curlH \(= \nabla_b \times H\)

meanas.fdmath.operators ¶

Matrix operators for finite difference simulations

Basic discrete calculus etc.

shift_circ ¶

shift_circ(
    axis: int, shape: Sequence[int], shift_distance: int = 1
) -> sparse.sparray

Utility operator for performing a circular shift along a specified axis by a specified number of elements.

Parameters:

Name Type Description Default
axis int

Axis to shift along. x=0, y=1, z=2

 required
shape Sequence[int]

Shape of the grid being shifted

 required
shift_distance int

Number of cells to shift by. May be negative. Default 1.

 1

Returns:

Type Description
sparray

Sparse matrix for performing the circular shift.

shift_with_mirror ¶

shift_with_mirror(
    axis: int, shape: Sequence[int], shift_distance: int = 1
) -> sparse.sparray

Utility operator for performing an n-element shift along a specified axis, with mirror boundary conditions applied to the cells beyond the receding edge.

Parameters:

Name Type Description Default
axis int

Axis to shift along. x=0, y=1, z=2

 required
shape Sequence[int]

Shape of the grid being shifted

 required
shift_distance int

Number of cells to shift by. May be negative. Default 1.

 1

Returns:

Type Description
sparray

Sparse matrix for performing the shift-with-mirror.

deriv_forward ¶

deriv_forward(
    dx_e: Sequence[NDArray[floating | complexfloating]],
) -> list[sparse.sparray]

Utility operators for taking discretized derivatives (forward variant).

Parameters:

Name Type Description Default
dx_e Sequence[NDArray[floating | complexfloating]]

Lists of cell sizes for all axes [[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...].

 required

Returns:

Type Description
list[sparray]

List of operators for taking forward derivatives along each axis.

deriv_back ¶

deriv_back(
    dx_h: Sequence[NDArray[floating | complexfloating]],
) -> list[sparse.sparray]

Utility operators for taking discretized derivatives (backward variant).

Parameters:

Name Type Description Default
dx_h Sequence[NDArray[floating | complexfloating]]

Lists of cell sizes for all axes [[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...].

 required

Returns:

Type Description
list[sparray]

List of operators for taking forward derivatives along each axis.

cross ¶

cross(B: Sequence[sparray]) -> sparse.sparray

Cross product operator

Parameters:

Name Type Description Default
B Sequence[sparray]

List [Bx, By, Bz] of sparse matrices corresponding to the x, y, z portions of the operator on the left side of the cross product.

 required

Returns:

Type Description
sparray

Sparse matrix corresponding to (B x), where x is the cross product.

vec_cross ¶

vec_cross(b: vfdfield_t) -> sparse.sparray

Vector cross product operator

Parameters:

Name Type Description Default
b vfdfield_t

Vector on the left side of the cross product.

 required

Returns:

Sparse matrix corresponding to (b x), where x is the cross product.

avg_forward ¶

avg_forward(
    axis: int, shape: Sequence[int]
) -> sparse.sparray

Forward average operator (x4 = (x4 + x5) / 2)

Parameters:

Name Type Description Default
axis int

Axis to average along (x=0, y=1, z=2)

 required
shape Sequence[int]

Shape of the grid to average

 required

Returns:

Type Description
sparray

Sparse matrix for forward average operation.

avg_back ¶

avg_back(axis: int, shape: Sequence[int]) -> sparse.sparray

Backward average operator (x4 = (x4 + x3) / 2)

Parameters:

Name Type Description Default
axis int

Axis to average along (x=0, y=1, z=2)

 required
shape Sequence[int]

Shape of the grid to average

 required

Returns:

Type Description
sparray

Sparse matrix for backward average operation.

curl_forward ¶

curl_forward(
    dx_e: Sequence[NDArray[floating | complexfloating]],
) -> sparse.sparray

Curl operator for use with the E field.

Parameters:

Name Type Description Default
dx_e Sequence[NDArray[floating | complexfloating]]

Lists of cell sizes for all axes [[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...].

 required

Returns:

Type Description
sparray

Sparse matrix for taking the discretized curl of the E-field

curl_back ¶

curl_back(
    dx_h: Sequence[NDArray[floating | complexfloating]],
) -> sparse.sparray

Curl operator for use with the H field.

Parameters:

Name Type Description Default
dx_h Sequence[NDArray[floating | complexfloating]]

Lists of cell sizes for all axes [[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...].

 required

Returns:

Type Description
sparray

Sparse matrix for taking the discretized curl of the H-field

meanas.fdmath.vectorization ¶

Functions for moving between a vector field (list of 3 ndarrays, [f_x, f_y, f_z]) and a 1D array representation of that field [f_x0, f_x1, f_x2,... f_y0,... f_z0,...]. Vectorized versions of the field use row-major (ie., C-style) ordering.

vec ¶

vec(f: None) -> None

vec(f: fdfield_t) -> vfdfield_t

vec(f: cfdfield_t) -> vcfdfield_t

vec(f: fdfield2_t) -> vfdfield2_t

vec(f: cfdfield2_t) -> vcfdfield2_t

vec(f: fdslice_t) -> vfdslice_t

vec(f: cfdslice_t) -> vcfdslice_t

vec(f: ArrayLike) -> NDArray

vec(
    f: fdfield_t
    | cfdfield_t
    | fdfield2_t
    | cfdfield2_t
    | fdslice_t
    | cfdslice_t
    | ArrayLike
    | None,
) -> (
    vfdfield_t
    | vcfdfield_t
    | vfdfield2_t
    | vcfdfield2_t
    | vfdslice_t
    | vcfdslice_t
    | NDArray
    | None
)

Create a 1D ndarray from a vector field which spans a 1-3D region.

Returns None if called with f=None.

Parameters:

Name Type Description Default
f fdfield_t | cfdfield_t | fdfield2_t | cfdfield2_t | fdslice_t | cfdslice_t | ArrayLike | None

A vector field, e.g. [f_x, f_y, f_z] where each f_ component is a 1- to 3-D ndarray (f_* should all be the same size). Doesn't fail with f=None.

 required

Returns:

Type Description
vfdfield_t | vcfdfield_t | vfdfield2_t | vcfdfield2_t | vfdslice_t | vcfdslice_t | NDArray | None

1D ndarray containing the linearized field (or None)

unvec ¶

unvec(
    v: None, shape: Sequence[int], nvdim: int = 3
) -> None

unvec(
    v: vfdfield_t, shape: Sequence[int], nvdim: int = 3
) -> fdfield_t

unvec(
    v: vcfdfield_t, shape: Sequence[int], nvdim: int = 3
) -> cfdfield_t

unvec(
    v: vfdfield2_t, shape: Sequence[int], nvdim: int = 3
) -> fdfield2_t

unvec(
    v: vcfdfield2_t, shape: Sequence[int], nvdim: int = 3
) -> cfdfield2_t

unvec(
    v: vfdslice_t, shape: Sequence[int], nvdim: int = 3
) -> fdslice_t

unvec(
    v: vcfdslice_t, shape: Sequence[int], nvdim: int = 3
) -> cfdslice_t

unvec(
    v: ArrayLike, shape: Sequence[int], nvdim: int = 3
) -> NDArray

unvec(
    v: vfdfield_t
    | vcfdfield_t
    | vfdfield2_t
    | vcfdfield2_t
    | vfdslice_t
    | vcfdslice_t
    | ArrayLike
    | None,
    shape: Sequence[int],
    nvdim: int = 3,
) -> (
    fdfield_t
    | cfdfield_t
    | fdfield2_t
    | cfdfield2_t
    | fdslice_t
    | cfdslice_t
    | NDArray
    | None
)

Perform the inverse of vec(): take a 1D ndarray and output an nvdim-component field of form e.g. [f_x, f_y, f_z] (nvdim=3) where each of f_* is a len(shape)-dimensional ndarray.

Returns None if called with v=None.

Parameters:

Name Type Description Default
v vfdfield_t | vcfdfield_t | vfdfield2_t | vcfdfield2_t | vfdslice_t | vcfdslice_t | ArrayLike | None

1D ndarray representing a vector field of shape shape (or None)

 required
shape Sequence[int]

shape of the vector field

 required
nvdim int

Number of components in each vector

 3

Returns:

Type Description
fdfield_t | cfdfield_t | fdfield2_t | cfdfield2_t | fdslice_t | cfdslice_t | NDArray | None

[f_x, f_y, f_z] where each f_ is a len(shape) dimensional ndarray (or None)

meanas.fdmath.types ¶

Types shared across multiple submodules

dx_lists_t module-attribute ¶

dx_lists_t = Sequence[
    Sequence[NDArray[floating | complexfloating]]
]

'dxes' datastructure which contains grid cell width information in the following format:

[[[dx_e[0], dx_e[1], ...], [dy_e[0], ...], [dz_e[0], ...]],
 [[dx_h[0], dx_h[1], ...], [dy_h[0], ...], [dz_h[0], ...]]]

where dx_e[0] is the x-width of the x=0 cells, as used when calculating dE/dx, and dy_h[0] is the y-width of the y=0 cells, as used when calculating dH/dy, etc.

dx_lists2_t module-attribute ¶

dx_lists2_t = Sequence[
    Sequence[NDArray[floating | complexfloating]]
]

2D 'dxes' datastructure which contains grid cell width information in the following format:

[[[dx_e[0], dx_e[1], ...], [dy_e[0], ...]],
 [[dx_h[0], dx_h[1], ...], [dy_h[0], ...]]]

where dx_e[0] is the x-width of the x=0 cells, as used when calculating dE/dx, and dy_h[0] is the y-width of the y=0 cells, as used when calculating dH/dy, etc.

dx_lists_mut module-attribute ¶

dx_lists_mut = MutableSequence[
    MutableSequence[NDArray[floating | complexfloating]]
]

Mutable version of dx_lists_t

dx_lists2_mut module-attribute ¶

dx_lists2_mut = MutableSequence[
    MutableSequence[NDArray[floating | complexfloating]]
]

Mutable version of dx_lists2_t

fdfield_updater_t module-attribute ¶

fdfield_updater_t = Callable[..., fdfield_t]

Convenience type for functions which take and return an fdfield_t

cfdfield_updater_t module-attribute ¶

cfdfield_updater_t = Callable[..., cfdfield_t]

Convenience type for functions which take and return an cfdfield_t

fdfield ¶

fdfield = fdfield_t | NDArray[floating]

Vector field with shape (3, X, Y, Z) (e.g. [E_x, E_y, E_z])

vfdfield ¶

vfdfield = vfdfield_t | NDArray[floating]

Linearized vector field (single vector of length 3XY*Z)

cfdfield ¶

cfdfield = cfdfield_t | NDArray[complexfloating]

Complex vector field with shape (3, X, Y, Z) (e.g. [E_x, E_y, E_z])

vcfdfield ¶

vcfdfield = vcfdfield_t | NDArray[complexfloating]

Linearized complex vector field (single vector of length 3XY*Z)

fdslice ¶

fdslice = fdslice_t | NDArray[floating]

Vector field slice with shape (3, X, Y) (e.g. [E_x, E_y, E_z] at a single Z position)

vfdslice ¶

vfdslice = vfdslice_t | NDArray[floating]

Linearized vector field slice (single vector of length 3XY)

cfdslice ¶

cfdslice = cfdslice_t | NDArray[complexfloating]

Complex vector field slice with shape (3, X, Y) (e.g. [E_x, E_y, E_z] at a single Z position)

vcfdslice ¶

vcfdslice = vcfdslice_t | NDArray[complexfloating]

Linearized complex vector field slice (single vector of length 3XY)

fdfield2 ¶

fdfield2 = fdfield2_t | NDArray[floating]

2D Vector field with shape (2, X, Y) (e.g. [E_x, E_y])

vfdfield2 ¶

vfdfield2 = vfdfield2_t | NDArray[floating]

2D Linearized vector field (single vector of length 2XY)

cfdfield2 ¶

cfdfield2 = cfdfield2_t | NDArray[complexfloating]

2D Complex vector field with shape (2, X, Y) (e.g. [E_x, E_y])

vcfdfield2 ¶

vcfdfield2 = vcfdfield2_t | NDArray[complexfloating]

2D Linearized complex vector field (single vector of length 2XY)