Solvers for eigenvalue / eigenvector problems
+ + + + + + + + + + +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(
+ 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 |
+
+ 1e-13
+ |
+
+ solver
+ |
+
+ Callable[..., NDArray[complex128]] | None
+ |
+
+
+
+ Solver function of the form |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ tuple[complex, NDArray[complex128]]
+ |
+
+
+
+ (eigenvalues, eigenvectors) + |
+
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]
+ |
+
+
+
+
|
+
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
resulting in
+Maxwell's equations are then
+With \(\Delta_t \to 0\), this simplifies to
+and then
+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(
+ 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 |
+ + required + | +
+ epsilon
+ |
+
+ fdfield
+ |
+
+
+
+ Dielectric constant + |
+ + required + | +
+ mu
+ |
+
+ fdfield | None
+ |
+
+
+
+ Magnetic permeability (default 1 everywhere) + |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ cfdfield_updater_t
+ |
+
+
+
+ Function |
+
+ cfdfield_updater_t
+ |
+
+
+
+
|
+
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 |
+ + 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 |
+
+ Callable[[cfdfield, cfdfield], tuple[cfdfield_t, cfdfield_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 |
+ + required + | +
+ mu
+ |
+
+ fdfield | None
+ |
+
+
+
+ Magnetic permeability (default 1 everywhere) + |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ cfdfield_updater_t
+ |
+
+
+
+ Function |
+
+ 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 |
+ + required + | +
+ mu
+ |
+
+ fdfield | None
+ |
+
+
+
+ Magnetic permeability (default 1 everywhere) + |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ cfdfield_updater_t
+ |
+
+
+
+ Function |
+
+ cfdfield_updater_t
+ |
+
+
+
+
|
+
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
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. |
+ + required + | +
+ omega
+ |
+
+ complex
+ |
+
+
+
+ Angular frequency of the simulation + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists_t
+ |
+
+
+
+ Grid parameters |
+ + required + | +
+ epsilon
+ |
+
+ fdfield
+ |
+
+
+
+ Dielectric constant distribution + |
+ + required + | +
+ mu
+ |
+
+ fdfield | None
+ |
+
+
+
+ Magnetic permeability (default 1 everywhere) + |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ cfdfield_updater_t
+ |
+
+
+
+ Function |
+
+ cfdfield_updater_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(...).
This function also shifts the input E field by one cell as required
+for computing the Poynting cross product (see meanas.fdfd module docs).
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 |
+ + required + | +
Returns:
+| Type | +Description | +
|---|---|
+ Callable[[cfdfield, cfdfield], cfdfield_t]
+ |
+
+
+
+ Function |
+
+ Callable[[cfdfield, cfdfield], cfdfield_t]
+ |
+
+
+
+ For time-average power, call it as |
+
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_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 + +
del x (1/mu * del x) - omega**2 * epsilon
+for use with the E-field, with wave equation + +
(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 |
+ + 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 |
+
+ None
+ |
+
+ pmc
+ |
+
+ vfdfield | None
+ |
+
+
+
+ Vectorized mask specifying PMC cells. Any cells where |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ sparray
+ |
+
+
+
+ Sparse matrix containing the wave operator. + |
+
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 |
+ + required + | +
Returns:
+| Type | +Description | +
|---|---|
+ tuple[sparray, sparray]
+ |
+
+
+
+ Preconditioner matrices |
+
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 + +
del x (1/epsilon * del x) - omega**2 * mu
+for use with the H-field, with wave equation + +
(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 |
+ + 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 |
+
+ None
+ |
+
+ pmc
+ |
+
+ vfdfield | None
+ |
+
+
+
+ Vectorized mask specifying PMC cells. Any cells where |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ sparray
+ |
+
+
+
+ Sparse matrix containing the wave operator. + |
+
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
[[-i * omega * epsilon, del x ],
+ [del x, i * omega * mu]]
+for use with a field vector of the form cat(vec(E), vec(H)):
Parameters:
+| Name | +Type | +Description | +Default | +
|---|---|---|---|
+ omega
+ |
+
+ complex
+ |
+
+
+
+ Angular frequency of the simulation + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists_t
+ |
+
+
+
+ Grid parameters |
+ + 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 |
+
+ None
+ |
+
+ pmc
+ |
+
+ vfdfield | None
+ |
+
+
+
+ Vectorized mask specifying PMC cells. Any cells where |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ sparray
+ |
+
+
+
+ Sparse matrix containing the wave operator. + |
+
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 |
+ + required + | +
+ mu
+ |
+
+ vfdfield | None
+ |
+
+
+
+ Vectorized magnetic permeability (default 1 everywhere) + |
+
+ None
+ |
+
+ pmc
+ |
+
+ vfdfield | None
+ |
+
+
+
+ Vectorized mask specifying PMC cells. Any cells where |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ sparray
+ |
+
+
+
+ Sparse matrix for converting E to H. + |
+
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 |
+ + required + | +
+ mu
+ |
+
+ vfdfield | None
+ |
+
+
+
+ Vectorized magnetic permeability (default 1 everywhere) + |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ sparray
+ |
+
+
+
+ Sparse matrix for converting M to J. + |
+
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 |
+ + required + | +
Returns:
+| Type | +Description | +
|---|---|
+ sparray
+ |
+
+
+
+ Sparse matrix containing the |
+
+ sparray
+ |
+
+
+
+ cross product. + |
+
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,
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 |
+ + required + | +
Returns:
+| Type | +Description | +
|---|---|
+ sparray
+ |
+
+
+
+ Sparse matrix containing the |
+
+ sparray
+ |
+
+
+
+ cross product. + |
+
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
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 |
+ + 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(
+ 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 |
+ + 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. + |
+
Solvers and solver interface for FDFD problems.
+ + + + + + + + + + +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
+ |
+
+
+
+
|
+ + 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 |
+
+ _scipy_qmr
+ |
+
+ matrix_solver_opts
+ |
+
+ dict[str, Any] | None
+ |
+
+
+
+ Passed as kwargs to |
+
+ None
+ |
+
+ E_guess
+ |
+
+ vcfdfield | None
+ |
+
+
+
+ Guess at the solution E-field. |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ vcfdfield_t
+ |
+
+
+
+ E-field which solves the system. + |
+
Functions for creating stretched coordinate perfectly matched layer (PML) absorbers.
+ + + + + + + + + + +module-attribute
+
+
+¶Typedef for s-functions, see prepare_s_function()
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 |
+
+ s_function_t
+ |
+
+
+
+ before use. + |
+
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]
+ |
+
+
+
+
|
+ + 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 |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ list[list[NDArray[float64]]]
+ |
+
+
+
+ Complex cell widths (dx_lists_mut) as discussed in |
+
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 |
+ + 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 |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ list[list[NDArray[float64]]]
+ |
+
+
+
+ Complex cell widths (dx_lists_mut) as discussed in |
+
+ list[list[NDArray[float64]]]
+ |
+
+
+
+ Multiple calls to this function may be necessary if multiple absorpbing boundaries are needed. + |
+
Functions for performing near-to-farfield transformation (and the reverse).
+ + + + + + + + + + +near_to_farfield(
+ E_near: transverse_slice_pair,
+ H_near: transverse_slice_pair,
+ dx: float,
+ dy: float,
+ padded_size: Sequence[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
+ |
+
+ transverse_slice_pair
+ |
+
+
+
+ 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
+ |
+
+ transverse_slice_pair
+ |
+
+
+
+ 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
+ |
+
+ Sequence[int] | int | None
+ |
+
+
+
+ Shape of the output. A single integer |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ dict[str, Any]
+ |
+
+
+
+ Dict with keys + |
+
+ dict[str, Any]
+ |
+
+
+
+
|
+
+ dict[str, Any]
+ |
+
+
+
+
|
+
+ dict[str, Any]
+ |
+
+
+
+
|
+
+ dict[str, Any]
+ |
+
+
+
+
|
+
+ dict[str, Any]
+ |
+
+
+
+
|
+
+ dict[str, Any]
+ |
+
+
+
+
|
+
far_to_nearfield(
+ E_far: transverse_slice_pair,
+ H_far: transverse_slice_pair,
+ dkx: float,
+ dky: float,
+ padded_size: Sequence[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
+ |
+
+ transverse_slice_pair
+ |
+
+
+
+ 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
+ |
+
+ transverse_slice_pair
+ |
+
+
+
+ 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
+ |
+
+ Sequence[int] | int | None
+ |
+
+
+
+ Shape of the output. A single integer |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ dict[str, Any]
+ |
+
+
+
+ Dict with keys + |
+
+ dict[str, Any]
+ |
+
+
+
+
|
+
+ dict[str, Any]
+ |
+
+
+
+
|
+
+ dict[str, Any]
+ |
+
+
+
+
|
+
Basic discrete calculus for finite difference (fd) simulations.
+Discrete fields are stored in one of two forms:
+fdfield_t form is a multidimensional numpy.NDArrayU[m, n, p], where m, n, and p are
+ discrete indices referring to positions on the x, y, and z axes respectively.E[:, m, n, p] = [Ex[m, n, p], Ey[m, n, p], Ez[m, n, p]].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)meanas.fdmath.functional.
+ The generated functions act on fields in the fdfield_t form.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).
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.
+Define the discrete forward derivative as + +
where \(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 + +
or
+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 + +
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.
Expanding to three dimensions, we can define two gradients
+
+
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
+
+
[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)
+There are also two divergences,
+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
+The two curls are then
+and
+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)
+If we discretize both space (m,n,p) and time (l), Maxwell's equations become
+with
+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, + +
implying that the discrete Maxwell's equations do not produce spurious charges.
+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:
+We can substitute in a time-harmonic fields
+resulting in
+This gives the frequency-domain wave equation,
+With uniform material distribution and no sources
+the frequency domain wave equation simplifies to
+Since \(\hat{\nabla} \cdot \tilde{E}_{\vec{r}} = 0\), we can simplify
+and we get
+We can convert this to three scalar-wave equations of the form
+with \(K^2 = \Omega^2 \mu \epsilon\). Now we let
+resulting in
+with similar expressions for the y and z dimnsions (and \(K_y, K_z\)).
+This implies
+where \(c = \sqrt{\mu \epsilon}\).
+Assuming real \((k_x, k_y, k_z), \omega\) will be real only if
+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).
+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).
+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.
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
+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.
+ + + + + + + + + + +Math functions for finite difference simulations
+Basic discrete calculus etc.
+ + + + + + + + + + +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
+ |
+
+ 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(
+ 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
+ |
+
+ 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(
+ 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
+ |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ Callable[[TT], TT]
+ |
+
+
+
+ Function |
+
+ Callable[[TT], TT]
+ |
+
+
+
+
|
+
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
+ |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ Callable[[TT], TT]
+ |
+
+
+
+ Function |
+
+ Callable[[TT], TT]
+ |
+
+
+
+
|
+
Matrix operators for finite difference simulations
+Basic discrete calculus etc.
+ + + + + + + + + + +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(
+ 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. + |
+
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
+ |
+ + required + | +
Returns:
+| Type | +Description | +
|---|---|
+ list[sparray]
+ |
+
+
+
+ List of operators for taking forward derivatives along each axis. + |
+
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
+ |
+ + required + | +
Returns:
+| Type | +Description | +
|---|---|
+ list[sparray]
+ |
+
+
+
+ List of operators for taking forward derivatives along each axis. + |
+
Cross product operator
+ + +Parameters:
+| Name | +Type | +Description | +Default | +
|---|---|---|---|
+ B
+ |
+
+ Sequence[sparray]
+ |
+
+
+
+ List |
+ + required + | +
Returns:
+| Type | +Description | +
|---|---|
+ sparray
+ |
+
+
+
+ Sparse matrix corresponding to (B x), where x is the cross product. + |
+
Vector cross product operator
+ + +Parameters:
+| Name | +Type | +Description | +Default | +
|---|---|---|---|
+ b
+ |
+
+ vfdfield
+ |
+
+
+
+ Vector on the left side of the cross product. + |
+ + required + | +
Returns:
+Sparse matrix corresponding to (b x), where x is the cross product.
+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. + |
+
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 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
+ |
+ + required + | +
Returns:
+| Type | +Description | +
|---|---|
+ sparray
+ |
+
+
+
+ Sparse matrix for taking the discretized curl of the E-field + |
+
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
+ |
+ + required + | +
Returns:
+| Type | +Description | +
|---|---|
+ sparray
+ |
+
+
+
+ Sparse matrix for taking the discretized curl of the H-field + |
+
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(
+ 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. |
+ + 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 |
+
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
+ |
+
+
+
+
|
+
Types shared across multiple submodules
+ + + + + + + + + + +module-attribute
+
+
+¶'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.
module-attribute
+
+
+¶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.
module-attribute
+
+
+¶Mutable version of dx_lists_t
module-attribute
+
+
+¶Mutable version of dx_lists2_t
module-attribute
+
+
+¶Convenience type for functions which take and return a real fdfield
module-attribute
+
+
+¶Convenience type for functions which take and return a complex cfdfield
Vector field with shape (3, X, Y, Z) (e.g. [E_x, E_y, E_z])
Linearized vector field (single vector of length 3XY*Z)
+Complex vector field with shape (3, X, Y, Z) (e.g. [E_x, E_y, E_z])
Linearized complex vector field (single vector of length 3XY*Z)
+Vector field slice with shape (3, X, Y) (e.g. [E_x, E_y, E_z] at a single Z position)
Linearized vector field slice (single vector of length 3XY)
+Complex vector field slice with shape (3, X, Y) (e.g. [E_x, E_y, E_z] at a single Z position)
Linearized complex vector field slice (single vector of length 3XY)
+2D Vector field with shape (2, X, Y) (e.g. [E_x, E_y])
2D Linearized vector field (single vector of length 2XY)
+2D Complex vector field with shape (2, X, Y) (e.g. [E_x, E_y])
2D Linearized complex vector field (single vector of length 2XY)
+Utilities for running finite-difference time-domain (FDTD) simulations
+See the discussion of Maxwell's Equations in meanas.fdmath for basic
+mathematical background.
From the discussion of "Plane waves and the Dispersion relation" in meanas.fdmath,
+we have
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.
Let
+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
+where in the last line the spatial subscripts have been dropped to emphasize +the time subscripts \(l, l'\), i.e.
+etc. +For \(l' = l + \frac{1}{2}\) we get
+and for \(l' = l - \frac{1}{2}\),
+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:
+Rewriting the Poynting theorem in terms of the energy expressions,
+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\)).
+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. temporal_phasor(...) and temporal_phasor_scale(...)
+apply the same Fourier sum to a scalar waveform, which is useful when a pulsed
+source envelope must be normalized before being applied to a point source or
+mode source. real_injection_scale(...) is the corresponding helper for the
+common real-valued injection pattern numpy.real(scale * waveform). Convenience
+wrappers accumulate_phasor_e, accumulate_phasor_h, and accumulate_phasor_j
+apply the usual Yee time offsets. reconstruct_real(...) and the corresponding
+reconstruct_real_e/h/j(...) wrappers perform the inverse operation, converting
+phasors back into real-valued field snapshots at explicit Yee-aligned times.
+For scalar omega, the reconstruction helpers accept either a plain field
+phasor or the batched (1, *sample_shape) form used by the accumulator helpers.
+The helpers accumulate
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.
For FDTD/FDFD equivalence, this phasor path is the primary comparison workflow.
+It isolates the guided +\omega response that the frequency-domain solver
+targets directly, regardless of whether the underlying FDTD run used real- or
+complex-valued fields.
For exact pulsed FDTD/FDFD equivalence it is often simplest to keep the
+injected source, fields, and CPML auxiliary state complex-valued. The
+real_injection_scale(...) helper is instead for the more ordinary one-run
+real-valued source path, where the intended positive-frequency waveform is
+injected through numpy.real(scale * waveform) and any remaining negative-
+frequency leakage is controlled by the pulse bandwidth and run length.
reconstruct_real(...) is for a different question: given a phasor, what late
+real-valued field snapshot should it produce? That raw-snapshot comparison is
+stricter and noisier because a monitor plane generally contains both the guided
+field and the remaining orthogonal content,
Phasor/modal comparisons mostly validate the guided +\omega term. Raw
+real-field comparisons expose both terms at once, so they should be treated as
+secondary diagnostics rather than the main solver-equivalence benchmark.
The Ricker wavelet (normalized second derivative of a Gaussian) is commonly used for the pulse +shape. It can be written
+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).
+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:
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:
+dxes with the matching SCPML transform.dxes.The electric-current sign convention used throughout the examples and tests is
+which matches the FDFD right-hand side
+Basic FDTD field updates
+ + + + + + + + + + +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
dxes: "Datastructure: dx_lists_t" sectionepsilon: "Permittivity and Permeability" sectionAlso see the "Timestep" section of meanas.fdtd for a discussion of
+the dt parameter.
Parameters:
+| Name | +Type | +Description | +Default | +
|---|---|---|---|
+ dt
+ |
+
+ float
+ |
+
+
+
+ Timestep. See |
+ + required + | +
+ dxes
+ |
+
+ dx_lists_t | None
+ |
+
+
+
+ Grid description; see |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ fdfield_updater_t
+ |
+
+
+
+ Function |
+
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
dxes: "Datastructure: dx_lists_t" sectionmu: "Permittivity and Permeability" sectionAlso see the "Timestep" section of meanas.fdtd for a discussion of
+the dt parameter.
Parameters:
+| Name | +Type | +Description | +Default | +
|---|---|---|---|
+ dt
+ |
+
+ float
+ |
+
+
+
+ Timestep. See |
+ + required + | +
+ dxes
+ |
+
+ dx_lists_t | None
+ |
+
+
+
+ Grid description; see |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ fdfield_updater_t
+ |
+
+
+
+ Function |
+
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:
+cpml_params[axis][polarity_index] table with cpml_params(...).updates_with_cpml(...) together with dt, dxes, and
+ epsilon.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(
+ 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 ( |
+ + required + | +
+ polarity
+ |
+
+ int
+ |
+
+
+
+ Which face along that axis ( |
+ + 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 |
+
+ 3.5
+ |
+
+ ma
+ |
+
+ float
+ |
+
+
+
+ Polynomial grading exponent for the complex-frequency shift |
+
+ 1
+ |
+
+ cfs_alpha
+ |
+
+ float
+ |
+
+
+
+ Maximum complex-frequency shift parameter. + |
+
+ 0
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ dict[str, Any]
+ |
+
+
+
+ Dictionary with: + |
+
+ dict[str, Any]
+ |
+
+
+
+
|
+
+ dict[str, Any]
+ |
+
+
+
+
|
+
+ dict[str, Any]
+ |
+
+
+
+
|
+
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[..., None], Callable[..., 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 |
+ + required + | +
+ dt
+ |
+
+ float
+ |
+
+
+
+ Timestep. + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists_t
+ |
+
+
+
+ Yee-grid spacing lists |
+ + 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[..., None]
+ |
+
+
+
+
|
+
+ Callable[..., None]
+ |
+
+
+
+ Yee updates: + |
+
+ tuple[Callable[..., None], Callable[..., None]]
+ |
+
+
+
+
|
+
+ tuple[Callable[..., None], Callable[..., None]]
+ |
+
+
+
+
|
+
+ tuple[Callable[..., None], Callable[..., None]]
+ |
+
+
+
+ The closures retain the CPML auxiliary state internally. + |
+
Boundary conditions
+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
+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 |
+ + required + | +
+ dxes
+ |
+
+ dx_lists_t | None
+ |
+
+
+
+ Grid description; see |
+
+ 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(
+ 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(
+ 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 |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ fdfield_t
+ |
+
+
+
+ Energy, at the time of the H-field |
+
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 |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ fdfield_t
+ |
+
+
+
+ Energy, at the time of the E-field |
+
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 |
+ + 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 |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ fdfield_t
+ |
+
+
+
+ Change in energy from the time of |
+
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 |
+ + 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 |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ fdfield_t
+ |
+
+
+
+ Change in energy from the time of |
+
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(
+ 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 |
+ + required + | +
+ hh
+ |
+
+ fdfield
+ |
+
+
+
+ Three-component magnetic-field product, such as |
+ + required + | +
+ epsilon
+ |
+
+ fdfield | float | None
+ |
+
+
+
+ Electric material weight; defaults to |
+
+ None
+ |
+
+ mu
+ |
+
+ fdfield | float | None
+ |
+
+
+
+ Magnetic material weight; defaults to |
+
+ 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. + |
+
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_laccumulate_phasor_h(..., step=l) for H_{l + 1/2}accumulate_phasor_j(..., step=l) for J_{l + 1/2}temporal_phasor(...) and temporal_phasor_scale(...) apply the same Fourier
+sum to a 1D scalar waveform. They are useful for normalizing a pulsed source
+before that scalar waveform is applied to a point source or spatial mode source.
+real_injection_scale(...) is a companion helper for the common real-valued
+injection pattern numpy.real(scale * waveform), where waveform is the
+analytic positive-frequency signal and the injected real current should still
+produce a desired phasor response.
+reconstruct_real(...) and its E/H/J wrappers perform the inverse operation:
+they turn one or more phasors back into real-valued field snapshots at explicit
+Yee-aligned sample times. For a scalar target frequency they accept either a
+plain field phasor or the batched (1, *sample_shape) form used elsewhere in
+this module.
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(
+ 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.
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.
temporal_phasor(
+ samples: ArrayLike,
+ omegas: float
+ | complex
+ | Sequence[float | complex]
+ | NDArray,
+ dt: float,
+ *,
+ start_step: int = 0,
+ offset_steps: float = 0.0,
+) -> NDArray[numpy.complexfloating]
+Fourier-project a 1D temporal waveform onto one or more angular frequencies.
+The returned quantity is
+dt * sum(exp(-1j * omega * t_step) * samples[step_index])
+where t_step = (start_step + step_index + offset_steps) * dt.
temporal_phasor_scale(
+ samples: ArrayLike,
+ omegas: float
+ | complex
+ | Sequence[float | complex]
+ | NDArray,
+ dt: float,
+ *,
+ start_step: int = 0,
+ offset_steps: float = 0.0,
+ target: ArrayLike = 1.0,
+) -> NDArray[numpy.complexfloating]
+Return the scalar multiplier that gives a desired temporal phasor response.
+The returned scale satisfies
+temporal_phasor(scale * samples, omegas, dt, ...) == target
+for each target frequency. The result keeps a leading frequency axis even
+when omegas is scalar.
real_injection_scale(
+ samples: ArrayLike,
+ omegas: float
+ | complex
+ | Sequence[float | complex]
+ | NDArray,
+ dt: float,
+ *,
+ start_step: int = 0,
+ offset_steps: float = 0.0,
+ target: ArrayLike = 1.0,
+) -> NDArray[numpy.complexfloating]
+Return the scale for a real-valued injection built from an analytic waveform.
+If the time-domain source is applied as
+numpy.real(scale * samples[step])
+then the desired positive-frequency phasor is obtained by compensating for +the 1/2 factor between the real-valued source and its analytic component:
+scale = 2 * target / temporal_phasor(samples, ...)
+This helper normalizes only the intended positive-frequency component. Any +residual negative-frequency leakage is controlled by the waveform design and +the accumulation window.
+ + +reconstruct_real(
+ phasors: ArrayLike,
+ omegas: float
+ | complex
+ | Sequence[float | complex]
+ | NDArray,
+ dt: float,
+ step: int,
+ *,
+ offset_steps: float = 0.0,
+) -> NDArray[numpy.floating]
+Reconstruct a real-valued field snapshot from one or more phasors.
+The returned quantity is
+real(phasor * exp(1j * omega * t_step))
+where t_step = (step + offset_steps) * dt.
For multi-frequency inputs, the leading frequency axis is preserved. For a
+scalar omega, callers may pass either (1, *sample_shape) or
+sample_shape; the return shape matches that choice.
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(
+ 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(
+ 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}.
reconstruct_real_e(
+ phasors: ArrayLike,
+ omegas: float
+ | complex
+ | Sequence[float | complex]
+ | NDArray,
+ dt: float,
+ step: int,
+) -> NDArray[numpy.floating]
+Reconstruct a real E-field snapshot taken at integer timestep step.
reconstruct_real_h(
+ phasors: ArrayLike,
+ omegas: float
+ | complex
+ | Sequence[float | complex]
+ | NDArray,
+ dt: float,
+ step: int,
+) -> NDArray[numpy.floating]
+Reconstruct a real H-field snapshot corresponding to H_{step + 1/2}.
reconstruct_real_j(
+ phasors: ArrayLike,
+ omegas: float
+ | complex
+ | Sequence[float | complex]
+ | NDArray,
+ dt: float,
+ step: int,
+) -> NDArray[numpy.floating]
+Reconstruct a real current snapshot corresponding to J_{step + 1/2}.
The package is documented directly from its docstrings. The most useful entry +points are:
+The waveguide and FDTD pages are the best places to start if you want the +mathematical derivations rather than just the callable reference.
+ + + + + + + + + + + + + +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
+Expanding the first two equations into vector components, we get
+Substituting in our expressions for \(\vec{E}\), \(\vec{H}\) and discretizing:
+Rewrite the last three equations as
+Now apply \(\imath \beta \tilde{\partial}_x\) to the last equation, +then substitute in for \(\imath \beta H_x\) and \(\imath \beta H_y\):
+With a similar approach (but using \(\imath \beta \tilde{\partial}_y\) instead), we can get
+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
+and
+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
+and, similarly,
+By combining both pairs of expressions, we get
+Using these, we can construct the eigenvalue problem
+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(
+ 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
+\(\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 |
+ + 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(
+ 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
+\(\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 |
+ + 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(
+ 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 |
+ + required + | +
+ omega
+ |
+
+ complex
+ |
+
+
+
+ The angular frequency of the system + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + required + | +
+ epsilon
+ |
+
+ vfdslice
+ |
+
+
+
+ Vectorized dielectric constant grid + |
+ + required + | +
+ mu
+ |
+
+ vfdslice | None
+ |
+
+
+
+ Vectorized magnetic permeability grid (default 1 everywhere) + |
+
+ None
+ |
+
+ prop_phase
+ |
+
+ float
+ |
+
+
+
+ Phase shift |
+
+ 0
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ vcfdslice_t
+ |
+
+
+
+
|
+
+ vcfdslice_t
+ |
+
+
+
+ and contains all three vector components. + |
+
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
+so the returned fields represent a unit-power forward mode under the +discrete Yee-grid Poynting inner product.
+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 |
+ + required + | +
+ omega
+ |
+
+ complex
+ |
+
+
+
+ The angular frequency of the system + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + required + | +
+ epsilon
+ |
+
+ vfdslice
+ |
+
+
+
+ Vectorized dielectric constant grid + |
+ + required + | +
+ mu
+ |
+
+ vfdslice | None
+ |
+
+
+
+ Vectorized magnetic permeability grid (default 1 everywhere) + |
+
+ None
+ |
+
+ prop_phase
+ |
+
+ float
+ |
+
+
+
+ Phase shift |
+
+ 0
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ vcfdslice_t
+ |
+
+
+
+
|
+
+ vcfdslice_t
+ |
+
+
+
+ and contains all three vector components. + |
+
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(
+ 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 |
+ + required + | +
+ omega
+ |
+
+ complex
+ |
+
+
+
+ The angular frequency of the system + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + 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(
+ 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 |
+ + required + | +
+ omega
+ |
+
+ complex
+ |
+
+
+
+ The angular frequency of the system + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + 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(
+ 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 |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + required + | +
+ mu
+ |
+
+ vfdslice | None
+ |
+
+
+
+ Vectorized magnetic permeability grid (default 1 everywhere) + |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ sparray
+ |
+
+
+
+ Sparse matrix representing the operator. + |
+
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
+as well as the intermediate equations
+Combining these, we get
+Parameters:
+| Name | +Type | +Description | +Default | +
|---|---|---|---|
+ wavenumber
+ |
+
+ complex
+ |
+
+
+
+ Wavenumber assuming fields have z-dependence of |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + required + | +
+ epsilon
+ |
+
+ vfdslice
+ |
+
+
+
+ Vectorized dielectric constant grid + |
+ + required + | +
Returns:
+| Type | +Description | +
|---|---|
+ sparray
+ |
+
+
+
+ Sparse matrix representing the operator. + |
+
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 |
+ + required + | +
+ omega
+ |
+
+ complex
+ |
+
+
+
+ The angular frequency of the system + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + required + | +
+ mu
+ |
+
+ vfdslice | None
+ |
+
+
+
+ Vectorized magnetic permeability grid (default 1 everywhere) + |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ sparray
+ |
+
+
+
+ Sparse matrix representation of the operator. + |
+
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 |
+ + required + | +
+ omega
+ |
+
+ complex
+ |
+
+
+
+ The angular frequency of the system + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + required + | +
+ epsilon
+ |
+
+ vfdslice
+ |
+
+
+
+ Vectorized dielectric constant grid + |
+ + required + | +
Returns:
+| Type | +Description | +
|---|---|
+ sparray
+ |
+
+
+
+ Sparse matrix representation of the operator. + |
+
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 |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + required + | +
Returns:
+| Type | +Description | +
|---|---|
+ sparray
+ |
+
+
+
+ Sparse matrix representation of the operator. + |
+
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 |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + required + | +
Returns:
+| Type | +Description | +
|---|---|
+ sparray
+ |
+
+
+
+ Sparse matrix representation of the operator. + |
+
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 |
+ + required + | +
+ omega
+ |
+
+ complex
+ |
+
+
+
+ The angular frequency of the system + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + 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 |
+
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 |
+ + required + | +
+ omega
+ |
+
+ complex
+ |
+
+
+
+ The angular frequency of the system + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + 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 |
+
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.
An adjoint approach is used to calculate the sensitivity; the derivation is provided here:
+Starting with the eigenvalue equation
+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\):
We then multiply by \(H_{yx}^\star = \begin{bmatrix}H_y^\star \\ -H_x^\star \end{bmatrix}\) from the left:
+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,
and we can simplify to
+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 |
+ + required + | +
+ h_norm
+ |
+
+ vcfdslice
+ |
+
+
+
+ Normalized, vectorized H_xyz field for the mode. E.g. as returned by |
+ + 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 |
+ + 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(
+ 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 |
+ + required + | +
+ epsilon
+ |
+
+ vfdslice
+ |
+
+
+
+ Dielectric constant + |
+ + required + | +
+ mu
+ |
+
+ vfdslice | None
+ |
+
+
+
+ Magnetic permeability (default 1 everywhere) + |
+
+ None
+ |
+
+ mode_margin
+ |
+
+ int
+ |
+
+
+
+ The eigensolver will actually solve for |
+
+ 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 + |
+
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 |
+
+ ()
+ |
+
+ **kwargs
+ |
+
+ Any
+ |
+
+
+
+ passed to |
+
+ {}
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ tuple[vcfdfield2_t, complex]
+ |
+
+
+
+ (e_xy, wavenumber) + |
+
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,
+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 |
+ + required + | +
+ h2
+ |
+
+ vcfdfield2
+ |
+
+
+
+ Vectorized magnetic field, typically from |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Two-dimensional Yee-grid spacing lists |
+ + 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 |
+
+ 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. + |
+
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:
+solve_mode(...).compute_source(...).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 axisThat same convention controls which side of the selected slice is used for the +overlap window and how the expanded field is phased.
+ + + + + + + + + + +solve_mode(
+ mode_number: int,
+ omega: complex,
+ dxes: dx_lists_t,
+ axis: int,
+ polarity: int,
+ slices: Sequence[slice],
+ epsilon: fdfield,
+ mu: fdfield | None = None,
+) -> Waveguide3DMode
+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 |
+ + 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]
+ |
+
+
+
+
|
+ + required + | +
+ epsilon
+ |
+
+ fdfield
+ |
+
+
+
+ Dielectric constant + |
+ + required + | +
+ mu
+ |
+
+ fdfield | None
+ |
+
+
+
+ Magnetic permeability (default 1 everywhere) + |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ Waveguide3DMode
+ |
+
+
+
+ Dictionary containing: + |
+
+ Waveguide3DMode
+ |
+
+
+
+
|
+
+ Waveguide3DMode
+ |
+
+
+
+
|
+
+ Waveguide3DMode
+ |
+
+
+
+
|
+
+ Waveguide3DMode
+ |
+
+
+
+
|
+
The returned fields are normalized through the waveguide_2d
+normalization convention before being expanded back to 3D.
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 |
+ + 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]
+ |
+
+
+
+
|
+ + required + | +
+ mu
+ |
+
+ fdfield | None
+ |
+
+
+
+ Magnetic permeability (default 1 everywhere) + |
+
+ None
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ cfdfield_t
+ |
+
+
+
+
|
+
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(
+ E: cfdfield,
+ 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
+where the sum is over the full Yee-grid field storage.
+The construction uses a two-cell window immediately upstream of the selected +slice:
+polarity=+1, the two cells just before slices[axis].startpolarity=-1, the two cells just after slices[axis].stopThe 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
+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
+ |
+
+
+
+ E-field of the mode + |
+ + required + | +
+ wavenumber
+ |
+
+ complex
+ |
+
+
+
+ Wavenumber of the mode + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists_t
+ |
+
+
+
+ Grid parameters |
+ + 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]
+ |
+
+
+
+
|
+ + required + | +
Returns:
+| Type | +Description | +
|---|---|
+ cfdfield_t
+ |
+
+
+
+
|
+
+ cfdfield_t
+ |
+
+
+
+ over the retained overlap window. + |
+
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 |
+ + 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]
+ |
+
+
+
+
|
+ + required + | +
Returns:
+| Type | +Description | +
|---|---|
+ cfdfield_t
+ |
+
+
+
+
|
+
This helper assumes that the waveguide cross-section remains constant +along the propagation axis and applies the phase factor
+to each copied slice.
+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
where m is the angular wavenumber returned by solve_mode(s). It is often
+convenient to introduce the corresponding linear wavenumber
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
and from them the diagonal metric matrices
+With the same forward/backward derivative notation used in waveguide_2d, the
+coordinate-transformed discrete curl equations used here are
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
This yields the transverse electric eigenproblem implemented by
+cylindrical_operator(...):
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(
+ omega: float,
+ dxes: dx_lists2_t,
+ epsilon: vfdslice,
+ rmin: float,
+) -> sparse.sparray
+Cylindrical coordinate waveguide operator of the form
+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 |
+ + 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(
+ mode_numbers: Sequence[int],
+ omega: float,
+ dxes: dx_lists2_t,
+ epsilon: vfdslice,
+ rmin: float,
+ mode_margin: int = 2,
+) -> tuple[
+ NDArray[numpy.complex128], NDArray[numpy.complex128]
+]
+Given a 2d (r, y) slice of epsilon, attempts to solve for the eigenmode + of the bent waveguide with the specified mode number.
+ + +Parameters:
+| Name | +Type | +Description | +Default | +
|---|---|---|---|
+ mode_numbers
+ |
+
+ Sequence[int]
+ |
+
+
+
+ Mode numbers to solve, 0-indexed. + |
+ + required + | +
+ omega
+ |
+
+ float
+ |
+
+
+
+ Angular frequency of the simulation + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types. +The first coordinate is assumed to be r, the second is y. + |
+ + required + | +
+ epsilon
+ |
+
+ vfdslice
+ |
+
+
+
+ Dielectric constant + |
+ + required + | +
+ rmin
+ |
+
+ float
+ |
+
+
+
+ Radius of curvature for the simulation. This should be the minimum value of +r within the simulation domain. + |
+ + required + | +
Returns:
+| Name | Type | +Description | +
|---|---|---|
e_xys |
+ NDArray[complex128]
+ |
+
+
+
+ NDArray of vfdfield_t specifying fields. First dimension is mode number. + |
+
angular_wavenumbers |
+ NDArray[complex128]
+ |
+
+
+
+ list of wavenumbers in 1/rad units. + |
+
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 |
+
+ ()
+ |
+
+ **kwargs
+ |
+
+ Any
+ |
+
+
+
+ passed to |
+
+ {}
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ tuple[vcfdfield2, complex]
+ |
+
+
+
+ (e_xy, angular_wavenumber) + |
+
linear_wavenumbers(
+ e_xys: Sequence[vcfdfield2] | NDArray[complex128],
+ angular_wavenumbers: ArrayLike,
+ epsilon: vfdslice,
+ dxes: dx_lists2_t,
+ rmin: float,
+) -> NDArray[numpy.complex128]
+Calculate linear wavenumbers (1/distance) based on angular wavenumbers (1/rad) + and the mode's energy distribution.
+ + +Parameters:
+| Name | +Type | +Description | +Default | +
|---|---|---|---|
+ e_xys
+ |
+
+ Sequence[vcfdfield2] | NDArray[complex128]
+ |
+
+
+
+ Vectorized mode fields with shape (num_modes, 2 * x *y) + |
+ + required + | +
+ angular_wavenumbers
+ |
+
+ ArrayLike
+ |
+
+
+
+ Wavenumbers assuming fields have theta-dependence of
+ |
+ + required + | +
+ epsilon
+ |
+
+ vfdslice
+ |
+
+
+
+ Vectorized dielectric constant grid with shape (3, x, y) + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + 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(
+ 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
+ |
+ + required + | +
+ omega
+ |
+
+ float
+ |
+
+
+
+ The angular frequency of the system + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + 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(
+ 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
+ |
+ + required + | +
+ omega
+ |
+
+ float
+ |
+
+
+
+ The angular frequency of the system + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + 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(
+ 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:
+Parameters:
+| Name | +Type | +Description | +Default | +
|---|---|---|---|
+ angular_wavenumber
+ |
+
+ complex
+ |
+
+
+
+ Wavenumber assuming fields have theta-dependence of
+ |
+ + required + | +
+ omega
+ |
+
+ float
+ |
+
+
+
+ The angular frequency of the system + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + 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(
+ 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 |
+ + 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 |
+
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
+ |
+ + required + | +
+ omega
+ |
+
+ float
+ |
+
+
+
+ The angular frequency of the system + |
+ + required + | +
+ dxes
+ |
+
+ dx_lists2_t
+ |
+
+
+
+ Grid parameters |
+ + 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 |
+
+ 0
+ |
+
Returns:
+| Type | +Description | +
|---|---|
+ vcfdslice_t
+ |
+
+
+
+
|
+
+ vcfdslice_t
+ |
+
+
+
+ and contains all three vector components. + |
+
The normalization step is delegated to _normalized_fields(...), which
+enforces unit forward power under the discrete inner product
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.