diff --git a/make_docs.sh b/make_docs.sh
new file mode 100755
index 0000000..f2f2fe2
--- /dev/null
+++ b/make_docs.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+cd ~/projects/meanas
+pdoc3 --html --force --template-dir pdoc_templates -o doc .
diff --git a/meanas/eigensolvers.py b/meanas/eigensolvers.py
index 5ca9962..960df60 100644
--- a/meanas/eigensolvers.py
+++ b/meanas/eigensolvers.py
@@ -15,10 +15,13 @@ def power_iteration(operator: sparse.spmatrix,
"""
Use power iteration to estimate the dominant eigenvector of a matrix.
- :param operator: Matrix to analyze.
- :param guess_vector: Starting point for the eigenvector. Default is a randomly chosen vector.
- :param iterations: Number of iterations to perform. Default 20.
- :return: (Largest-magnitude eigenvalue, Corresponding eigenvector estimate)
+ Args:
+ operator: Matrix to analyze.
+ guess_vector: Starting point for the eigenvector. Default is a randomly chosen vector.
+ iterations: Number of iterations to perform. Default 20.
+
+ Returns:
+ (Largest-magnitude eigenvalue, Corresponding eigenvector estimate)
"""
if numpy.any(numpy.equal(guess_vector, None)):
v = numpy.random.rand(operator.shape[0])
@@ -91,12 +94,15 @@ def signed_eigensolve(operator: sparse.spmatrix or spalg.LinearOperator,
Find the largest-magnitude positive-only (or negative-only) eigenvalues and
eigenvectors of the provided matrix.
- :param operator: Matrix to analyze.
- :param how_many: How many eigenvalues to find.
- :param negative: Whether to find negative-only eigenvalues.
- Default False (positive only).
- :return: (sorted list of eigenvalues, 2D ndarray of corresponding eigenvectors)
- eigenvectors[:, k] corresponds to the k-th eigenvalue
+ Args:
+ operator: Matrix to analyze.
+ how_many: How many eigenvalues to find.
+ negative: Whether to find negative-only eigenvalues.
+ Default False (positive only).
+
+ Returns:
+ (sorted list of eigenvalues, 2D ndarray of corresponding eigenvectors)
+ `eigenvectors[:, k]` corresponds to the k-th eigenvalue
"""
# Use power iteration to estimate the dominant eigenvector
lm_eigval, _ = power_iteration(operator)
diff --git a/meanas/fdfd/__init__.py b/meanas/fdfd/__init__.py
index 80ba55e..d0d58a9 100644
--- a/meanas/fdfd/__init__.py
+++ b/meanas/fdfd/__init__.py
@@ -1,2 +1,17 @@
-from . import solvers, operators, functional, scpml, waveguide, waveguide_mode
+"""
+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.
+"""
+from . import solvers, operators, functional, scpml, waveguide_2d, waveguide_3d
# from . import farfield, bloch TODO
diff --git a/meanas/fdfd/bloch.py b/meanas/fdfd/bloch.py
index 9c252e7..0ca64a7 100644
--- a/meanas/fdfd/bloch.py
+++ b/meanas/fdfd/bloch.py
@@ -83,7 +83,7 @@ import scipy.optimize
from scipy.linalg import norm
import scipy.sparse.linalg as spalg
-from . import field_t
+from .. import field_t
logger = logging.getLogger(__name__)
diff --git a/meanas/fdfd/farfield.py b/meanas/fdfd/farfield.py
index faa25b5..665e70f 100644
--- a/meanas/fdfd/farfield.py
+++ b/meanas/fdfd/farfield.py
@@ -1,7 +1,7 @@
"""
Functions for performing near-to-farfield transformation (and the reverse).
"""
-from typing import Dict, List
+from typing import Dict, List, Any
import numpy
from numpy.fft import fft2, fftshift, fftfreq, ifft2, ifftshift
from numpy import pi
@@ -14,7 +14,7 @@ def near_to_farfield(E_near: field_t,
dx: float,
dy: float,
padded_size: List[int] = None
- ) -> Dict[str]:
+ ) -> Dict[str, Any]:
"""
Compute the farfield, i.e. the distribution of the fields after propagation
through several wavelengths of uniform medium.
@@ -122,7 +122,7 @@ def far_to_nearfield(E_far: field_t,
dkx: float,
dky: float,
padded_size: List[int] = None
- ) -> Dict[str]:
+ ) -> Dict[str, Any]:
"""
Compute the farfield, i.e. the distribution of the fields after propagation
through several wavelengths of uniform medium.
diff --git a/meanas/fdfd/functional.py b/meanas/fdfd/functional.py
index bd31192..74b66da 100644
--- a/meanas/fdfd/functional.py
+++ b/meanas/fdfd/functional.py
@@ -2,88 +2,41 @@
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 field inputs with shape (3, X, Y, Z),
+The functions generated here expect `field_t` inputs with shape (3, X, Y, Z),
e.g. E = [E_x, E_y, E_z] where each component has shape (X, Y, Z)
"""
-from typing import List, Callable
+from typing import List, Callable, Tuple
import numpy
from .. import dx_lists_t, field_t
+from ..fdmath.functional import curl_forward, curl_back
__author__ = 'Jan Petykiewicz'
-functional_matrix = Callable[[field_t], field_t]
-
-
-def curl_h(dxes: dx_lists_t) -> functional_matrix:
- """
- Curl operator for use with the H field.
-
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :return: Function for taking the discretized curl of the H-field, F(H) -> curlH
- """
- dxyz_b = numpy.meshgrid(*dxes[1], indexing='ij')
-
- def dh(f, ax):
- return (f - numpy.roll(f, 1, axis=ax)) / dxyz_b[ax]
-
- def ch_fun(h: field_t) -> field_t:
- e = numpy.empty_like(h)
- e[0] = dh(h[2], 1)
- e[0] -= dh(h[1], 2)
- e[1] = dh(h[0], 2)
- e[1] -= dh(h[2], 0)
- e[2] = dh(h[1], 0)
- e[2] -= dh(h[0], 1)
- return e
-
- return ch_fun
-
-
-def curl_e(dxes: dx_lists_t) -> functional_matrix:
- """
- Curl operator for use with the E field.
-
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :return: Function for taking the discretized curl of the E-field, F(E) -> curlE
- """
- dxyz_a = numpy.meshgrid(*dxes[0], indexing='ij')
-
- def de(f, ax):
- return (numpy.roll(f, -1, axis=ax) - f) / dxyz_a[ax]
-
- def ce_fun(e: field_t) -> field_t:
- h = numpy.empty_like(e)
- h[0] = de(e[2], 1)
- h[0] -= de(e[1], 2)
- h[1] = de(e[0], 2)
- h[1] -= de(e[2], 0)
- h[2] = de(e[1], 0)
- h[2] -= de(e[0], 1)
- return h
-
- return ce_fun
+field_transform_t = Callable[[field_t], field_t]
def e_full(omega: complex,
dxes: dx_lists_t,
epsilon: field_t,
mu: field_t = None
- ) -> functional_matrix:
+ ) -> field_transform_t:
"""
- Wave operator del x (1/mu * del x) - omega**2 * epsilon, for use with E-field,
- with wave equation
- (del x (1/mu * del x) - omega**2 * epsilon) E = -i * omega * J
+ Wave operator for use with E-field. See `operators.e_full` for details.
- :param omega: Angular frequency of the simulation
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :param epsilon: Dielectric constant
- :param mu: Magnetic permeability (default 1 everywhere)
- :return: Function implementing the wave operator A(E) -> E
+ Args:
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
+ epsilon: Dielectric constant
+ mu: Magnetic permeability (default 1 everywhere)
+
+ Return:
+ Function `f` implementing the wave operator
+ `f(E)` -> `-i * omega * J`
"""
- ch = curl_h(dxes)
- ce = curl_e(dxes)
+ ch = curl_back(dxes[1])
+ ce = curl_forward(dxes[0])
def op_1(e):
curls = ch(ce(e))
@@ -103,18 +56,23 @@ def eh_full(omega: complex,
dxes: dx_lists_t,
epsilon: field_t,
mu: field_t = None
- ) -> functional_matrix:
+ ) -> Callable[[field_t, field_t], Tuple[field_t, field_t]]:
"""
Wave operator for full (both E and H) field representation.
+ See `operators.eh_full`.
- :param omega: Angular frequency of the simulation
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :param epsilon: Dielectric constant
- :param mu: Magnetic permeability (default 1 everywhere)
- :return: Function implementing the wave operator A(E, H) -> (E, H)
+ Args:
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
+ epsilon: Dielectric constant
+ mu: Magnetic permeability (default 1 everywhere)
+
+ Returns:
+ Function `f` implementing the wave operator
+ `f(E, H)` -> `(J, -M)`
"""
- ch = curl_h(dxes)
- ce = curl_e(dxes)
+ ch = curl_back(dxes[1])
+ ce = curl_forward(dxes[0])
def op_1(e, h):
return (ch(h) - 1j * omega * epsilon * e,
@@ -133,23 +91,27 @@ def eh_full(omega: complex,
def e2h(omega: complex,
dxes: dx_lists_t,
mu: field_t = None,
- ) -> functional_matrix:
+ ) -> field_transform_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.
+ Utility operator for converting the `E` field into the `H` field.
+ For use with `e_full` -- assumes that there is no magnetic current `M`.
- :param omega: Angular frequency of the simulation
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :param mu: Magnetic permeability (default 1 everywhere)
- :return: Function for converting E to H
- """
- A2 = curl_e(dxes)
+ Args:
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+ mu: Magnetic permeability (default 1 everywhere)
+
+ Return:
+ Function `f` for converting `E` to `H`,
+ `f(E)` -> `H`
+ """
+ ce = curl_forward(dxes[0])
def e2h_1_1(e):
- return A2(e) / (-1j * omega)
+ return ce(e) / (-1j * omega)
def e2h_mu(e):
- return A2(e) / (-1j * omega * mu)
+ return ce(e) / (-1j * omega * mu)
if numpy.any(numpy.equal(mu, None)):
return e2h_1_1
@@ -160,18 +122,22 @@ def e2h(omega: complex,
def m2j(omega: complex,
dxes: dx_lists_t,
mu: field_t = None,
- ) -> functional_matrix:
+ ) -> field_transform_t:
"""
- Utility operator for converting magnetic current (M) distribution
- into equivalent electric current distribution (J).
- For use with e.g. e_full().
+ Utility operator for converting magnetic current `M` distribution
+ into equivalent electric current distribution `J`.
+ For use with e.g. `e_full`.
- :param omega: Angular frequency of the simulation
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :param mu: Magnetic permeability (default 1 everywhere)
- :return: Function for converting M to J
- """
- ch = curl_h(dxes)
+ Args:
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+ mu: Magnetic permeability (default 1 everywhere)
+
+ Returns:
+ Function `f` for converting `M` to `J`,
+ `f(M)` -> `J`
+ """
+ ch = curl_back(dxes[1])
def m2j_mu(m):
J = ch(m / mu) / (-1j * omega)
@@ -192,10 +158,23 @@ def e_tfsf_source(TF_region: field_t,
dxes: dx_lists_t,
epsilon: field_t,
mu: field_t = None,
- ) -> functional_matrix:
+ ) -> field_transform_t:
"""
- Operator that turuns an E-field distribution into a total-field/scattered-field
+ Operator that turns an E-field distribution into a total-field/scattered-field
(TFSF) source.
+
+ Args:
+ TF_region: mask which is set to 1 in the total-field region, and 0 elsewhere
+ (i.e. in the scattered-field region).
+ Should have the same shape as the simulation grid, e.g. `epsilon[0].shape`.
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+ epsilon: Dielectric constant distribution
+ mu: Magnetic permeability (default 1 everywhere)
+
+ Returns:
+ Function `f` which takes an E field and returns a current distribution,
+ `f(E)` -> `J`
"""
# TODO documentation
A = e_full(omega, dxes, epsilon, mu)
@@ -205,7 +184,28 @@ def e_tfsf_source(TF_region: field_t,
return neg_iwj / (-1j * omega)
-def poynting_e_cross_h(dxes: dx_lists_t):
+def poynting_e_cross_h(dxes: dx_lists_t) -> Callable[[field_t, field_t], field_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 \\)
+
+ Note:
+ This function also shifts the input `E` field by one cell as required
+ for computing the Poynting cross product (see `meanas.fdfd` module docs).
+
+ Note:
+ If `E` and `H` are peak amplitudes as assumed elsewhere in this code,
+ the time-average of the poynting vector is ` = Re(S)/2 = Re(E x H) / 2`.
+ The factor of `1/2` can be omitted if root-mean-square quantities are used
+ instead.
+
+ Args:
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+
+ Returns:
+ Function `f` that returns E x H as required for the poynting vector.
+ """
def exh(e: field_t, h: field_t):
s = numpy.empty_like(e)
ex = e[0] * dxes[0][0][:, None, None]
diff --git a/meanas/fdfd/operators.py b/meanas/fdfd/operators.py
index 905c23e..20cf96a 100644
--- a/meanas/fdfd/operators.py
+++ b/meanas/fdfd/operators.py
@@ -1,18 +1,19 @@
"""
Sparse matrix operators for use with electromagnetic wave equations.
-These functions return sparse-matrix (scipy.sparse.spmatrix) representations of
+These functions return sparse-matrix (`scipy.sparse.spmatrix`) representations of
a variety of operators, intended for use with E and H fields vectorized using the
- meanas.vec() and .unvec() functions (column-major/Fortran ordering).
+ `meanas.vec()` and `meanas.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).
+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 meanas.dx_lists_type; see
-the meanas.types submodule for details.
+Many of these functions require a `dxes` parameter, of type `dx_lists_t`; see
+the `meanas.types` submodule for details.
The following operators are included:
+
- E-only wave operator
- H-only wave operator
- EH wave operator
@@ -20,8 +21,6 @@ The following operators are included:
- E to H conversion
- M to J conversion
- Poynting cross products
-
-Also available:
- Circular shifts
- Discrete derivatives
- Averaging operators
@@ -33,6 +32,7 @@ import numpy
import scipy.sparse as sparse
from .. import vec, dx_lists_t, vfield_t
+from ..fdmath.operators import shift_with_mirror, rotation, curl_forward, curl_back
__author__ = 'Jan Petykiewicz'
@@ -46,26 +46,35 @@ def e_full(omega: complex,
pmc: vfield_t = None,
) -> sparse.spmatrix:
"""
- Wave operator del x (1/mu * del x) - omega**2 * epsilon, for use with E-field,
- with wave equation
- (del x (1/mu * del x) - omega**2 * epsilon) E = -i * omega * J
+ Wave operator
+ $$ \\nabla \\times (\\frac{1}{\\mu} \\nabla \\times) - \\omega^2 \\epsilon $$
- To make this matrix symmetric, use the preconditions from e_full_preconditioners().
+ del x (1/mu * del x) - omega**2 * epsilon
- :param omega: Angular frequency of the simulation
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :param epsilon: Vectorized dielectric constant
- :param mu: Vectorized magnetic permeability (default 1 everywhere).
- :param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted
- as containing a perfect electrical conductor (PEC).
- The PEC is applied per-field-component (ie, pec.size == epsilon.size)
- :param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted
- as containing a perfect magnetic conductor (PMC).
- The PMC is applied per-field-component (ie, pmc.size == epsilon.size)
- :return: Sparse matrix containing the wave operator
+ for use with the E-field, with wave equation
+ $$ (\\nabla \\times (\\frac{1}{\\mu} \\nabla \\times) - \\omega^2 \\epsilon) E = -\\imath \\omega J $$
+
+ (del x (1/mu * del x) - omega**2 * epsilon) E = -i * omega * J
+
+ To make this matrix symmetric, use the preconditioners from `e_full_preconditioners()`.
+
+ Args:
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+ epsilon: Vectorized dielectric constant
+ mu: Vectorized magnetic permeability (default 1 everywhere).
+ pec: Vectorized mask specifying PEC cells. Any cells where `pec != 0` are interpreted
+ as containing a perfect electrical conductor (PEC).
+ The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`)
+ pmc: Vectorized mask specifying PMC cells. Any cells where `pmc != 0` are interpreted
+ as containing a perfect magnetic conductor (PMC).
+ The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`)
+
+ Returns:
+ Sparse matrix containing the wave operator.
"""
- ce = curl_e(dxes)
- ch = curl_h(dxes)
+ ch = curl_back(dxes[1])
+ ce = curl_forward(dxes[0])
if numpy.any(numpy.equal(pec, None)):
pe = sparse.eye(epsilon.size)
@@ -90,15 +99,18 @@ def e_full(omega: complex,
def e_full_preconditioners(dxes: dx_lists_t
) -> Tuple[sparse.spmatrix, sparse.spmatrix]:
"""
- Left and right preconditioners (Pl, Pr) for symmetrizing the e_full 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
+ 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
+ The preconditioner matrices are diagonal and complex, with `Pr = 1 / Pl`
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :return: Preconditioner matrices (Pl, Pr)
+ Args:
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+
+ Returns:
+ Preconditioner matrices `(Pl, Pr)`.
"""
p_squared = [dxes[0][0][:, None, None] * dxes[1][1][None, :, None] * dxes[1][2][None, None, :],
dxes[1][0][:, None, None] * dxes[0][1][None, :, None] * dxes[1][2][None, None, :],
@@ -118,24 +130,33 @@ def h_full(omega: complex,
pmc: vfield_t = None,
) -> sparse.spmatrix:
"""
- Wave operator del x (1/epsilon * del x) - omega**2 * mu, for use with H-field,
- with wave equation
- (del x (1/epsilon * del x) - omega**2 * mu) H = i * omega * M
+ Wave operator
+ $$ \\nabla \\times (\\frac{1}{\\epsilon} \\nabla \\times) - \\omega^2 \\mu $$
- :param omega: Angular frequency of the simulation
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :param epsilon: Vectorized dielectric constant
- :param mu: Vectorized magnetic permeability (default 1 everywhere)
- :param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted
- as containing a perfect electrical conductor (PEC).
- The PEC is applied per-field-component (ie, pec.size == epsilon.size)
- :param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted
- as containing a perfect magnetic conductor (PMC).
- The PMC is applied per-field-component (ie, pmc.size == epsilon.size)
- :return: Sparse matrix containing the wave operator
+ del x (1/epsilon * del x) - omega**2 * mu
+
+ for use with the H-field, with wave equation
+ $$ (\\nabla \\times (\\frac{1}{\\epsilon} \\nabla \\times) - \\omega^2 \\mu) E = \\imath \\omega M $$
+
+ (del x (1/epsilon * del x) - omega**2 * mu) E = i * omega * M
+
+ Args:
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+ epsilon: Vectorized dielectric constant
+ mu: Vectorized magnetic permeability (default 1 everywhere)
+ pec: Vectorized mask specifying PEC cells. Any cells where `pec != 0` are interpreted
+ as containing a perfect electrical conductor (PEC).
+ The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`)
+ pmc: Vectorized mask specifying PMC cells. Any cells where `pmc != 0` are interpreted
+ as containing a perfect magnetic conductor (PMC).
+ The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`)
+
+ Returns:
+ Sparse matrix containing the wave operator.
"""
- ec = curl_e(dxes)
- hc = curl_h(dxes)
+ ch = curl_back(dxes[1])
+ ce = curl_forward(dxes[0])
if numpy.any(numpy.equal(pec, None)):
pe = sparse.eye(epsilon.size)
@@ -153,7 +174,7 @@ def h_full(omega: complex,
else:
m = sparse.diags(mu)
- A = pm @ (ec @ pe @ e_div @ hc - omega**2 * m) @ pm
+ A = pm @ (ce @ pe @ e_div @ ch - omega**2 * m) @ pm
return A
@@ -165,24 +186,42 @@ def eh_full(omega: complex,
pmc: vfield_t = None
) -> sparse.spmatrix:
"""
- Wave operator for [E, H] field representation. This operator implements Maxwell's
+ 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]]
+ $$ \\begin{bmatrix}
+ -\\imath \\omega \\epsilon & \\nabla \\times \\\\
+ \\nabla \\times & \\imath \\omega \\mu
+ \\end{bmatrix} $$
- for use with a field vector of the form hstack(vec(E), vec(H)).
+ [[-i * omega * epsilon, del x ],
+ [del x, i * omega * mu]]
- :param omega: Angular frequency of the simulation
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :param epsilon: Vectorized dielectric constant
- :param mu: Vectorized magnetic permeability (default 1 everywhere)
- :param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted
- as containing a perfect electrical conductor (PEC).
- The PEC is applied per-field-component (i.e., pec.size == epsilon.size)
- :param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted
- as containing a perfect magnetic conductor (PMC).
- The PMC is applied per-field-component (i.e., pmc.size == epsilon.size)
- :return: Sparse matrix containing the wave operator
+ for use with a field vector of the form `cat(vec(E), vec(H))`:
+ $$ \\begin{bmatrix}
+ -\\imath \\omega \\epsilon & \\nabla \\times \\\\
+ \\nabla \\times & \\imath \\omega \\mu
+ \\end{bmatrix}
+ \\begin{bmatrix} E \\\\
+ H
+ \\end{bmatrix}
+ = \\begin{bmatrix} J \\\\
+ -M
+ \\end{bmatrix} $$
+
+ Args:
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+ epsilon: Vectorized dielectric constant
+ mu: Vectorized magnetic permeability (default 1 everywhere)
+ pec: Vectorized mask specifying PEC cells. Any cells where `pec != 0` are interpreted
+ as containing a perfect electrical conductor (PEC).
+ The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`)
+ pmc: Vectorized mask specifying PMC cells. Any cells where `pmc != 0` are interpreted
+ as containing a perfect magnetic conductor (PMC).
+ The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`)
+
+ Returns:
+ Sparse matrix containing the wave operator.
"""
if numpy.any(numpy.equal(pec, None)):
pe = sparse.eye(epsilon.size)
@@ -200,34 +239,14 @@ def eh_full(omega: complex,
iwm *= sparse.diags(mu)
iwm = pm @ iwm @ pm
- A1 = pe @ curl_h(dxes) @ pm
- A2 = pm @ curl_e(dxes) @ pe
+ A1 = pe @ curl_back(dxes[1]) @ pm
+ A2 = pm @ curl_forward(dxes[0]) @ pe
A = sparse.bmat([[-iwe, A1],
[A2, iwm]])
return A
-def curl_h(dxes: dx_lists_t) -> sparse.spmatrix:
- """
- Curl operator for use with the H field.
-
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :return: Sparse matrix for taking the discretized curl of the H-field
- """
- return cross(deriv_back(dxes[1]))
-
-
-def curl_e(dxes: dx_lists_t) -> sparse.spmatrix:
- """
- Curl operator for use with the E field.
-
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :return: Sparse matrix for taking the discretized curl of the E-field
- """
- return cross(deriv_forward(dxes[0]))
-
-
def e2h(omega: complex,
dxes: dx_lists_t,
mu: vfield_t = None,
@@ -235,17 +254,20 @@ def e2h(omega: complex,
) -> sparse.spmatrix:
"""
Utility operator for converting the E field into the H field.
- For use with e_full -- assumes that there is no magnetic current M.
+ For use with `e_full()` -- assumes that there is no magnetic current M.
- :param omega: Angular frequency of the simulation
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :param mu: Vectorized magnetic permeability (default 1 everywhere)
- :param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted
- as containing a perfect magnetic conductor (PMC).
- The PMC is applied per-field-component (ie, pmc.size == epsilon.size)
- :return: Sparse matrix for converting E to H
+ Args:
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+ mu: Vectorized magnetic permeability (default 1 everywhere)
+ pmc: Vectorized mask specifying PMC cells. Any cells where `pmc != 0` are interpreted
+ as containing a perfect magnetic conductor (PMC).
+ The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`)
+
+ Returns:
+ Sparse matrix for converting E to H.
"""
- op = curl_e(dxes) / (-1j * omega)
+ op = curl_forward(dxes[0]) / (-1j * omega)
if not numpy.any(numpy.equal(mu, None)):
op = sparse.diags(1 / mu) @ op
@@ -261,16 +283,18 @@ def m2j(omega: complex,
mu: vfield_t = None
) -> sparse.spmatrix:
"""
- Utility operator for converting M field into J.
- Converts a magnetic current M into an electric current J.
- For use with eg. e_full.
+ Operator for converting a magnetic current M into an electric current J.
+ For use with eg. `e_full()`.
- :param omega: Angular frequency of the simulation
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :param mu: Vectorized magnetic permeability (default 1 everywhere)
- :return: Sparse matrix for converting E to H
+ Args:
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+ mu: Vectorized magnetic permeability (default 1 everywhere)
+
+ Returns:
+ Sparse matrix for converting M to J.
"""
- op = curl_h(dxes) / (1j * omega)
+ op = curl_back(dxes[1]) / (1j * omega)
if not numpy.any(numpy.equal(mu, None)):
op = op @ sparse.diags(1 / mu)
@@ -278,178 +302,17 @@ def m2j(omega: complex,
return op
-def rotation(axis: int, shape: List[int], shift_distance: int=1) -> sparse.spmatrix:
- """
- Utility operator for performing a circular shift along a specified axis by a
- specified number of elements.
-
- :param axis: Axis to shift along. x=0, y=1, z=2
- :param shape: Shape of the grid being shifted
- :param shift_distance: Number of cells to shift by. May be negative. Default 1.
- :return: Sparse matrix for performing the circular shift
- """
- if len(shape) not in (2, 3):
- raise Exception('Invalid shape: {}'.format(shape))
- if axis not in range(len(shape)):
- raise Exception('Invalid direction: {}, shape is {}'.format(axis, shape))
-
- shifts = [abs(shift_distance) if a == axis else 0 for a in range(3)]
- shifted_diags = [(numpy.arange(n) + s) % n for n, s in zip(shape, shifts)]
- ijk = numpy.meshgrid(*shifted_diags, indexing='ij')
-
- n = numpy.prod(shape)
- i_ind = numpy.arange(n)
- j_ind = numpy.ravel_multi_index(ijk, shape, order='C')
-
- vij = (numpy.ones(n), (i_ind, j_ind.ravel(order='C')))
-
- d = sparse.csr_matrix(vij, shape=(n, n))
-
- if shift_distance < 0:
- d = d.T
-
- return d
-
-
-def shift_with_mirror(axis: int, shape: List[int], shift_distance: int=1) -> sparse.spmatrix:
- """
- Utility operator for performing an n-element shift along a specified axis, with mirror
- boundary conditions applied to the cells beyond the receding edge.
-
- :param axis: Axis to shift along. x=0, y=1, z=2
- :param shape: Shape of the grid being shifted
- :param shift_distance: Number of cells to shift by. May be negative. Default 1.
- :return: Sparse matrix for performing the circular shift
- """
- if len(shape) not in (2, 3):
- raise Exception('Invalid shape: {}'.format(shape))
- if axis not in range(len(shape)):
- raise Exception('Invalid direction: {}, shape is {}'.format(axis, shape))
- if shift_distance >= shape[axis]:
- raise Exception('Shift ({}) is too large for axis {} of size {}'.format(
- shift_distance, axis, shape[axis]))
-
- def mirrored_range(n, s):
- v = numpy.arange(n) + s
- v = numpy.where(v >= n, 2 * n - v - 1, v)
- v = numpy.where(v < 0, - 1 - v, v)
- return v
-
- shifts = [shift_distance if a == axis else 0 for a in range(3)]
- shifted_diags = [mirrored_range(n, s) for n, s in zip(shape, shifts)]
- ijk = numpy.meshgrid(*shifted_diags, indexing='ij')
-
- n = numpy.prod(shape)
- i_ind = numpy.arange(n)
- j_ind = numpy.ravel_multi_index(ijk, shape, order='C')
-
- vij = (numpy.ones(n), (i_ind, j_ind.ravel(order='C')))
-
- d = sparse.csr_matrix(vij, shape=(n, n))
- return d
-
-
-def deriv_forward(dx_e: List[numpy.ndarray]) -> List[sparse.spmatrix]:
- """
- Utility operators for taking discretized derivatives (forward variant).
-
- :param dx_e: Lists of cell sizes for all axes [[dx_0, dx_1, ...], ...].
- :return: List of operators for taking forward derivatives along each axis.
- """
- shape = [s.size for s in dx_e]
- n = numpy.prod(shape)
-
- dx_e_expanded = numpy.meshgrid(*dx_e, indexing='ij')
-
- def deriv(axis):
- return rotation(axis, shape, 1) - sparse.eye(n)
-
- Ds = [sparse.diags(+1 / dx.ravel(order='C')) @ deriv(a)
- for a, dx in enumerate(dx_e_expanded)]
-
- return Ds
-
-
-def deriv_back(dx_h: List[numpy.ndarray]) -> List[sparse.spmatrix]:
- """
- Utility operators for taking discretized derivatives (backward variant).
-
- :param dx_h: Lists of cell sizes for all axes [[dx_0, dx_1, ...], ...].
- :return: List of operators for taking forward derivatives along each axis.
- """
- shape = [s.size for s in dx_h]
- n = numpy.prod(shape)
-
- dx_h_expanded = numpy.meshgrid(*dx_h, indexing='ij')
-
- def deriv(axis):
- return rotation(axis, shape, -1) - sparse.eye(n)
-
- Ds = [sparse.diags(-1 / dx.ravel(order='C')) @ deriv(a)
- for a, dx in enumerate(dx_h_expanded)]
-
- return Ds
-
-
-def cross(B: List[sparse.spmatrix]) -> sparse.spmatrix:
- """
- Cross product operator
-
- :param B: List [Bx, By, Bz] of sparse matrices corresponding to the x, y, z
- portions of the operator on the left side of the cross product.
- :return: Sparse matrix corresponding to (B x), where x is the cross product
- """
- n = B[0].shape[0]
- zero = sparse.csr_matrix((n, n))
- return sparse.bmat([[zero, -B[2], B[1]],
- [B[2], zero, -B[0]],
- [-B[1], B[0], zero]])
-
-
-def vec_cross(b: vfield_t) -> sparse.spmatrix:
- """
- Vector cross product operator
-
- :param b: Vector on the left side of the cross product
- :return: Sparse matrix corresponding to (b x), where x is the cross product
- """
- B = [sparse.diags(c) for c in numpy.split(b, 3)]
- return cross(B)
-
-
-def avgf(axis: int, shape: List[int]) -> sparse.spmatrix:
- """
- Forward average operator (x4 = (x4 + x5) / 2)
-
- :param axis: Axis to average along (x=0, y=1, z=2)
- :param shape: Shape of the grid to average
- :return: Sparse matrix for forward average operation
- """
- if len(shape) not in (2, 3):
- raise Exception('Invalid shape: {}'.format(shape))
-
- n = numpy.prod(shape)
- return 0.5 * (sparse.eye(n) + rotation(axis, shape))
-
-
-def avgb(axis: int, shape: List[int]) -> sparse.spmatrix:
- """
- Backward average operator (x4 = (x4 + x3) / 2)
-
- :param axis: Axis to average along (x=0, y=1, z=2)
- :param shape: Shape of the grid to average
- :return: Sparse matrix for backward average operation
- """
- return avgf(axis, shape).T
-
-
def poynting_e_cross(e: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix:
"""
- Operator for computing the Poynting vector, containing the (E x) portion of the Poynting vector.
+ Operator for computing the Poynting vector, containing the
+ (E x) portion of the Poynting vector.
- :param e: Vectorized E-field for the ExH cross product
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :return: Sparse matrix containing (E x) portion of Poynting cross product
+ Args:
+ e: Vectorized E-field for the ExH cross product
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+
+ Returns:
+ Sparse matrix containing (E x) portion of Poynting cross product.
"""
shape = [len(dx) for dx in dxes[0]]
@@ -472,9 +335,12 @@ def poynting_h_cross(h: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix:
"""
Operator for computing the Poynting vector, containing the (H x) portion of the Poynting vector.
- :param h: Vectorized H-field for the HxE cross product
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :return: Sparse matrix containing (H x) portion of Poynting cross product
+ Args:
+ h: Vectorized H-field for the HxE cross product
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+
+ Returns:
+ Sparse matrix containing (H x) portion of Poynting cross product.
"""
shape = [len(dx) for dx in dxes[0]]
@@ -499,8 +365,22 @@ def e_tfsf_source(TF_region: vfield_t,
mu: vfield_t = None,
) -> sparse.spmatrix:
"""
- Operator that turns an E-field distribution into a total-field/scattered-field
- (TFSF) source.
+ Operator that turns a desired E-field distribution into a
+ total-field/scattered-field (TFSF) source.
+
+ TODO: Reference Rumpf paper
+
+ Args:
+ TF_region: Mask, which is set to 1 inside the total-field region and 0 in the
+ scattered-field region
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+ epsilon: Vectorized dielectric constant
+ mu: Vectorized magnetic permeability (default 1 everywhere).
+
+ Returns:
+ Sparse matrix that turns an E-field into a current (J) distribution.
+
"""
# TODO documentation
A = e_full(omega, dxes, epsilon, mu)
@@ -518,7 +398,19 @@ def e_boundary_source(mask: vfield_t,
"""
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.
+ `e_tfsf_source()` with an additional masking step.
+
+ Args:
+ mask: 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.
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+ epsilon: Vectorized dielectric constant
+ mu: Vectorized magnetic permeability (default 1 everywhere).
+
+ Returns:
+ Sparse matrix that turns an E-field into a current (J) distribution.
"""
full = e_tfsf_source(TF_region=mask, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu)
diff --git a/meanas/fdfd/scpml.py b/meanas/fdfd/scpml.py
index 897d43a..c9a93e6 100644
--- a/meanas/fdfd/scpml.py
+++ b/meanas/fdfd/scpml.py
@@ -1,5 +1,5 @@
"""
-Functions for creating stretched coordinate PMLs.
+Functions for creating stretched coordinate perfectly matched layer (PML) absorbers.
"""
from typing import List, Callable
@@ -7,24 +7,29 @@ import numpy
from .. import dx_lists_t
+
__author__ = 'Jan Petykiewicz'
-s_function_type = Callable[[float], float]
+s_function_t = Callable[[float], float]
+"""Typedef for s-functions"""
def prepare_s_function(ln_R: float = -16,
m: float = 4
- ) -> s_function_type:
+ ) -> s_function_t:
"""
Create an s_function to pass to the SCPML functions. This is used when you would like to
customize the PML parameters.
- :param ln_R: Natural logarithm of the desired reflectance
- :param m: Polynomial order for the PML (imaginary part increases as distance ** m)
- :return: An s_function, which takes an ndarray (distances) and returns an ndarray (complex part
- of the cell width; needs to be divided by sqrt(epilon_effective) * real(omega))
- before use.
+ Args:
+ ln_R: Natural logarithm of the desired reflectance
+ m: Polynomial order for the PML (imaginary part increases as distance ** m)
+
+ Returns:
+ An s_function, which takes an ndarray (distances) and returns an ndarray (complex part
+ of the cell width; needs to be divided by `sqrt(epilon_effective) * real(omega))`
+ before use.
"""
def s_factor(distance: numpy.ndarray) -> numpy.ndarray:
s_max = (m + 1) * ln_R / 2 # / 2 because we assume periodic boundaries
@@ -36,26 +41,29 @@ def uniform_grid_scpml(shape: numpy.ndarray or List[int],
thicknesses: numpy.ndarray or List[int],
omega: float,
epsilon_effective: float = 1.0,
- s_function: s_function_type = None,
+ s_function: s_function_t = None,
) -> dx_lists_t:
"""
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(...).
+ If you want something more fine-grained, check out `stretch_with_scpml(...)`.
- :param shape: Shape of the grid, including the PMLs (which are 2*thicknesses thick)
- :param thicknesses: [th_x, th_y, th_z] Thickness of the PML in each direction.
- Both polarities are added.
- Each th_ of pml is applied twice, once on each edge of the grid along the given axis.
- th_* may be zero, in which case no pml is added.
- :param omega: Angular frequency for the simulation
- :param epsilon_effective: Effective epsilon of the PML. Match this to the material
- at the edge of your grid.
- Default 1.
- :param s_function: s_function created by prepare_s_function(...), allowing
- customization of pml parameters.
- Default uses prepare_s_function() with no parameters.
- :return: Complex cell widths (dx_lists)
+ Args:
+ shape: Shape of the grid, including the PMLs (which are 2*thicknesses thick)
+ thicknesses: `[th_x, th_y, th_z]`
+ Thickness of the PML in each direction.
+ Both polarities are added.
+ Each th_ of pml is applied twice, once on each edge of the grid along the given axis.
+ `th_*` may be zero, in which case no pml is added.
+ omega: Angular frequency for the simulation
+ epsilon_effective: Effective epsilon of the PML. Match this to the material
+ at the edge of your grid.
+ Default 1.
+ s_function: created by `prepare_s_function(...)`, allowing customization of pml parameters.
+ Default uses `prepare_s_function()` with no parameters.
+
+ Returns:
+ Complex cell widths (dx_lists_t) as discussed in `meanas.types`.
"""
if s_function is None:
s_function = prepare_s_function()
@@ -88,21 +96,25 @@ def stretch_with_scpml(dxes: dx_lists_t,
omega: float,
epsilon_effective: float = 1.0,
thickness: int = 10,
- s_function: s_function_type = None,
+ s_function: s_function_t = None,
) -> dx_lists_t:
"""
Stretch dxes to contain a stretched-coordinate PML (SCPML) in one direction along one axis.
- :param dxes: dx_tuple with coordinates to stretch
- :param axis: axis to stretch (0=x, 1=y, 2=z)
- :param polarity: direction to stretch (-1 for -ve, +1 for +ve)
- :param omega: Angular frequency for the simulation
- :param epsilon_effective: Effective epsilon of the PML. Match this to the material at the
- edge of your grid. Default 1.
- :param thickness: number of cells to use for pml (default 10)
- :param s_function: s_function created by prepare_s_function(...), allowing customization
- of pml parameters. Default uses prepare_s_function() with no parameters.
- :return: Complex cell widths
+ Args:
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+ axis: axis to stretch (0=x, 1=y, 2=z)
+ polarity: direction to stretch (-1 for -ve, +1 for +ve)
+ omega: Angular frequency for the simulation
+ epsilon_effective: Effective epsilon of the PML. Match this to the material at the
+ edge of your grid. Default 1.
+ thickness: number of cells to use for pml (default 10)
+ s_function: Created by `prepare_s_function(...)`, allowing customization
+ of pml parameters. Default uses `prepare_s_function()` with no parameters.
+
+ Returns:
+ Complex cell widths (dx_lists_t) as discussed in `meanas.types`.
+ Multiple calls to this function may be necessary if multiple absorpbing boundaries are needed.
"""
if s_function is None:
s_function = prepare_s_function()
@@ -147,25 +159,3 @@ def stretch_with_scpml(dxes: dx_lists_t,
dxes[1][axis] = dx_bi
return dxes
-
-
-def generate_periodic_dx(pos: List[numpy.ndarray]) -> dx_lists_t:
- """
- Given a list of 3 ndarrays cell centers, creates the cell width parameters for a periodic grid.
-
- :param pos: List of 3 ndarrays of cell centers
- :return: (dx_a, dx_b) cell widths (no pml)
- """
- if len(pos) != 3:
- raise Exception('Must have len(pos) == 3')
-
- dx_a = [numpy.array(numpy.inf)] * 3
- dx_b = [numpy.array(numpy.inf)] * 3
-
- for i, p_orig in enumerate(pos):
- p = numpy.array(p_orig, dtype=float)
- if p.size != 1:
- p_shifted = numpy.hstack((p[1:], p[-1] + (p[1] - p[0])))
- dx_a[i] = numpy.diff(p)
- dx_b[i] = numpy.diff((p + p_shifted) / 2)
- return dx_a, dx_b
diff --git a/meanas/fdfd/solvers.py b/meanas/fdfd/solvers.py
index a0ce403..aa9633a 100644
--- a/meanas/fdfd/solvers.py
+++ b/meanas/fdfd/solvers.py
@@ -1,5 +1,5 @@
"""
-Solvers for FDFD problems.
+Solvers and solver interface for FDFD problems.
"""
from typing import List, Callable, Dict, Any
@@ -22,10 +22,13 @@ def _scipy_qmr(A: scipy.sparse.csr_matrix,
"""
Wrapper for scipy.sparse.linalg.qmr
- :param A: Sparse matrix
- :param b: Right-hand-side vector
- :param kwargs: Passed as **kwargs to the wrapped function
- :return: Guess for solution (returned even if didn't converge)
+ Args:
+ A: Sparse matrix
+ b: Right-hand-side vector
+ kwargs: Passed as **kwargs to the wrapped function
+
+ Returns:
+ Guess for solution (returned even if didn't converge)
"""
'''
@@ -70,27 +73,31 @@ def generic(omega: complex,
"""
Conjugate gradient FDFD solver using CSR sparse matrices.
- All ndarray arguments should be 1D array, as returned by meanas.vec().
+ All ndarray arguments should be 1D arrays, as returned by `meanas.vec()`.
- :param omega: Complex frequency to solve at.
- :param dxes: [[dx_e, dy_e, dz_e], [dx_h, dy_h, dz_h]] (complex cell sizes)
- :param J: Electric current distribution (at E-field locations)
- :param epsilon: Dielectric constant distribution (at E-field locations)
- :param mu: Magnetic permeability distribution (at H-field locations)
- :param pec: Perfect electric conductor distribution
- (at E-field locations; non-zero value indicates PEC is present)
- :param pmc: Perfect magnetic conductor distribution
- (at H-field locations; non-zero value indicates PMC is present)
- :param adjoint: If true, solves the adjoint problem.
- :param matrix_solver: Called as matrix_solver(A, b, **matrix_solver_opts) -> x
- Where A: scipy.sparse.csr_matrix
- b: numpy.ndarray
- x: numpy.ndarray
- Default is a wrapped version of scipy.sparse.linalg.qmr()
- which doesn't return convergence info and logs the residual
- every 100 iterations.
- :param matrix_solver_opts: Passed as kwargs to matrix_solver(...)
- :return: E-field which solves the system.
+ Args:
+ omega: Complex frequency to solve at.
+ dxes: `[[dx_e, dy_e, dz_e], [dx_h, dy_h, dz_h]]` (complex cell sizes) as
+ discussed in `meanas.types`
+ J: Electric current distribution (at E-field locations)
+ epsilon: Dielectric constant distribution (at E-field locations)
+ mu: Magnetic permeability distribution (at H-field locations)
+ pec: Perfect electric conductor distribution
+ (at E-field locations; non-zero value indicates PEC is present)
+ pmc: Perfect magnetic conductor distribution
+ (at H-field locations; non-zero value indicates PMC is present)
+ adjoint: If true, solves the adjoint problem.
+ matrix_solver: Called as `matrix_solver(A, b, **matrix_solver_opts) -> x`,
+ where `A`: `scipy.sparse.csr_matrix`;
+ `b`: `numpy.ndarray`;
+ `x`: `numpy.ndarray`;
+ Default is a wrapped version of `scipy.sparse.linalg.qmr()`
+ which doesn't return convergence info and logs the residual
+ every 100 iterations.
+ matrix_solver_opts: Passed as kwargs to `matrix_solver(...)`
+
+ Returns:
+ E-field which solves the system.
"""
if matrix_solver_opts is None:
diff --git a/meanas/fdfd/waveguide.py b/meanas/fdfd/waveguide.py
deleted file mode 100644
index c24a471..0000000
--- a/meanas/fdfd/waveguide.py
+++ /dev/null
@@ -1,492 +0,0 @@
-"""
-Various operators and helper functions for solving for waveguide modes.
-
-Assuming a z-dependence of the from exp(-i * wavenumber * z), we can simplify Maxwell's
- equations in the absence of sources to the form
-
-A @ [H_x, H_y] = wavenumber**2 * [H_x, H_y]
-
-with A =
-omega**2 * epsilon * mu +
-epsilon * [[-Dy], [Dx]] / epsilon * [-Dy, Dx] +
-[[Dx], [Dy]] / mu * [Dx, Dy] * mu
-
-which is the form used in this file.
-
-As the z-dependence is known, all the functions in this file assume a 2D grid
- (ie. dxes = [[[dx_e_0, dx_e_1, ...], [dy_e_0, ...]], [[dx_h_0, ...], [dy_h_0, ...]]])
- with propagation along the z axis.
-"""
-# TODO update module docs
-
-from typing import List, Tuple
-import numpy
-from numpy.linalg import norm
-import scipy.sparse as sparse
-
-from .. import vec, unvec, dx_lists_t, field_t, vfield_t
-from . import operators
-
-
-__author__ = 'Jan Petykiewicz'
-
-
-def operator_e(omega: complex,
- dxes: dx_lists_t,
- epsilon: vfield_t,
- mu: vfield_t = None,
- ) -> sparse.spmatrix:
- if numpy.any(numpy.equal(mu, None)):
- mu = numpy.ones_like(epsilon)
-
- Dfx, Dfy = operators.deriv_forward(dxes[0])
- Dbx, Dby = operators.deriv_back(dxes[1])
-
- eps_parts = numpy.split(epsilon, 3)
- eps_xy = sparse.diags(numpy.hstack((eps_parts[0], eps_parts[1])))
- eps_z_inv = sparse.diags(1 / eps_parts[2])
-
- mu_parts = numpy.split(mu, 3)
- mu_yx = sparse.diags(numpy.hstack((mu_parts[1], mu_parts[0])))
- mu_z_inv = sparse.diags(1 / mu_parts[2])
-
- op = omega * omega * mu_yx @ eps_xy + \
- mu_yx @ sparse.vstack((-Dby, Dbx)) @ mu_z_inv @ sparse.hstack((-Dfy, Dfx)) + \
- sparse.vstack((Dfx, Dfy)) @ eps_z_inv @ sparse.hstack((Dbx, Dby)) @ eps_xy
- return op
-
-
-def operator_h(omega: complex,
- dxes: dx_lists_t,
- epsilon: vfield_t,
- mu: vfield_t = None,
- ) -> sparse.spmatrix:
- """
- 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 [H_x, H_y].
-
- This operator can be used to form an eigenvalue problem of the form
- A @ [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).
-
- :param omega: The angular frequency of the system
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
- :param epsilon: Vectorized dielectric constant grid
- :param mu: Vectorized magnetic permeability grid (default 1 everywhere)
- :return: Sparse matrix representation of the operator
- """
- if numpy.any(numpy.equal(mu, None)):
- mu = numpy.ones_like(epsilon)
-
- Dfx, Dfy = operators.deriv_forward(dxes[0])
- Dbx, Dby = operators.deriv_back(dxes[1])
-
- eps_parts = numpy.split(epsilon, 3)
- eps_yx = sparse.diags(numpy.hstack((eps_parts[1], eps_parts[0])))
- eps_z_inv = sparse.diags(1 / eps_parts[2])
-
- mu_parts = numpy.split(mu, 3)
- mu_xy = sparse.diags(numpy.hstack((mu_parts[0], mu_parts[1])))
- mu_z_inv = sparse.diags(1 / mu_parts[2])
-
- op = omega * omega * eps_yx @ mu_xy + \
- eps_yx @ sparse.vstack((-Dfy, Dfx)) @ eps_z_inv @ sparse.hstack((-Dby, Dbx)) + \
- sparse.vstack((Dbx, Dby)) @ mu_z_inv @ sparse.hstack((Dfx, Dfy)) @ mu_xy
-
- return op
-
-
-def normalized_fields_e(e_xy: numpy.ndarray,
- wavenumber: complex,
- omega: complex,
- dxes: dx_lists_t,
- epsilon: vfield_t,
- mu: vfield_t = None,
- prop_phase: float = 0,
- ) -> Tuple[vfield_t, vfield_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.
-
- :param e_xy: Vector containing E_x and E_y fields
- :param wavenumber: Wavenumber satisfying `operator_e(...) @ e_xy == wavenumber**2 * e_xy`
- :param omega: The angular frequency of the system
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
- :param epsilon: Vectorized dielectric constant grid
- :param mu: Vectorized magnetic permeability grid (default 1 everywhere)
- :param prop_phase: Phase shift (dz * corrected_wavenumber) over 1 cell in propagation direction.
- Default 0 (continuous propagation direction, i.e. dz->0).
- :return: Normalized, vectorized (e, h) containing all vector components.
- """
- e = exy2e(wavenumber=wavenumber, dxes=dxes, epsilon=epsilon) @ e_xy
- h = exy2h(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) @ e_xy
- e_norm, h_norm = _normalized_fields(e=e, h=h, omega=omega, dxes=dxes, epsilon=epsilon,
- mu=mu, prop_phase=prop_phase)
- return e_norm, h_norm
-
-
-def normalized_fields_h(h_xy: numpy.ndarray,
- wavenumber: complex,
- omega: complex,
- dxes: dx_lists_t,
- epsilon: vfield_t,
- mu: vfield_t = None,
- prop_phase: float = 0,
- ) -> Tuple[vfield_t, vfield_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.
-
- :param e_xy: Vector containing E_x and E_y fields
- :param wavenumber: Wavenumber satisfying `operator_e(...) @ e_xy == wavenumber**2 * e_xy`
- :param omega: The angular frequency of the system
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
- :param epsilon: Vectorized dielectric constant grid
- :param mu: Vectorized magnetic permeability grid (default 1 everywhere)
- :param dxes_prop: Grid cell width in the propagation direction. Default 0 (continuous).
- :return: Normalized, vectorized (e, h) containing all vector components.
- """
- e = hxy2e(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) @ h_xy
- h = hxy2h(wavenumber=wavenumber, dxes=dxes, mu=mu) @ h_xy
- e_norm, h_norm = _normalized_fields(e=e, h=h, omega=omega, dxes=dxes, epsilon=epsilon,
- mu=mu, prop_phase=prop_phase)
- return e_norm, h_norm
-
-
-def _normalized_fields(e: numpy.ndarray,
- h: numpy.ndarray,
- omega: complex,
- dxes: dx_lists_t,
- epsilon: vfield_t,
- mu: vfield_t = None,
- prop_phase: float = 0,
- ) -> Tuple[vfield_t, vfield_t]:
- # TODO documentation
- shape = [s.size for s in dxes[0]]
- dxes_real = [[numpy.real(d) for d in numpy.meshgrid(*dxes[v], indexing='ij')] for v in (0, 1)]
-
- E = unvec(e, shape)
- H = unvec(h, shape)
-
- # Find time-averaged Sz and normalize to it
- # H phase is adjusted by a half-cell forward shift for Yee cell, and 1-cell reverse shift for Poynting
- phase = numpy.exp(-1j * -prop_phase / 2)
- Sz_a = E[0] * numpy.conj(H[1] * phase) * dxes_real[0][1] * dxes_real[1][0]
- Sz_b = E[1] * numpy.conj(H[0] * phase) * dxes_real[0][0] * dxes_real[1][1]
- Sz_tavg = numpy.real(Sz_a.sum() - Sz_b.sum()) * 0.5 # 0.5 since E, H are assumed to be peak (not RMS) amplitudes
- assert Sz_tavg > 0, 'Found a mode propagating in the wrong direction! Sz_tavg={}'.format(Sz_tavg)
-
- energy = epsilon * e.conj() * e
-
- norm_amplitude = 1 / numpy.sqrt(Sz_tavg)
- norm_angle = -numpy.angle(e[energy.argmax()]) # Will randomly add a negative sign when mode is symmetric
-
- # Try to break symmetry to assign a consistent sign [experimental TODO]
- E_weighted = unvec(e * energy * numpy.exp(1j * norm_angle), shape)
- sign = numpy.sign(E_weighted[:, :max(shape[0]//2, 1), :max(shape[1]//2, 1)].real.sum())
-
- norm_factor = sign * norm_amplitude * numpy.exp(1j * norm_angle)
-
- e *= norm_factor
- h *= norm_factor
-
- return e, h
-
-
-def exy2h(wavenumber: complex,
- omega: complex,
- dxes: dx_lists_t,
- epsilon: vfield_t,
- mu: vfield_t = None
- ) -> sparse.spmatrix:
- """
- 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
-
- :param wavenumber: Wavenumber satisfying `operator_e(...) @ e_xy == wavenumber**2 * e_xy`
- :param omega: The angular frequency of the system
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
- :param epsilon: Vectorized dielectric constant grid
- :param mu: Vectorized magnetic permeability grid (default 1 everywhere)
- :return: Sparse matrix representing the operator
- """
- e2hop = e2h(wavenumber=wavenumber, omega=omega, dxes=dxes, mu=mu)
- return e2hop @ exy2e(wavenumber=wavenumber, dxes=dxes, epsilon=epsilon)
-
-
-def hxy2e(wavenumber: complex,
- omega: complex,
- dxes: dx_lists_t,
- epsilon: vfield_t,
- mu: vfield_t = None
- ) -> sparse.spmatrix:
- """
- 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
-
- :param wavenumber: Wavenumber satisfying `operator_h(...) @ h_xy == wavenumber**2 * h_xy`
- :param omega: The angular frequency of the system
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
- :param epsilon: Vectorized dielectric constant grid
- :param mu: Vectorized magnetic permeability grid (default 1 everywhere)
- :return: Sparse matrix representing the operator
- """
- h2eop = h2e(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon)
- return h2eop @ hxy2h(wavenumber=wavenumber, dxes=dxes, mu=mu)
-
-
-def hxy2h(wavenumber: complex,
- dxes: dx_lists_t,
- mu: vfield_t = None
- ) -> sparse.spmatrix:
- """
- 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
-
- :param wavenumber: Wavenumber satisfying `operator_h(...) @ h_xy == wavenumber**2 * h_xy`
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
- :param mu: Vectorized magnetic permeability grid (default 1 everywhere)
- :return: Sparse matrix representing the operator
- """
- Dfx, Dfy = operators.deriv_forward(dxes[0])
- hxy2hz = sparse.hstack((Dfx, Dfy)) / (1j * wavenumber)
-
- if not numpy.any(numpy.equal(mu, None)):
- mu_parts = numpy.split(mu, 3)
- mu_xy = sparse.diags(numpy.hstack((mu_parts[0], mu_parts[1])))
- mu_z_inv = sparse.diags(1 / mu_parts[2])
-
- hxy2hz = mu_z_inv @ hxy2hz @ mu_xy
-
- n_pts = dxes[1][0].size * dxes[1][1].size
- op = sparse.vstack((sparse.eye(2 * n_pts),
- hxy2hz))
- return op
-
-
-def exy2e(wavenumber: complex,
- dxes: dx_lists_t,
- epsilon: vfield_t,
- ) -> sparse.spmatrix:
- """
- 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
-
- :param wavenumber: Wavenumber satisfying `operator_e(...) @ e_xy == wavenumber**2 * e_xy`
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
- :param epsilon: Vectorized dielectric constant grid
- :return: Sparse matrix representing the operator
- """
- Dbx, Dby = operators.deriv_back(dxes[1])
- exy2ez = sparse.hstack((Dbx, Dby)) / (1j * wavenumber)
-
- if not numpy.any(numpy.equal(epsilon, None)):
- epsilon_parts = numpy.split(epsilon, 3)
- epsilon_xy = sparse.diags(numpy.hstack((epsilon_parts[0], epsilon_parts[1])))
- epsilon_z_inv = sparse.diags(1 / epsilon_parts[2])
-
- exy2ez = epsilon_z_inv @ exy2ez @ epsilon_xy
-
- n_pts = dxes[0][0].size * dxes[0][1].size
- op = sparse.vstack((sparse.eye(2 * n_pts),
- exy2ez))
- return op
-
-
-def e2h(wavenumber: complex,
- omega: complex,
- dxes: dx_lists_t,
- mu: vfield_t = None
- ) -> sparse.spmatrix:
- """
- Returns an operator which, when applied to a vectorized E eigenfield, produces
- the vectorized H eigenfield.
-
- :param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v
- :param omega: The angular frequency of the system
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
- :param mu: Vectorized magnetic permeability grid (default 1 everywhere)
- :return: Sparse matrix representation of the operator
- """
- op = curl_e(wavenumber, dxes) / (-1j * omega)
- if not numpy.any(numpy.equal(mu, None)):
- op = sparse.diags(1 / mu) @ op
- return op
-
-
-def h2e(wavenumber: complex,
- omega: complex,
- dxes: dx_lists_t,
- epsilon: vfield_t
- ) -> sparse.spmatrix:
- """
- Returns an operator which, when applied to a vectorized H eigenfield, produces
- the vectorized E eigenfield.
-
- :param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v
- :param omega: The angular frequency of the system
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
- :param epsilon: Vectorized dielectric constant grid
- :return: Sparse matrix representation of the operator
- """
- op = sparse.diags(1 / (1j * omega * epsilon)) @ curl_h(wavenumber, dxes)
- return op
-
-
-def curl_e(wavenumber: complex, dxes: dx_lists_t) -> sparse.spmatrix:
- """
- Discretized curl operator for use with the waveguide E field.
-
- :param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
- :return: Sparse matrix representation of the operator
- """
- n = 1
- for d in dxes[0]:
- n *= len(d)
-
- Bz = -1j * wavenumber * sparse.eye(n)
- Dfx, Dfy = operators.deriv_forward(dxes[0])
- return operators.cross([Dfx, Dfy, Bz])
-
-
-def curl_h(wavenumber: complex, dxes: dx_lists_t) -> sparse.spmatrix:
- """
- Discretized curl operator for use with the waveguide H field.
-
- :param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
- :return: Sparse matrix representation of the operator
- """
- n = 1
- for d in dxes[1]:
- n *= len(d)
-
- Bz = -1j * wavenumber * sparse.eye(n)
- Dbx, Dby = operators.deriv_back(dxes[1])
- return operators.cross([Dbx, Dby, Bz])
-
-
-def h_err(h: vfield_t,
- wavenumber: complex,
- omega: complex,
- dxes: dx_lists_t,
- epsilon: vfield_t,
- mu: vfield_t = None
- ) -> float:
- """
- Calculates the relative error in the H field
-
- :param h: Vectorized H field
- :param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v
- :param omega: The angular frequency of the system
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
- :param epsilon: Vectorized dielectric constant grid
- :param mu: Vectorized magnetic permeability grid (default 1 everywhere)
- :return: Relative error norm(OP @ h) / norm(h)
- """
- ce = curl_e(wavenumber, dxes)
- ch = curl_h(wavenumber, dxes)
-
- eps_inv = sparse.diags(1 / epsilon)
-
- if numpy.any(numpy.equal(mu, None)):
- op = ce @ eps_inv @ ch @ h - omega ** 2 * h
- else:
- op = ce @ eps_inv @ ch @ h - omega ** 2 * (mu * h)
-
- return norm(op) / norm(h)
-
-
-def e_err(e: vfield_t,
- wavenumber: complex,
- omega: complex,
- dxes: dx_lists_t,
- epsilon: vfield_t,
- mu: vfield_t = None
- ) -> float:
- """
- Calculates the relative error in the E field
-
- :param e: Vectorized E field
- :param wavenumber: Wavenumber satisfying A @ v == wavenumber**2 * v
- :param omega: The angular frequency of the system
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
- :param epsilon: Vectorized dielectric constant grid
- :param mu: Vectorized magnetic permeability grid (default 1 everywhere)
- :return: Relative error norm(OP @ e) / norm(e)
- """
- ce = curl_e(wavenumber, dxes)
- ch = curl_h(wavenumber, dxes)
-
- if numpy.any(numpy.equal(mu, None)):
- op = ch @ ce @ e - omega ** 2 * (epsilon * e)
- else:
- mu_inv = sparse.diags(1 / mu)
- op = ch @ mu_inv @ ce @ e - omega ** 2 * (epsilon * e)
-
- return norm(op) / norm(e)
-
-
-def cylindrical_operator(omega: complex,
- dxes: dx_lists_t,
- epsilon: vfield_t,
- r0: float,
- ) -> sparse.spmatrix:
- """
- Cylindrical coordinate waveguide operator of the form
-
- TODO
-
- 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] = wavenumber**2 * [E_r, E_y]
-
- which can then be solved for the eigenmodes of the system (an exp(-i * wavenumber * theta)
- theta-dependence is assumed for the fields).
-
- :param omega: The angular frequency of the system
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types (2D)
- :param epsilon: Vectorized dielectric constant grid
- :param r0: Radius of curvature for the simulation. This should be the minimum value of
- r within the simulation domain.
- :return: Sparse matrix representation of the operator
- """
-
- Dfx, Dfy = operators.deriv_forward(dxes[0])
- Dbx, Dby = operators.deriv_back(dxes[1])
-
- rx = r0 + numpy.cumsum(dxes[0][0])
- ry = r0 + dxes[0][0]/2.0 + numpy.cumsum(dxes[1][0])
- tx = rx/r0
- ty = ry/r0
-
- Tx = sparse.diags(vec(tx[:, None].repeat(dxes[0][1].size, axis=1)))
- Ty = sparse.diags(vec(ty[:, None].repeat(dxes[1][1].size, axis=1)))
-
- eps_parts = numpy.split(epsilon, 3)
- eps_x = sparse.diags(eps_parts[0])
- eps_y = sparse.diags(eps_parts[1])
- eps_z_inv = sparse.diags(1 / eps_parts[2])
-
- pa = sparse.vstack((Dfx, Dfy)) @ Tx @ eps_z_inv @ sparse.hstack((Dbx, Dby))
- pb = sparse.vstack((Dfx, Dfy)) @ Tx @ eps_z_inv @ sparse.hstack((Dby, Dbx))
- a0 = Ty @ eps_x + omega**-2 * Dby @ Ty @ Dfy
- a1 = Tx @ eps_y + omega**-2 * Dbx @ Ty @ Dfx
- b0 = Dbx @ Ty @ Dfy
- b1 = Dby @ Ty @ Dfx
-
- diag = sparse.block_diag
- op = (omega**2 * diag((Tx, Ty)) + pa) @ diag((a0, a1)) + \
- - (sparse.bmat(((None, Ty), (Tx, None))) + omega**-2 * pb) @ diag((b0, b1))
-
- return op
-
diff --git a/meanas/fdfd/waveguide_2d.py b/meanas/fdfd/waveguide_2d.py
new file mode 100644
index 0000000..6a7f24f
--- /dev/null
+++ b/meanas/fdfd/waveguide_2d.py
@@ -0,0 +1,607 @@
+"""
+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, ...]]]`).
+"""
+# TODO update module docs
+
+from typing import List, Tuple
+import numpy
+from numpy.linalg import norm
+import scipy.sparse as sparse
+
+from .. import vec, unvec, dx_lists_t, field_t, vfield_t
+from ..eigensolvers import signed_eigensolve, rayleigh_quotient_iteration
+from . import operators
+
+
+__author__ = 'Jan Petykiewicz'
+
+
+def operator_e(omega: complex,
+ dxes: dx_lists_t,
+ epsilon: vfield_t,
+ mu: vfield_t = None,
+ ) -> sparse.spmatrix:
+ """
+ Waveguide operator of the form
+
+ omega**2 * mu * epsilon +
+ mu * [[-Dy], [Dx]] / mu * [-Dy, Dx] +
+ [[Dx], [Dy]] / epsilon * [Dx, Dy] * epsilon
+
+ for use with a field vector of the form `cat([E_x, E_y])`.
+
+ More precisely, the operator is
+ $$ \\omega^2 \\mu_{yx} \\epsilon_{xy} +
+ \\mu_{yx} \\begin{bmatrix} -D_{by} \\\\
+ D_{bx} \\end{bmatrix} \\mu_z^{-1}
+ \\begin{bmatrix} -D_{fy} & D_{fx} \\end{bmatrix} +
+ \\begin{bmatrix} D_{fx} \\\\
+ D_{fy} \\end{bmatrix} \\epsilon_z^{-1}
+ \\begin{bmatrix} D_{bx} & D_{by} \\end{bmatrix} \\epsilon_{xy} $$
+
+ where
+ \\( \\epsilon_{xy} = \\begin{bmatrix}
+ \\epsilon_x & 0 \\\\
+ 0 & \\epsilon_y
+ \\end{bmatrix} \\),
+ \\( \\mu_{yx} = \\begin{bmatrix}
+ \\mu_y & 0 \\\\
+ 0 & \\mu_x
+ \\end{bmatrix} \\),
+ \\( D_{fx} \\) and \\( D_{bx} \\) are the forward and backward derivatives along x,
+ and each \\( \\epsilon_x, \\mu_y, \\) etc. is a diagonal matrix representing
+
+
+ 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).
+
+ Args:
+ omega: The angular frequency of the system.
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D)
+ epsilon: Vectorized dielectric constant grid
+ mu: Vectorized magnetic permeability grid (default 1 everywhere)
+
+ Returns:
+ Sparse matrix representation of the operator.
+ """
+ if numpy.any(numpy.equal(mu, None)):
+ mu = numpy.ones_like(epsilon)
+
+ Dfx, Dfy = operators.deriv_forward(dxes[0])
+ Dbx, Dby = operators.deriv_back(dxes[1])
+
+ eps_parts = numpy.split(epsilon, 3)
+ eps_xy = sparse.diags(numpy.hstack((eps_parts[0], eps_parts[1])))
+ eps_z_inv = sparse.diags(1 / eps_parts[2])
+
+ mu_parts = numpy.split(mu, 3)
+ mu_yx = sparse.diags(numpy.hstack((mu_parts[1], mu_parts[0])))
+ mu_z_inv = sparse.diags(1 / mu_parts[2])
+
+ op = omega * omega * mu_yx @ eps_xy + \
+ mu_yx @ sparse.vstack((-Dby, Dbx)) @ mu_z_inv @ sparse.hstack((-Dfy, Dfx)) + \
+ sparse.vstack((Dfx, Dfy)) @ eps_z_inv @ sparse.hstack((Dbx, Dby)) @ eps_xy
+ return op
+
+
+def operator_h(omega: complex,
+ dxes: dx_lists_t,
+ epsilon: vfield_t,
+ mu: vfield_t = None,
+ ) -> sparse.spmatrix:
+ """
+ Waveguide operator of the form
+
+ omega**2 * epsilon * mu +
+ epsilon * [[-Dy], [Dx]] / epsilon * [-Dy, Dx] +
+ [[Dx], [Dy]] / mu * [Dx, Dy] * mu
+
+ for use with a field vector of the form `cat([H_x, H_y])`.
+
+ More precisely, the operator is
+ $$ \\omega^2 \\epsilon_{yx} \\mu_{xy} +
+ \\epsilon_{yx} \\begin{bmatrix} -D_{fy} \\\\
+ D_{fx} \\end{bmatrix} \\epsilon_z^{-1}
+ \\begin{bmatrix} -D_{by} & D_{bx} \\end{bmatrix} +
+ \\begin{bmatrix} D_{bx} \\\\
+ D_{by} \\end{bmatrix} \\mu_z^{-1}
+ \\begin{bmatrix} D_{fx} & D_{fy} \\end{bmatrix} \\mu_{xy} $$
+
+ where
+ \\( \\epsilon_{yx} = \\begin{bmatrix}
+ \\epsilon_y & 0 \\\\
+ 0 & \\epsilon_x
+ \\end{bmatrix} \\),
+ \\( \\mu_{xy} = \\begin{bmatrix}
+ \\mu_x & 0 \\\\
+ 0 & \\mu_y
+ \\end{bmatrix} \\),
+ \\( D_{fx} \\) and \\( D_{bx} \\) are the forward and backward derivatives along x,
+ and each \\( \\epsilon_x, \\mu_y, \\) etc. is a diagonal matrix.
+
+
+ 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).
+
+ Args:
+ omega: The angular frequency of the system.
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D)
+ epsilon: Vectorized dielectric constant grid
+ mu: Vectorized magnetic permeability grid (default 1 everywhere)
+
+ Returns:
+ Sparse matrix representation of the operator.
+ """
+ if numpy.any(numpy.equal(mu, None)):
+ mu = numpy.ones_like(epsilon)
+
+ Dfx, Dfy = operators.deriv_forward(dxes[0])
+ Dbx, Dby = operators.deriv_back(dxes[1])
+
+ eps_parts = numpy.split(epsilon, 3)
+ eps_yx = sparse.diags(numpy.hstack((eps_parts[1], eps_parts[0])))
+ eps_z_inv = sparse.diags(1 / eps_parts[2])
+
+ mu_parts = numpy.split(mu, 3)
+ mu_xy = sparse.diags(numpy.hstack((mu_parts[0], mu_parts[1])))
+ mu_z_inv = sparse.diags(1 / mu_parts[2])
+
+ op = omega * omega * eps_yx @ mu_xy + \
+ eps_yx @ sparse.vstack((-Dfy, Dfx)) @ eps_z_inv @ sparse.hstack((-Dby, Dbx)) + \
+ sparse.vstack((Dbx, Dby)) @ mu_z_inv @ sparse.hstack((Dfx, Dfy)) @ mu_xy
+
+ return op
+
+
+def normalized_fields_e(e_xy: numpy.ndarray,
+ wavenumber: complex,
+ omega: complex,
+ dxes: dx_lists_t,
+ epsilon: vfield_t,
+ mu: vfield_t = None,
+ prop_phase: float = 0,
+ ) -> Tuple[vfield_t, vfield_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.
+
+ Args:
+ e_xy: Vector containing E_x and E_y fields
+ wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`.
+ It should satisfy `operator_e() @ e_xy == wavenumber**2 * e_xy`
+ omega: The angular frequency of the system
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D)
+ epsilon: Vectorized dielectric constant grid
+ mu: Vectorized magnetic permeability grid (default 1 everywhere)
+ prop_phase: Phase shift `(dz * corrected_wavenumber)` over 1 cell in propagation direction.
+ Default 0 (continuous propagation direction, i.e. dz->0).
+
+ Returns:
+ `(e, h)`, where each field is vectorized, normalized,
+ and contains all three vector components.
+ """
+ e = exy2e(wavenumber=wavenumber, dxes=dxes, epsilon=epsilon) @ e_xy
+ h = exy2h(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) @ e_xy
+ e_norm, h_norm = _normalized_fields(e=e, h=h, omega=omega, dxes=dxes, epsilon=epsilon,
+ mu=mu, prop_phase=prop_phase)
+ return e_norm, h_norm
+
+
+def normalized_fields_h(h_xy: numpy.ndarray,
+ wavenumber: complex,
+ omega: complex,
+ dxes: dx_lists_t,
+ epsilon: vfield_t,
+ mu: vfield_t = None,
+ prop_phase: float = 0,
+ ) -> Tuple[vfield_t, vfield_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.
+
+ Args:
+ h_xy: Vector containing H_x and H_y fields
+ wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`.
+ It should satisfy `operator_h() @ h_xy == wavenumber**2 * h_xy`
+ omega: The angular frequency of the system
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D)
+ epsilon: Vectorized dielectric constant grid
+ mu: Vectorized magnetic permeability grid (default 1 everywhere)
+ prop_phase: Phase shift `(dz * corrected_wavenumber)` over 1 cell in propagation direction.
+ Default 0 (continuous propagation direction, i.e. dz->0).
+
+ Returns:
+ `(e, h)`, where each field is vectorized, normalized,
+ and contains all three vector components.
+ """
+ e = hxy2e(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) @ h_xy
+ h = hxy2h(wavenumber=wavenumber, dxes=dxes, mu=mu) @ h_xy
+ e_norm, h_norm = _normalized_fields(e=e, h=h, omega=omega, dxes=dxes, epsilon=epsilon,
+ mu=mu, prop_phase=prop_phase)
+ return e_norm, h_norm
+
+
+def _normalized_fields(e: numpy.ndarray,
+ h: numpy.ndarray,
+ omega: complex,
+ dxes: dx_lists_t,
+ epsilon: vfield_t,
+ mu: vfield_t = None,
+ prop_phase: float = 0,
+ ) -> Tuple[vfield_t, vfield_t]:
+ # TODO documentation
+ shape = [s.size for s in dxes[0]]
+ dxes_real = [[numpy.real(d) for d in numpy.meshgrid(*dxes[v], indexing='ij')] for v in (0, 1)]
+
+ E = unvec(e, shape)
+ H = unvec(h, shape)
+
+ # Find time-averaged Sz and normalize to it
+ # H phase is adjusted by a half-cell forward shift for Yee cell, and 1-cell reverse shift for Poynting
+ phase = numpy.exp(-1j * -prop_phase / 2)
+ Sz_a = E[0] * numpy.conj(H[1] * phase) * dxes_real[0][1] * dxes_real[1][0]
+ Sz_b = E[1] * numpy.conj(H[0] * phase) * dxes_real[0][0] * dxes_real[1][1]
+ Sz_tavg = numpy.real(Sz_a.sum() - Sz_b.sum()) * 0.5 # 0.5 since E, H are assumed to be peak (not RMS) amplitudes
+ assert Sz_tavg > 0, 'Found a mode propagating in the wrong direction! Sz_tavg={}'.format(Sz_tavg)
+
+ energy = epsilon * e.conj() * e
+
+ norm_amplitude = 1 / numpy.sqrt(Sz_tavg)
+ norm_angle = -numpy.angle(e[energy.argmax()]) # Will randomly add a negative sign when mode is symmetric
+
+ # Try to break symmetry to assign a consistent sign [experimental TODO]
+ E_weighted = unvec(e * energy * numpy.exp(1j * norm_angle), shape)
+ sign = numpy.sign(E_weighted[:, :max(shape[0]//2, 1), :max(shape[1]//2, 1)].real.sum())
+
+ norm_factor = sign * norm_amplitude * numpy.exp(1j * norm_angle)
+
+ e *= norm_factor
+ h *= norm_factor
+
+ return e, h
+
+
+def exy2h(wavenumber: complex,
+ omega: complex,
+ dxes: dx_lists_t,
+ epsilon: vfield_t,
+ mu: vfield_t = None
+ ) -> sparse.spmatrix:
+ """
+ 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
+
+ Args:
+ wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`.
+ It should satisfy `operator_e() @ e_xy == wavenumber**2 * e_xy`
+ omega: The angular frequency of the system
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D)
+ epsilon: Vectorized dielectric constant grid
+ mu: Vectorized magnetic permeability grid (default 1 everywhere)
+
+ Returns:
+ Sparse matrix representing the operator.
+ """
+ e2hop = e2h(wavenumber=wavenumber, omega=omega, dxes=dxes, mu=mu)
+ return e2hop @ exy2e(wavenumber=wavenumber, dxes=dxes, epsilon=epsilon)
+
+
+def hxy2e(wavenumber: complex,
+ omega: complex,
+ dxes: dx_lists_t,
+ epsilon: vfield_t,
+ mu: vfield_t = None
+ ) -> sparse.spmatrix:
+ """
+ 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
+
+ Args:
+ wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`.
+ It should satisfy `operator_h() @ h_xy == wavenumber**2 * h_xy`
+ omega: The angular frequency of the system
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D)
+ epsilon: Vectorized dielectric constant grid
+ mu: Vectorized magnetic permeability grid (default 1 everywhere)
+
+ Returns:
+ Sparse matrix representing the operator.
+ """
+ h2eop = h2e(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon)
+ return h2eop @ hxy2h(wavenumber=wavenumber, dxes=dxes, mu=mu)
+
+
+def hxy2h(wavenumber: complex,
+ dxes: dx_lists_t,
+ mu: vfield_t = None
+ ) -> sparse.spmatrix:
+ """
+ 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
+
+ Args:
+ wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`.
+ It should satisfy `operator_h() @ h_xy == wavenumber**2 * h_xy`
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D)
+ mu: Vectorized magnetic permeability grid (default 1 everywhere)
+
+ Returns:
+ Sparse matrix representing the operator.
+ """
+ Dfx, Dfy = operators.deriv_forward(dxes[0])
+ hxy2hz = sparse.hstack((Dfx, Dfy)) / (1j * wavenumber)
+
+ if not numpy.any(numpy.equal(mu, None)):
+ mu_parts = numpy.split(mu, 3)
+ mu_xy = sparse.diags(numpy.hstack((mu_parts[0], mu_parts[1])))
+ mu_z_inv = sparse.diags(1 / mu_parts[2])
+
+ hxy2hz = mu_z_inv @ hxy2hz @ mu_xy
+
+ n_pts = dxes[1][0].size * dxes[1][1].size
+ op = sparse.vstack((sparse.eye(2 * n_pts),
+ hxy2hz))
+ return op
+
+
+def exy2e(wavenumber: complex,
+ dxes: dx_lists_t,
+ epsilon: vfield_t,
+ ) -> sparse.spmatrix:
+ """
+ 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
+
+ Args:
+ wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`
+ It should satisfy `operator_e() @ e_xy == wavenumber**2 * e_xy`
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D)
+ epsilon: Vectorized dielectric constant grid
+
+ Returns:
+ Sparse matrix representing the operator.
+ """
+ Dbx, Dby = operators.deriv_back(dxes[1])
+ exy2ez = sparse.hstack((Dbx, Dby)) / (1j * wavenumber)
+
+ if not numpy.any(numpy.equal(epsilon, None)):
+ epsilon_parts = numpy.split(epsilon, 3)
+ epsilon_xy = sparse.diags(numpy.hstack((epsilon_parts[0], epsilon_parts[1])))
+ epsilon_z_inv = sparse.diags(1 / epsilon_parts[2])
+
+ exy2ez = epsilon_z_inv @ exy2ez @ epsilon_xy
+
+ n_pts = dxes[0][0].size * dxes[0][1].size
+ op = sparse.vstack((sparse.eye(2 * n_pts),
+ exy2ez))
+ return op
+
+
+def e2h(wavenumber: complex,
+ omega: complex,
+ dxes: dx_lists_t,
+ mu: vfield_t = None
+ ) -> sparse.spmatrix:
+ """
+ Returns an operator which, when applied to a vectorized E eigenfield, produces
+ the vectorized H eigenfield.
+
+ Args:
+ wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`
+ omega: The angular frequency of the system
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D)
+ mu: Vectorized magnetic permeability grid (default 1 everywhere)
+
+ Returns:
+ Sparse matrix representation of the operator.
+ """
+ op = curl_e(wavenumber, dxes) / (-1j * omega)
+ if not numpy.any(numpy.equal(mu, None)):
+ op = sparse.diags(1 / mu) @ op
+ return op
+
+
+def h2e(wavenumber: complex,
+ omega: complex,
+ dxes: dx_lists_t,
+ epsilon: vfield_t
+ ) -> sparse.spmatrix:
+ """
+ Returns an operator which, when applied to a vectorized H eigenfield, produces
+ the vectorized E eigenfield.
+
+ Args:
+ wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`
+ omega: The angular frequency of the system
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D)
+ epsilon: Vectorized dielectric constant grid
+
+ Returns:
+ Sparse matrix representation of the operator.
+ """
+ op = sparse.diags(1 / (1j * omega * epsilon)) @ curl_h(wavenumber, dxes)
+ return op
+
+
+def curl_e(wavenumber: complex, dxes: dx_lists_t) -> sparse.spmatrix:
+ """
+ Discretized curl operator for use with the waveguide E field.
+
+ Args:
+ wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D)
+
+ Return:
+ Sparse matrix representation of the operator.
+ """
+ n = 1
+ for d in dxes[0]:
+ n *= len(d)
+
+ Bz = -1j * wavenumber * sparse.eye(n)
+ Dfx, Dfy = operators.deriv_forward(dxes[0])
+ return operators.cross([Dfx, Dfy, Bz])
+
+
+def curl_h(wavenumber: complex, dxes: dx_lists_t) -> sparse.spmatrix:
+ """
+ Discretized curl operator for use with the waveguide H field.
+
+ Args:
+ wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D)
+
+ Return:
+ Sparse matrix representation of the operator.
+ """
+ n = 1
+ for d in dxes[1]:
+ n *= len(d)
+
+ Bz = -1j * wavenumber * sparse.eye(n)
+ Dbx, Dby = operators.deriv_back(dxes[1])
+ return operators.cross([Dbx, Dby, Bz])
+
+
+def h_err(h: vfield_t,
+ wavenumber: complex,
+ omega: complex,
+ dxes: dx_lists_t,
+ epsilon: vfield_t,
+ mu: vfield_t = None
+ ) -> float:
+ """
+ Calculates the relative error in the H field
+
+ Args:
+ h: Vectorized H field
+ wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`
+ omega: The angular frequency of the system
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D)
+ epsilon: Vectorized dielectric constant grid
+ mu: Vectorized magnetic permeability grid (default 1 everywhere)
+
+ Returns:
+ Relative error `norm(A_h @ h) / norm(h)`.
+ """
+ ce = curl_e(wavenumber, dxes)
+ ch = curl_h(wavenumber, dxes)
+
+ eps_inv = sparse.diags(1 / epsilon)
+
+ if numpy.any(numpy.equal(mu, None)):
+ op = ce @ eps_inv @ ch @ h - omega ** 2 * h
+ else:
+ op = ce @ eps_inv @ ch @ h - omega ** 2 * (mu * h)
+
+ return norm(op) / norm(h)
+
+
+def e_err(e: vfield_t,
+ wavenumber: complex,
+ omega: complex,
+ dxes: dx_lists_t,
+ epsilon: vfield_t,
+ mu: vfield_t = None
+ ) -> float:
+ """
+ Calculates the relative error in the E field
+
+ Args:
+ e: Vectorized E field
+ wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`
+ omega: The angular frequency of the system
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D)
+ epsilon: Vectorized dielectric constant grid
+ mu: Vectorized magnetic permeability grid (default 1 everywhere)
+
+ Returns:
+ Relative error `norm(A_e @ e) / norm(e)`.
+ """
+ ce = curl_e(wavenumber, dxes)
+ ch = curl_h(wavenumber, dxes)
+
+ if numpy.any(numpy.equal(mu, None)):
+ op = ch @ ce @ e - omega ** 2 * (epsilon * e)
+ else:
+ mu_inv = sparse.diags(1 / mu)
+ op = ch @ mu_inv @ ce @ e - omega ** 2 * (epsilon * e)
+
+ return norm(op) / norm(e)
+
+
+def solve_modes(mode_numbers: List[int],
+ omega: complex,
+ dxes: dx_lists_t,
+ epsilon: vfield_t,
+ mu: vfield_t = None,
+ mode_margin: int = 2,
+ ) -> Tuple[List[vfield_t], List[complex]]:
+ """
+ Given a 2D region, attempts to solve for the eigenmode with the specified mode numbers.
+
+ Args:
+ mode_numbers: List of 0-indexed mode numbers to solve for
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+ epsilon: Dielectric constant
+ mu: Magnetic permeability (default 1 everywhere)
+ mode_margin: The eigensolver will actually solve for `(max(mode_number) + mode_margin)`
+ modes, but only return the target mode. Increasing this value can improve the solver's
+ ability to find the correct mode. Default 2.
+
+ Returns:
+ (e_xys, wavenumbers)
+ """
+
+ '''
+ Solve for the largest-magnitude eigenvalue of the real operator
+ '''
+ dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes]
+ A_r = waveguide.operator_e(numpy.real(omega), dxes_real, numpy.real(epsilon), numpy.real(mu))
+
+ eigvals, eigvecs = signed_eigensolve(A_r, max(mode_number) + mode_margin)
+ e_xys = eigvecs[:, -(numpy.array(mode_number) + 1)]
+
+ '''
+ Now solve for the eigenvector of the full operator, using the real operator's
+ eigenvector as an initial guess for Rayleigh quotient iteration.
+ '''
+ A = waveguide.operator_e(omega, dxes, epsilon, mu)
+ eigvals, e_xys = rayleigh_quotient_iteration(A, e_xys)
+
+ # Calculate the wave-vector (force the real part to be positive)
+ wavenumbers = numpy.sqrt(eigvals)
+ wavenumbers *= numpy.sign(numpy.real(wavenumbers))
+
+ return e_xys, wavenumbers
+
+
+def solve_mode(mode_number: int,
+ *args,
+ **kwargs
+ ) -> Tuple[vfield_t, complex]:
+ """
+ Wrapper around `solve_modes()` that solves for a single mode.
+
+ Args:
+ mode_number: 0-indexed mode number to solve for
+ *args: passed to `solve_modes()`
+ **kwargs: passed to `solve_modes()`
+
+ Returns:
+ (e_xy, wavenumber)
+ """
+ return solve_modes(mode_numbers=[mode_number], *args, **kwargs)
diff --git a/meanas/fdfd/waveguide_3d.py b/meanas/fdfd/waveguide_3d.py
new file mode 100644
index 0000000..18727ce
--- /dev/null
+++ b/meanas/fdfd/waveguide_3d.py
@@ -0,0 +1,236 @@
+"""
+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.
+"""
+from typing import Dict, List, Tuple
+import numpy
+import scipy.sparse as sparse
+
+from .. import vec, unvec, dx_lists_t, vfield_t, field_t
+from . import operators, waveguide_2d, functional
+
+
+def solve_mode(mode_number: int,
+ omega: complex,
+ dxes: dx_lists_t,
+ axis: int,
+ polarity: int,
+ slices: List[slice],
+ epsilon: field_t,
+ mu: field_t = None,
+ ) -> Dict[str, complex or numpy.ndarray]:
+ """
+ Given a 3D grid, selects a slice from the grid and attempts to
+ solve for an eigenmode propagating through that slice.
+
+ Args:
+ mode_number: Number of the mode, 0-indexed
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+ axis: Propagation axis (0=x, 1=y, 2=z)
+ polarity: Propagation direction (+1 for +ve, -1 for -ve)
+ slices: `epsilon[tuple(slices)]` is used to select the portion of the grid to use
+ as the waveguide cross-section. `slices[axis]` should select only one item.
+ epsilon: Dielectric constant
+ mu: Magnetic permeability (default 1 everywhere)
+
+ Returns:
+ `{'E': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex}`
+ """
+ if mu is None:
+ mu = numpy.ones_like(epsilon)
+
+ slices = tuple(slices)
+
+ '''
+ Solve the 2D problem in the specified plane
+ '''
+ # Define rotation to set z as propagation direction
+ order = numpy.roll(range(3), 2 - axis)
+ reverse_order = numpy.roll(range(3), axis - 2)
+
+ # Find dx in propagation direction
+ dxab_forward = numpy.array([dx[order[2]][slices[order[2]]] for dx in dxes])
+ dx_prop = 0.5 * sum(dxab_forward)[0]
+
+ # Reduce to 2D and solve the 2D problem
+ args_2d = {
+ 'omega': omega,
+ 'dxes': [[dx[i][slices[i]] for i in order[:2]] for dx in dxes],
+ 'epsilon': vec([epsilon[i][slices].transpose(order) for i in order]),
+ 'mu': vec([mu[i][slices].transpose(order) for i in order]),
+ }
+ e_xy, wavenumber_2d = waveguide_2d.solve_mode(mode_number, **args_2d)
+
+ '''
+ Apply corrections and expand to 3D
+ '''
+ # Correct wavenumber to account for numerical dispersion.
+ wavenumber = 2/dx_prop * numpy.arcsin(wavenumber_2d * dx_prop/2)
+
+ shape = [d.size for d in args_2d['dxes'][0]]
+ ve, vh = waveguide.normalized_fields_e(e_xy, wavenumber=wavenumber_2d, **args_2d, prop_phase=dx_prop * wavenumber)
+ e = unvec(ve, shape)
+ h = unvec(vh, shape)
+
+ # Adjust for propagation direction
+ h *= polarity
+
+ # Apply phase shift to H-field
+ h[:2] *= numpy.exp(-1j * polarity * 0.5 * wavenumber * dx_prop)
+ e[2] *= numpy.exp(-1j * polarity * 0.5 * wavenumber * dx_prop)
+
+ # Expand E, H to full epsilon space we were given
+ E = numpy.zeros_like(epsilon, dtype=complex)
+ H = numpy.zeros_like(epsilon, dtype=complex)
+ for a, o in enumerate(reverse_order):
+ E[(a, *slices)] = e[o][:, :, None].transpose(reverse_order)
+ H[(a, *slices)] = h[o][:, :, None].transpose(reverse_order)
+
+ results = {
+ 'wavenumber': wavenumber,
+ 'wavenumber_2d': wavenumber_2d,
+ 'H': H,
+ 'E': E,
+ }
+ return results
+
+
+def compute_source(E: field_t,
+ wavenumber: complex,
+ omega: complex,
+ dxes: dx_lists_t,
+ axis: int,
+ polarity: int,
+ slices: List[slice],
+ epsilon: field_t,
+ mu: field_t = None,
+ ) -> field_t:
+ """
+ Given an eigenmode obtained by `solve_mode`, returns the current source distribution
+ necessary to position a unidirectional source at the slice location.
+
+ Args:
+ E: E-field of the mode
+ wavenumber: Wavenumber of the mode
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+ axis: Propagation axis (0=x, 1=y, 2=z)
+ polarity: Propagation direction (+1 for +ve, -1 for -ve)
+ slices: `epsilon[tuple(slices)]` is used to select the portion of the grid to use
+ as the waveguide cross-section. `slices[axis]` should select only one item.
+ mu: Magnetic permeability (default 1 everywhere)
+
+ Returns:
+ J distribution for the unidirectional source
+ """
+ E_expanded = expand_e(E=E, dxes=dxes, wavenumber=wavenumber, axis=axis,
+ polarity=polarity, slices=slices)
+
+ smask = [slice(None)] * 4
+ if polarity > 0:
+ smask[axis + 1] = slice(slices[axis].start, None)
+ else:
+ smask[axis + 1] = slice(None, slices[axis].stop)
+
+ mask = numpy.zeros_like(E_expanded, dtype=int)
+ mask[tuple(smask)] = 1
+
+ masked_e2j = operators.e_boundary_source(mask=vec(mask), omega=omega, dxes=dxes, epsilon=vec(epsilon), mu=vec(mu))
+ J = unvec(masked_e2j @ vec(E_expanded), E.shape[1:])
+ return J
+
+
+def compute_overlap_e(E: field_t,
+ wavenumber: complex,
+ dxes: dx_lists_t,
+ axis: int,
+ polarity: int,
+ slices: List[slice],
+ ) -> field_t: # TODO DOCS
+ """
+ Given an eigenmode obtained by `solve_mode`, calculates an overlap_e for the
+ mode orthogonality relation Integrate(((E x H_mode) + (E_mode x H)) dot dn)
+ [assumes reflection symmetry].
+
+ Args:
+ E: E-field of the mode
+ H: H-field of the mode (advanced by half of a Yee cell from E)
+ wavenumber: Wavenumber of the mode
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+ axis: Propagation axis (0=x, 1=y, 2=z)
+ polarity: Propagation direction (+1 for +ve, -1 for -ve)
+ slices: `epsilon[tuple(slices)]` is used to select the portion of the grid to use
+ as the waveguide cross-section. slices[axis] should select only one item.
+ mu: Magnetic permeability (default 1 everywhere)
+
+ Returns:
+ overlap_e such that `numpy.sum(overlap_e * other_e)` computes the overlap integral
+ """
+ slices = tuple(slices)
+
+ Ee = expand_e(E=E, wavenumber=wavenumber, dxes=dxes,
+ axis=axis, polarity=polarity, slices=slices)
+
+ start, stop = sorted((slices[axis].start, slices[axis].start - 2 * polarity))
+
+ slices2 = list(slices)
+ slices2[axis] = slice(start, stop)
+ slices2 = (slice(None), *slices2)
+
+ Etgt = numpy.zeros_like(Ee)
+ Etgt[slices2] = Ee[slices2]
+
+ Etgt /= (Etgt.conj() * Etgt).sum()
+ return Etgt
+
+
+def expand_e(E: field_t,
+ wavenumber: complex,
+ dxes: dx_lists_t,
+ axis: int,
+ polarity: int,
+ slices: List[slice],
+ ) -> field_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.
+
+ Args:
+ E: E-field of the mode
+ wavenumber: Wavenumber of the mode
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types`
+ axis: Propagation axis (0=x, 1=y, 2=z)
+ polarity: Propagation direction (+1 for +ve, -1 for -ve)
+ slices: `epsilon[tuple(slices)]` is used to select the portion of the grid to use
+ as the waveguide cross-section. slices[axis] should select only one item.
+
+ Returns:
+ `E`, with the original field expanded along the specified `axis`.
+ """
+ slices = tuple(slices)
+
+ # Determine phase factors for parallel slices
+ a_shape = numpy.roll([1, -1, 1, 1], axis)
+ a_E = numpy.real(dxes[0][axis]).cumsum()
+ r_E = a_E - a_E[slices[axis]]
+ iphi = polarity * -1j * wavenumber
+ phase_E = numpy.exp(iphi * r_E).reshape(a_shape)
+
+ # Expand our slice to the entire grid using the phase factors
+ E_expanded = numpy.zeros_like(E)
+
+ slices_exp = list(slices)
+ slices_exp[axis] = slice(E.shape[axis + 1])
+ slices_exp = (slice(None), *slices_exp)
+
+ slices_in = (slice(None), *slices)
+
+ E_expanded[slices_exp] = phase_E * numpy.array(E)[slices_in]
+ return E_expanded
diff --git a/meanas/fdfd/waveguide_cyl.py b/meanas/fdfd/waveguide_cyl.py
new file mode 100644
index 0000000..ebfb41d
--- /dev/null
+++ b/meanas/fdfd/waveguide_cyl.py
@@ -0,0 +1,138 @@
+"""
+Operators and helper functions for cylindrical waveguides with unchanging cross-section.
+
+WORK IN PROGRESS, CURRENTLY BROKEN
+
+As the z-dependence is known, all the functions in this file assume a 2D grid
+ (i.e. `dxes = [[[dr_e_0, dx_e_1, ...], [dy_e_0, ...]], [[dr_h_0, ...], [dy_h_0, ...]]]`).
+"""
+# TODO update module docs
+
+from typing import List, Tuple, Dict
+import numpy
+from numpy.linalg import norm
+import scipy.sparse as sparse
+
+from .. import vec, unvec, dx_lists_t, field_t, vfield_t
+from ..eigensolvers import signed_eigensolve, rayleigh_quotient_iteration
+from . import operators
+
+
+__author__ = 'Jan Petykiewicz'
+
+
+def cylindrical_operator(omega: complex,
+ dxes: dx_lists_t,
+ epsilon: vfield_t,
+ r0: float,
+ ) -> sparse.spmatrix:
+ """
+ Cylindrical coordinate waveguide operator of the form
+
+ TODO
+
+ 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] = wavenumber**2 * [E_r, E_y]
+
+ which can then be solved for the eigenmodes of the system
+ (an `exp(-i * wavenumber * theta)` theta-dependence is assumed for the fields).
+
+ Args:
+ omega: The angular frequency of the system
+ dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.types` (2D)
+ epsilon: Vectorized dielectric constant grid
+ r0: Radius of curvature for the simulation. This should be the minimum value of
+ r within the simulation domain.
+
+ Returns:
+ Sparse matrix representation of the operator
+ """
+
+ Dfx, Dfy = operators.deriv_forward(dxes[0])
+ Dbx, Dby = operators.deriv_back(dxes[1])
+
+ rx = r0 + numpy.cumsum(dxes[0][0])
+ ry = r0 + dxes[0][0]/2.0 + numpy.cumsum(dxes[1][0])
+ tx = rx/r0
+ ty = ry/r0
+
+ Tx = sparse.diags(vec(tx[:, None].repeat(dxes[0][1].size, axis=1)))
+ Ty = sparse.diags(vec(ty[:, None].repeat(dxes[1][1].size, axis=1)))
+
+ eps_parts = numpy.split(epsilon, 3)
+ eps_x = sparse.diags(eps_parts[0])
+ eps_y = sparse.diags(eps_parts[1])
+ eps_z_inv = sparse.diags(1 / eps_parts[2])
+
+ pa = sparse.vstack((Dfx, Dfy)) @ Tx @ eps_z_inv @ sparse.hstack((Dbx, Dby))
+ pb = sparse.vstack((Dfx, Dfy)) @ Tx @ eps_z_inv @ sparse.hstack((Dby, Dbx))
+ a0 = Ty @ eps_x + omega**-2 * Dby @ Ty @ Dfy
+ a1 = Tx @ eps_y + omega**-2 * Dbx @ Ty @ Dfx
+ b0 = Dbx @ Ty @ Dfy
+ b1 = Dby @ Ty @ Dfx
+
+ diag = sparse.block_diag
+ op = (omega**2 * diag((Tx, Ty)) + pa) @ diag((a0, a1)) + \
+ - (sparse.bmat(((None, Ty), (Tx, None))) + omega**-2 * pb) @ diag((b0, b1))
+
+ return op
+
+
+def solve_mode(mode_number: int,
+ omega: complex,
+ dxes: dx_lists_t,
+ epsilon: vfield_t,
+ r0: float,
+ ) -> Dict[str, complex or field_t]:
+ """
+ TODO: fixup
+ Given a 2d (r, y) slice of epsilon, attempts to solve for the eigenmode
+ of the bent waveguide with the specified mode number.
+
+ Args:
+ mode_number: Number of the mode, 0-indexed
+ omega: Angular frequency of the simulation
+ dxes: Grid parameters [dx_e, dx_h] as described in meanas.types.
+ The first coordinate is assumed to be r, the second is y.
+ epsilon: Dielectric constant
+ r0: Radius of curvature for the simulation. This should be the minimum value of
+ r within the simulation domain.
+
+ Returns:
+ `{'E': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex}`
+ """
+
+ '''
+ Solve for the largest-magnitude eigenvalue of the real operator
+ '''
+ dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes]
+
+ A_r = waveguide.cylindrical_operator(numpy.real(omega), dxes_real, numpy.real(epsilon), r0)
+ eigvals, eigvecs = signed_eigensolve(A_r, mode_number + 3)
+ e_xy = eigvecs[:, -(mode_number+1)]
+
+ '''
+ Now solve for the eigenvector of the full operator, using the real operator's
+ eigenvector as an initial guess for Rayleigh quotient iteration.
+ '''
+ A = waveguide.cylindrical_operator(omega, dxes, epsilon, r0)
+ eigval, e_xy = rayleigh_quotient_iteration(A, e_xy)
+
+ # Calculate the wave-vector (force the real part to be positive)
+ wavenumber = numpy.sqrt(eigval)
+ wavenumber *= numpy.sign(numpy.real(wavenumber))
+
+ # TODO: Perform correction on wavenumber to account for numerical dispersion.
+
+ shape = [d.size for d in dxes[0]]
+ e_xy = numpy.hstack((e_xy, numpy.zeros(shape[0] * shape[1])))
+ fields = {
+ 'wavenumber': wavenumber,
+ 'E': unvec(e_xy, shape),
+# 'E': unvec(e, shape),
+# 'H': unvec(h, shape),
+ }
+
+ return fields
diff --git a/meanas/fdfd/waveguide_mode.py b/meanas/fdfd/waveguide_mode.py
deleted file mode 100644
index 4318cc0..0000000
--- a/meanas/fdfd/waveguide_mode.py
+++ /dev/null
@@ -1,307 +0,0 @@
-from typing import Dict, List, Tuple
-import numpy
-import scipy.sparse as sparse
-
-from .. import vec, unvec, dx_lists_t, vfield_t, field_t
-from . import operators, waveguide, functional
-from ..eigensolvers import signed_eigensolve, rayleigh_quotient_iteration
-
-
-def vsolve_waveguide_mode_2d(mode_number: int,
- omega: complex,
- dxes: dx_lists_t,
- epsilon: vfield_t,
- mu: vfield_t = None,
- mode_margin: int = 2,
- ) -> Tuple[vfield_t, complex]:
- """
- Given a 2d region, attempts to solve for the eigenmode with the specified mode number.
-
- :param mode_number: Number of the mode, 0-indexed.
- :param omega: Angular frequency of the simulation
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :param epsilon: Dielectric constant
- :param mu: Magnetic permeability (default 1 everywhere)
- :param mode_margin: The eigensolver will actually solve for (mode_number + mode_margin)
- modes, but only return the target mode. Increasing this value can improve the solver's
- ability to find the correct mode. Default 2.
- :return: (e_xy, wavenumber)
- """
-
- '''
- Solve for the largest-magnitude eigenvalue of the real operator
- '''
- dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes]
- A_r = waveguide.operator_e(numpy.real(omega), dxes_real, numpy.real(epsilon), numpy.real(mu))
-
- eigvals, eigvecs = signed_eigensolve(A_r, mode_number + mode_margin)
- e_xy = eigvecs[:, -(mode_number + 1)]
-
- '''
- Now solve for the eigenvector of the full operator, using the real operator's
- eigenvector as an initial guess for Rayleigh quotient iteration.
- '''
- A = waveguide.operator_e(omega, dxes, epsilon, mu)
- eigval, e_xy = rayleigh_quotient_iteration(A, e_xy)
-
- # Calculate the wave-vector (force the real part to be positive)
- wavenumber = numpy.sqrt(eigval)
- wavenumber *= numpy.sign(numpy.real(wavenumber))
-
- return e_xy, wavenumber
-
-
-
-def solve_waveguide_mode(mode_number: int,
- omega: complex,
- dxes: dx_lists_t,
- axis: int,
- polarity: int,
- slices: List[slice],
- epsilon: field_t,
- mu: field_t = None,
- ) -> Dict[str, complex or numpy.ndarray]:
- """
- Given a 3D grid, selects a slice from the grid and attempts to
- solve for an eigenmode propagating through that slice.
-
- :param mode_number: Number of the mode, 0-indexed
- :param omega: Angular frequency of the simulation
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :param axis: Propagation axis (0=x, 1=y, 2=z)
- :param polarity: Propagation direction (+1 for +ve, -1 for -ve)
- :param slices: epsilon[tuple(slices)] is used to select the portion of the grid to use
- as the waveguide cross-section. slices[axis] should select only one
- :param epsilon: Dielectric constant
- :param mu: Magnetic permeability (default 1 everywhere)
- :return: {'E': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex}
- """
- if mu is None:
- mu = numpy.ones_like(epsilon)
-
- slices = tuple(slices)
-
- '''
- Solve the 2D problem in the specified plane
- '''
- # Define rotation to set z as propagation direction
- order = numpy.roll(range(3), 2 - axis)
- reverse_order = numpy.roll(range(3), axis - 2)
-
- # Find dx in propagation direction
- dxab_forward = numpy.array([dx[order[2]][slices[order[2]]] for dx in dxes])
- dx_prop = 0.5 * sum(dxab_forward)[0]
-
- # Reduce to 2D and solve the 2D problem
- args_2d = {
- 'omega': omega,
- 'dxes': [[dx[i][slices[i]] for i in order[:2]] for dx in dxes],
- 'epsilon': vec([epsilon[i][slices].transpose(order) for i in order]),
- 'mu': vec([mu[i][slices].transpose(order) for i in order]),
- }
- e_xy, wavenumber_2d = vsolve_waveguide_mode_2d(mode_number, **args_2d)
-
- '''
- Apply corrections and expand to 3D
- '''
- # Correct wavenumber to account for numerical dispersion.
- wavenumber = 2/dx_prop * numpy.arcsin(wavenumber_2d * dx_prop/2)
- print(wavenumber_2d / wavenumber)
-
- shape = [d.size for d in args_2d['dxes'][0]]
- ve, vh = waveguide.normalized_fields_e(e_xy, wavenumber=wavenumber_2d, **args_2d, prop_phase=dx_prop * wavenumber)
- e = unvec(ve, shape)
- h = unvec(vh, shape)
-
- # Adjust for propagation direction
- h *= polarity
-
- # Apply phase shift to H-field
- h[:2] *= numpy.exp(-1j * polarity * 0.5 * wavenumber * dx_prop)
- e[2] *= numpy.exp(-1j * polarity * 0.5 * wavenumber * dx_prop)
-
- # Expand E, H to full epsilon space we were given
- E = numpy.zeros_like(epsilon, dtype=complex)
- H = numpy.zeros_like(epsilon, dtype=complex)
- for a, o in enumerate(reverse_order):
- E[(a, *slices)] = e[o][:, :, None].transpose(reverse_order)
- H[(a, *slices)] = h[o][:, :, None].transpose(reverse_order)
-
- results = {
- 'wavenumber': wavenumber,
- 'wavenumber_2d': wavenumber_2d,
- 'H': H,
- 'E': E,
- }
- return results
-
-
-def compute_source(E: field_t,
- wavenumber: complex,
- omega: complex,
- dxes: dx_lists_t,
- axis: int,
- polarity: int,
- slices: List[slice],
- epsilon: field_t,
- mu: field_t = None,
- ) -> field_t:
- """
- Given an eigenmode obtained by solve_waveguide_mode, returns the current source distribution
- necessary to position a unidirectional source at the slice location.
-
- :param E: E-field of the mode
- :param wavenumber: Wavenumber of the mode
- :param omega: Angular frequency of the simulation
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :param axis: Propagation axis (0=x, 1=y, 2=z)
- :param polarity: Propagation direction (+1 for +ve, -1 for -ve)
- :param slices: epsilon[tuple(slices)] is used to select the portion of the grid to use
- as the waveguide cross-section. slices[axis] should select only one
- :param mu: Magnetic permeability (default 1 everywhere)
- :return: J distribution for the unidirectional source
- """
- E_expanded = expand_wgmode_e(E=E, dxes=dxes, wavenumber=wavenumber, axis=axis,
- polarity=polarity, slices=slices)
-
- smask = [slice(None)] * 4
- if polarity > 0:
- smask[axis + 1] = slice(slices[axis].start, None)
- else:
- smask[axis + 1] = slice(None, slices[axis].stop)
-
- mask = numpy.zeros_like(E_expanded, dtype=int)
- mask[tuple(smask)] = 1
-
- masked_e2j = operators.e_boundary_source(mask=vec(mask), omega=omega, dxes=dxes, epsilon=vec(epsilon), mu=vec(mu))
- J = unvec(masked_e2j @ vec(E_expanded), E.shape[1:])
- return J
-
-
-def compute_overlap_e(E: field_t,
- wavenumber: complex,
- dxes: dx_lists_t,
- axis: int,
- polarity: int,
- slices: List[slice],
- ) -> field_t: # TODO DOCS
- """
- Given an eigenmode obtained by solve_waveguide_mode, calculates overlap_e for the
- mode orthogonality relation Integrate(((E x H_mode) + (E_mode x H)) dot dn)
- [assumes reflection symmetry].i
-
- overlap_e makes use of the e2h operator to collapse the above expression into
- (vec(E) @ vec(overlap_e)), allowing for simple calculation of the mode overlap.
-
- :param E: E-field of the mode
- :param H: H-field of the mode (advanced by half of a Yee cell from E)
- :param wavenumber: Wavenumber of the mode
- :param omega: Angular frequency of the simulation
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :param axis: Propagation axis (0=x, 1=y, 2=z)
- :param polarity: Propagation direction (+1 for +ve, -1 for -ve)
- :param slices: epsilon[tuple(slices)] is used to select the portion of the grid to use
- as the waveguide cross-section. slices[axis] should select only one
- :param mu: Magnetic permeability (default 1 everywhere)
- :return: overlap_e for calculating the mode overlap
- """
- slices = tuple(slices)
-
- Ee = expand_wgmode_e(E=E, wavenumber=wavenumber, dxes=dxes,
- axis=axis, polarity=polarity, slices=slices)
-
- start, stop = sorted((slices[axis].start, slices[axis].start - 2 * polarity))
-
- slices2 = list(slices)
- slices2[axis] = slice(start, stop)
- slices2 = (slice(None), *slices2)
-
- Etgt = numpy.zeros_like(Ee)
- Etgt[slices2] = Ee[slices2]
-
- Etgt /= (Etgt.conj() * Etgt).sum()
- return Etgt
-
-
-def solve_waveguide_mode_cylindrical(mode_number: int,
- omega: complex,
- dxes: dx_lists_t,
- epsilon: vfield_t,
- r0: float,
- ) -> Dict[str, complex or field_t]:
- """
- TODO: fixup
- Given a 2d (r, y) slice of epsilon, attempts to solve for the eigenmode
- of the bent waveguide with the specified mode number.
-
- :param mode_number: Number of the mode, 0-indexed
- :param omega: Angular frequency of the simulation
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types.
- The first coordinate is assumed to be r, the second is y.
- :param epsilon: Dielectric constant
- :param r0: Radius of curvature for the simulation. This should be the minimum value of
- r within the simulation domain.
- :return: {'E': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex}
- """
-
- '''
- Solve for the largest-magnitude eigenvalue of the real operator
- '''
- dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes]
-
- A_r = waveguide.cylindrical_operator(numpy.real(omega), dxes_real, numpy.real(epsilon), r0)
- eigvals, eigvecs = signed_eigensolve(A_r, mode_number + 3)
- e_xy = eigvecs[:, -(mode_number+1)]
-
- '''
- Now solve for the eigenvector of the full operator, using the real operator's
- eigenvector as an initial guess for Rayleigh quotient iteration.
- '''
- A = waveguide.cylindrical_operator(omega, dxes, epsilon, r0)
- eigval, e_xy = rayleigh_quotient_iteration(A, e_xy)
-
- # Calculate the wave-vector (force the real part to be positive)
- wavenumber = numpy.sqrt(eigval)
- wavenumber *= numpy.sign(numpy.real(wavenumber))
-
- # TODO: Perform correction on wavenumber to account for numerical dispersion.
-
- shape = [d.size for d in dxes[0]]
- e_xy = numpy.hstack((e_xy, numpy.zeros(shape[0] * shape[1])))
- fields = {
- 'wavenumber': wavenumber,
- 'E': unvec(e_xy, shape),
-# 'E': unvec(e, shape),
-# 'H': unvec(h, shape),
- }
-
- return fields
-
-
-def expand_wgmode_e(E: field_t,
- wavenumber: complex,
- dxes: dx_lists_t,
- axis: int,
- polarity: int,
- slices: List[slice],
- ) -> field_t:
- slices = tuple(slices)
-
- # Determine phase factors for parallel slices
- a_shape = numpy.roll([1, -1, 1, 1], axis)
- a_E = numpy.real(dxes[0][axis]).cumsum()
- r_E = a_E - a_E[slices[axis]]
- iphi = polarity * -1j * wavenumber
- phase_E = numpy.exp(iphi * r_E).reshape(a_shape)
-
- # Expand our slice to the entire grid using the phase factors
- E_expanded = numpy.zeros_like(E)
-
- slices_exp = list(slices)
- slices_exp[axis] = slice(E.shape[axis + 1])
- slices_exp = (slice(None), *slices_exp)
-
- slices_in = (slice(None), *slices)
-
- E_expanded[slices_exp] = phase_E * numpy.array(E)[slices_in]
- return E_expanded
diff --git a/meanas/fdmath/functional.py b/meanas/fdmath/functional.py
new file mode 100644
index 0000000..8593d4a
--- /dev/null
+++ b/meanas/fdmath/functional.py
@@ -0,0 +1,109 @@
+"""
+Math functions for finite difference simulations
+
+Basic discrete calculus etc.
+"""
+from typing import List, Callable, Tuple, Dict
+import numpy
+
+from .. import field_t, field_updater
+
+
+def deriv_forward(dx_e: List[numpy.ndarray] = None) -> field_updater:
+ """
+ Utility operators for taking discretized derivatives (backward variant).
+
+ Args:
+ dx_e: Lists of cell sizes for all axes
+ `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`.
+
+ Returns:
+ List of functions for taking forward derivatives along each axis.
+ """
+ if dx_e:
+ derivs = [lambda f: (numpy.roll(f, -1, axis=0) - f) / dx_e[0][:, None, None],
+ lambda f: (numpy.roll(f, -1, axis=1) - f) / dx_e[1][None, :, None],
+ lambda f: (numpy.roll(f, -1, axis=2) - f) / dx_e[2][None, None, :]]
+ else:
+ derivs = [lambda f: numpy.roll(f, -1, axis=0) - f,
+ lambda f: numpy.roll(f, -1, axis=1) - f,
+ lambda f: numpy.roll(f, -1, axis=2) - f]
+ return derivs
+
+
+def deriv_back(dx_h: List[numpy.ndarray] = None) -> field_updater:
+ """
+ Utility operators for taking discretized derivatives (forward variant).
+
+ Args:
+ dx_h: Lists of cell sizes for all axes
+ `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`.
+
+ Returns:
+ List of functions for taking forward derivatives along each axis.
+ """
+ if dx_h:
+ derivs = [lambda f: (f - numpy.roll(f, 1, axis=0)) / dx_h[0][:, None, None],
+ lambda f: (f - numpy.roll(f, 1, axis=1)) / dx_h[1][None, :, None],
+ lambda f: (f - numpy.roll(f, 1, axis=2)) / dx_h[2][None, None, :]]
+ else:
+ derivs = [lambda f: f - numpy.roll(f, 1, axis=0),
+ lambda f: f - numpy.roll(f, 1, axis=1),
+ lambda f: f - numpy.roll(f, 1, axis=2)]
+ return derivs
+
+
+def curl_forward(dx_e: List[numpy.ndarray] = None) -> field_updater:
+ """
+ Curl operator for use with the E field.
+
+ Args:
+ dx_e: Lists of cell sizes for all axes
+ `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`.
+
+ Returns:
+ Function `f` for taking the discrete forward curl of a field,
+ `f(E)` -> curlE \\( = \\nabla_f \\times E \\)
+ """
+ Dx, Dy, Dz = deriv_forward(dx_e)
+
+ def ce_fun(e: field_t) -> field_t:
+ output = numpy.empty_like(e)
+ output[0] = Dy(e[2])
+ output[1] = Dz(e[0])
+ output[2] = Dx(e[1])
+ output[0] -= Dz(e[1])
+ output[1] -= Dx(e[2])
+ output[2] -= Dy(e[0])
+ return output
+
+ return ce_fun
+
+
+def curl_back(dx_h: List[numpy.ndarray] = None) -> field_updater:
+ """
+ Create a function which takes the backward curl of a field.
+
+ Args:
+ dx_h: Lists of cell sizes for all axes
+ `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`.
+
+ Returns:
+ Function `f` for taking the discrete backward curl of a field,
+ `f(H)` -> curlH \\( = \\nabla_b \\times H \\)
+ """
+ Dx, Dy, Dz = deriv_back(dx_h)
+
+ def ch_fun(h: field_t) -> field_t:
+ output = numpy.empty_like(h)
+ output[0] = Dy(h[2])
+ output[1] = Dz(h[0])
+ output[2] = Dx(h[1])
+ output[0] -= Dz(h[1])
+ output[1] -= Dx(h[2])
+ output[2] -= Dy(h[0])
+ return output
+
+ return ch_fun
+
+
diff --git a/meanas/fdmath/operators.py b/meanas/fdmath/operators.py
new file mode 100644
index 0000000..64f04aa
--- /dev/null
+++ b/meanas/fdmath/operators.py
@@ -0,0 +1,231 @@
+"""
+Matrix operators for finite difference simulations
+
+Basic discrete calculus etc.
+"""
+from typing import List, Callable, Tuple, Dict
+import numpy
+import scipy.sparse as sparse
+
+from .. import field_t, vfield_t
+
+
+def rotation(axis: int, shape: List[int], shift_distance: int=1) -> sparse.spmatrix:
+ """
+ Utility operator for performing a circular shift along a specified axis by a
+ specified number of elements.
+
+ Args:
+ axis: Axis to shift along. x=0, y=1, z=2
+ shape: Shape of the grid being shifted
+ shift_distance: Number of cells to shift by. May be negative. Default 1.
+
+ Returns:
+ Sparse matrix for performing the circular shift.
+ """
+ if len(shape) not in (2, 3):
+ raise Exception('Invalid shape: {}'.format(shape))
+ if axis not in range(len(shape)):
+ raise Exception('Invalid direction: {}, shape is {}'.format(axis, shape))
+
+ shifts = [abs(shift_distance) if a == axis else 0 for a in range(3)]
+ shifted_diags = [(numpy.arange(n) + s) % n for n, s in zip(shape, shifts)]
+ ijk = numpy.meshgrid(*shifted_diags, indexing='ij')
+
+ n = numpy.prod(shape)
+ i_ind = numpy.arange(n)
+ j_ind = numpy.ravel_multi_index(ijk, shape, order='C')
+
+ vij = (numpy.ones(n), (i_ind, j_ind.ravel(order='C')))
+
+ d = sparse.csr_matrix(vij, shape=(n, n))
+
+ if shift_distance < 0:
+ d = d.T
+
+ return d
+
+
+def shift_with_mirror(axis: int, shape: List[int], shift_distance: int=1) -> sparse.spmatrix:
+ """
+ Utility operator for performing an n-element shift along a specified axis, with mirror
+ boundary conditions applied to the cells beyond the receding edge.
+
+ Args:
+ axis: Axis to shift along. x=0, y=1, z=2
+ shape: Shape of the grid being shifted
+ shift_distance: Number of cells to shift by. May be negative. Default 1.
+
+ Returns:
+ Sparse matrix for performing the shift-with-mirror.
+ """
+ if len(shape) not in (2, 3):
+ raise Exception('Invalid shape: {}'.format(shape))
+ if axis not in range(len(shape)):
+ raise Exception('Invalid direction: {}, shape is {}'.format(axis, shape))
+ if shift_distance >= shape[axis]:
+ raise Exception('Shift ({}) is too large for axis {} of size {}'.format(
+ shift_distance, axis, shape[axis]))
+
+ def mirrored_range(n, s):
+ v = numpy.arange(n) + s
+ v = numpy.where(v >= n, 2 * n - v - 1, v)
+ v = numpy.where(v < 0, - 1 - v, v)
+ return v
+
+ shifts = [shift_distance if a == axis else 0 for a in range(3)]
+ shifted_diags = [mirrored_range(n, s) for n, s in zip(shape, shifts)]
+ ijk = numpy.meshgrid(*shifted_diags, indexing='ij')
+
+ n = numpy.prod(shape)
+ i_ind = numpy.arange(n)
+ j_ind = numpy.ravel_multi_index(ijk, shape, order='C')
+
+ vij = (numpy.ones(n), (i_ind, j_ind.ravel(order='C')))
+
+ d = sparse.csr_matrix(vij, shape=(n, n))
+ return d
+
+
+def deriv_forward(dx_e: List[numpy.ndarray]) -> List[sparse.spmatrix]:
+ """
+ Utility operators for taking discretized derivatives (forward variant).
+
+ Args:
+ dx_e: Lists of cell sizes for all axes
+ `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`.
+
+ Returns:
+ List of operators for taking forward derivatives along each axis.
+ """
+ shape = [s.size for s in dx_e]
+ n = numpy.prod(shape)
+
+ dx_e_expanded = numpy.meshgrid(*dx_e, indexing='ij')
+
+ def deriv(axis):
+ return rotation(axis, shape, 1) - sparse.eye(n)
+
+ Ds = [sparse.diags(+1 / dx.ravel(order='C')) @ deriv(a)
+ for a, dx in enumerate(dx_e_expanded)]
+
+ return Ds
+
+
+def deriv_back(dx_h: List[numpy.ndarray]) -> List[sparse.spmatrix]:
+ """
+ Utility operators for taking discretized derivatives (backward variant).
+
+ Args:
+ dx_h: Lists of cell sizes for all axes
+ `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`.
+
+ Returns:
+ List of operators for taking forward derivatives along each axis.
+ """
+ shape = [s.size for s in dx_h]
+ n = numpy.prod(shape)
+
+ dx_h_expanded = numpy.meshgrid(*dx_h, indexing='ij')
+
+ def deriv(axis):
+ return rotation(axis, shape, -1) - sparse.eye(n)
+
+ Ds = [sparse.diags(-1 / dx.ravel(order='C')) @ deriv(a)
+ for a, dx in enumerate(dx_h_expanded)]
+
+ return Ds
+
+
+def cross(B: List[sparse.spmatrix]) -> sparse.spmatrix:
+ """
+ Cross product operator
+
+ Args:
+ B: List `[Bx, By, Bz]` of sparse matrices corresponding to the x, y, z
+ portions of the operator on the left side of the cross product.
+
+ Returns:
+ Sparse matrix corresponding to (B x), where x is the cross product.
+ """
+ n = B[0].shape[0]
+ zero = sparse.csr_matrix((n, n))
+ return sparse.bmat([[zero, -B[2], B[1]],
+ [B[2], zero, -B[0]],
+ [-B[1], B[0], zero]])
+
+
+def vec_cross(b: vfield_t) -> sparse.spmatrix:
+ """
+ Vector cross product operator
+
+ Args:
+ b: Vector on the left side of the cross product.
+
+ Returns:
+
+ Sparse matrix corresponding to (b x), where x is the cross product.
+
+ """
+ B = [sparse.diags(c) for c in numpy.split(b, 3)]
+ return cross(B)
+
+
+def avg_forward(axis: int, shape: List[int]) -> sparse.spmatrix:
+ """
+ Forward average operator `(x4 = (x4 + x5) / 2)`
+
+ Args:
+ axis: Axis to average along (x=0, y=1, z=2)
+ shape: Shape of the grid to average
+
+ Returns:
+ Sparse matrix for forward average operation.
+ """
+ if len(shape) not in (2, 3):
+ raise Exception('Invalid shape: {}'.format(shape))
+
+ n = numpy.prod(shape)
+ return 0.5 * (sparse.eye(n) + rotation(axis, shape))
+
+
+def avg_back(axis: int, shape: List[int]) -> sparse.spmatrix:
+ """
+ Backward average operator `(x4 = (x4 + x3) / 2)`
+
+ Args:
+ axis: Axis to average along (x=0, y=1, z=2)
+ shape: Shape of the grid to average
+
+ Returns:
+ Sparse matrix for backward average operation.
+ """
+ return avg_forward(axis, shape).T
+
+
+def curl_forward(dx_e: List[numpy.ndarray]) -> sparse.spmatrix:
+ """
+ Curl operator for use with the E field.
+
+ Args:
+ dx_e: Lists of cell sizes for all axes
+ `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`.
+
+ Returns:
+ Sparse matrix for taking the discretized curl of the E-field
+ """
+ return cross(deriv_forward(dx_e))
+
+
+def curl_back(dx_h: List[numpy.ndarray]) -> sparse.spmatrix:
+ """
+ Curl operator for use with the H field.
+
+ Args:
+ dx_h: Lists of cell sizes for all axes
+ `[[dx_0, dx_1, ...], [dy_0, dy_1, ...], ...]`.
+
+ Returns:
+ Sparse matrix for taking the discretized curl of the H-field
+ """
+ return cross(deriv_back(dx_h))
diff --git a/meanas/fdtd/__init__.py b/meanas/fdtd/__init__.py
index a1d278a..02cfc05 100644
--- a/meanas/fdtd/__init__.py
+++ b/meanas/fdtd/__init__.py
@@ -1,5 +1,5 @@
"""
-Basic FDTD functionality
+Utilities for running finite-difference time-domain (FDTD) simulations
"""
from .base import maxwell_e, maxwell_h
diff --git a/meanas/fdtd/base.py b/meanas/fdtd/base.py
index 389a7b5..87a234b 100644
--- a/meanas/fdtd/base.py
+++ b/meanas/fdtd/base.py
@@ -5,70 +5,14 @@ from typing import List, Callable, Tuple, Dict
import numpy
from .. import dx_lists_t, field_t, field_updater
+from ..fdmath.functional import curl_forward, curl_back
+
__author__ = 'Jan Petykiewicz'
-def curl_h(dxes: dx_lists_t = None) -> field_updater:
- """
- Curl operator for use with the H field.
-
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :return: Function for taking the discretized curl of the H-field, F(H) -> curlH
- """
- if dxes:
- dxyz_b = numpy.meshgrid(*dxes[1], indexing='ij')
-
- def dh(f, ax):
- return (f - numpy.roll(f, 1, axis=ax)) / dxyz_b[ax]
- else:
- def dh(f, ax):
- return f - numpy.roll(f, 1, axis=ax)
-
- def ch_fun(h: field_t) -> field_t:
- output = numpy.empty_like(h)
- output[0] = dh(h[2], 1)
- output[1] = dh(h[0], 2)
- output[2] = dh(h[1], 0)
- output[0] -= dh(h[1], 2)
- output[1] -= dh(h[2], 0)
- output[2] -= dh(h[0], 1)
- return output
-
- return ch_fun
-
-
-def curl_e(dxes: dx_lists_t = None) -> field_updater:
- """
- Curl operator for use with the E field.
-
- :param dxes: Grid parameters [dx_e, dx_h] as described in meanas.types
- :return: Function for taking the discretized curl of the E-field, F(E) -> curlE
- """
- if dxes is not None:
- dxyz_a = numpy.meshgrid(*dxes[0], indexing='ij')
-
- def de(f, ax):
- return (numpy.roll(f, -1, axis=ax) - f) / dxyz_a[ax]
- else:
- def de(f, ax):
- return numpy.roll(f, -1, axis=ax) - f
-
- def ce_fun(e: field_t) -> field_t:
- output = numpy.empty_like(e)
- output[0] = de(e[2], 1)
- output[1] = de(e[0], 2)
- output[2] = de(e[1], 0)
- output[0] -= de(e[1], 2)
- output[1] -= de(e[2], 0)
- output[2] -= de(e[0], 1)
- return output
-
- return ce_fun
-
-
def maxwell_e(dt: float, dxes: dx_lists_t = None) -> field_updater:
- curl_h_fun = curl_h(dxes)
+ curl_h_fun = curl_back(dxes[1])
def me_fun(e: field_t, h: field_t, epsilon: field_t):
e += dt * curl_h_fun(h) / epsilon
@@ -78,7 +22,7 @@ def maxwell_e(dt: float, dxes: dx_lists_t = None) -> field_updater:
def maxwell_h(dt: float, dxes: dx_lists_t = None) -> field_updater:
- curl_e_fun = curl_e(dxes)
+ curl_e_fun = curl_forward(dxes[0])
def mh_fun(e: field_t, h: field_t):
h -= dt * curl_e_fun(e)
diff --git a/meanas/fdtd/energy.py b/meanas/fdtd/energy.py
index 8644646..d5ea1dc 100644
--- a/meanas/fdtd/energy.py
+++ b/meanas/fdtd/energy.py
@@ -35,6 +35,7 @@ def poynting_divergence(s: field_t = None,
if s is None:
s = poynting(e, h, dxes=dxes)
+ #TODO use deriv operators
ds = ((s[0] - numpy.roll(s[0], 1, axis=0)) +
(s[1] - numpy.roll(s[1], 1, axis=1)) +
(s[2] - numpy.roll(s[2], 1, axis=2)))
diff --git a/meanas/test/__init__.py b/meanas/test/__init__.py
index e69de29..e02b636 100644
--- a/meanas/test/__init__.py
+++ b/meanas/test/__init__.py
@@ -0,0 +1,3 @@
+"""
+Tests (run with `python3 -m pytest -rxPXs | tee results.txt`)
+"""
diff --git a/meanas/test/test_fdfd_pml.py b/meanas/test/test_fdfd_pml.py
index 112193d..eafa14f 100644
--- a/meanas/test/test_fdfd_pml.py
+++ b/meanas/test/test_fdfd_pml.py
@@ -101,8 +101,8 @@ def j_distribution(request, shape, epsilon, dxes, omega, src_polarity):
slices[dim] = slice(shape[dim + 1] // 2,
shape[dim + 1] // 2 + 1)
- j = fdfd.waveguide_mode.compute_source(E=e, wavenumber=wavenumber_corrected, omega=omega, dxes=dxes,
- axis=dim, polarity=src_polarity, slices=slices, epsilon=epsilon)
+ j = fdfd.waveguide_3d.compute_source(E=e, wavenumber=wavenumber_corrected, omega=omega, dxes=dxes,
+ axis=dim, polarity=src_polarity, slices=slices, epsilon=epsilon)
yield j
@@ -145,4 +145,3 @@ def sim(request, shape, epsilon, dxes, j_distribution, omega, pec, pmc):
)
return sim
-
diff --git a/meanas/types.py b/meanas/types.py
index 7dc5c1c..667a286 100644
--- a/meanas/types.py
+++ b/meanas/types.py
@@ -6,17 +6,20 @@ from typing import List, Callable
# Field types
-field_t = numpy.ndarray # vector field with shape (3, X, Y, Z) (e.g. [E_x, E_y, E_z])
-vfield_t = numpy.ndarray # linearized vector field (vector of length 3*X*Y*Z)
+field_t = numpy.ndarray
+"""vector field with shape (3, X, Y, Z) (e.g. `[E_x, E_y, E_z]`)"""
+vfield_t = numpy.ndarray
+"""Linearized vector field (vector of length 3*X*Y*Z)"""
+
+dx_lists_t = List[List[numpy.ndarray]]
'''
'dxes' datastructure which contains grid cell width information in the following format:
- [[[dx_e_0, dx_e_1, ...], [dy_e_0, ...], [dz_e_0, ...]],
- [[dx_h_0, dx_h_1, ...], [dy_h_0, ...], [dz_h_0, ...]]]
- where dx_e_0 is the x-width of the x=0 cells, as used when calculating dE/dx,
- and dy_h_0 is the y-width of the y=0 cells, as used when calculating dH/dy, etc.
+ `[[[dx_e_0, dx_e_1, ...], [dy_e_0, ...], [dz_e_0, ...]],
+ [[dx_h_0, dx_h_1, ...], [dy_h_0, ...], [dz_h_0, ...]]]`
+ where `dx_e_0` is the x-width of the `x=0` cells, as used when calculating dE/dx,
+ and `dy_h_0` is the y-width of the `y=0` cells, as used when calculating dH/dy, etc.
'''
-dx_lists_t = List[List[numpy.ndarray]]
field_updater = Callable[[field_t], field_t]
diff --git a/pdoc_templates/config.mako b/pdoc_templates/config.mako
new file mode 100644
index 0000000..0642566
--- /dev/null
+++ b/pdoc_templates/config.mako
@@ -0,0 +1,46 @@
+<%!
+ # Template configuration. Copy over in your template directory
+ # (used with --template-dir) and adapt as required.
+ html_lang = 'en'
+ show_inherited_members = False
+ extract_module_toc_into_sidebar = True
+ list_class_variables_in_index = True
+ sort_identifiers = True
+ show_type_annotations = True
+
+ # Show collapsed source code block next to each item.
+ # Disabling this can improve rendering speed of large modules.
+ show_source_code = True
+
+ # If set, format links to objects in online source code repository
+ # according to this template. Supported keywords for interpolation
+ # are: commit, path, start_line, end_line.
+ #git_link_template = 'https://github.com/USER/PROJECT/blob/{commit}/{path}#L{start_line}-L{end_line}'
+ #git_link_template = 'https://gitlab.com/USER/PROJECT/blob/{commit}/{path}#L{start_line}-L{end_line}'
+ #git_link_template = 'https://bitbucket.org/USER/PROJECT/src/{commit}/{path}#lines-{start_line}:{end_line}'
+ #git_link_template = 'https://CGIT_HOSTNAME/PROJECT/tree/{path}?id={commit}#n{start-line}'
+ git_link_template = None
+
+ # A prefix to use for every HTML hyperlink in the generated documentation.
+ # No prefix results in all links being relative.
+ link_prefix = ''
+
+ # Enable syntax highlighting for code/source blocks by including Highlight.js
+ syntax_highlighting = True
+
+ # Set the style keyword such as 'atom-one-light' or 'github-gist'
+ # Options: https://github.com/highlightjs/highlight.js/tree/master/src/styles
+ # Demo: https://highlightjs.org/static/demo/
+ hljs_style = 'github'
+
+ # If set, insert Google Analytics tracking code. Value is GA
+ # tracking id (UA-XXXXXX-Y).
+ google_analytics = ''
+
+ # If set, render LaTeX math syntax within \(...\) (inline equations),
+ # or within \[...\] or $$...$$ or `.. math::` (block equations)
+ # as nicely-formatted math formulas using MathJax.
+ # Note: in Python docstrings, either all backslashes need to be escaped (\\)
+ # or you need to use raw r-strings.
+ latex_math = True
+%>
diff --git a/pdoc_templates/css.mako b/pdoc_templates/css.mako
new file mode 100644
index 0000000..39a77ed
--- /dev/null
+++ b/pdoc_templates/css.mako
@@ -0,0 +1,389 @@
+<%!
+ from pdoc.html_helpers import minify_css
+%>
+
+<%def name="mobile()" filter="minify_css">
+ .flex {
+ display: flex !important;
+ }
+
+ body {
+ line-height: 1.5em;
+ background: black;
+ color: #DDD;
+ }
+
+ #content {
+ padding: 20px;
+ }
+
+ #sidebar {
+ padding: 30px;
+ overflow: hidden;
+ }
+
+ .http-server-breadcrumbs {
+ font-size: 130%;
+ margin: 0 0 15px 0;
+ }
+
+ #footer {
+ font-size: .75em;
+ padding: 5px 30px;
+ border-top: 1px solid #ddd;
+ text-align: right;
+ }
+ #footer p {
+ margin: 0 0 0 1em;
+ display: inline-block;
+ }
+ #footer p:last-child {
+ margin-right: 30px;
+ }
+
+ h1, h2, h3, h4, h5 {
+ font-weight: 300;
+ }
+ h1 {
+ font-size: 2.5em;
+ line-height: 1.1em;
+ }
+ h2 {
+ font-size: 1.75em;
+ margin: 1em 0 .50em 0;
+ }
+ h3 {
+ font-size: 1.4em;
+ margin: 25px 0 10px 0;
+ }
+ h4 {
+ margin: 0;
+ font-size: 105%;
+ }
+
+ a {
+ color: #999;
+ text-decoration: none;
+ transition: color .3s ease-in-out;
+ }
+ a:hover {
+ color: #18d;
+ }
+
+ .title code {
+ font-weight: bold;
+ }
+ h2[id^="header-"] {
+ margin-top: 2em;
+ }
+ .ident {
+ color: #7ff;
+ }
+
+ pre code {
+ background: transparent;
+ font-size: .8em;
+ line-height: 1.4em;
+ }
+ code {
+ background: #0d0d0e;
+ padding: 1px 4px;
+ overflow-wrap: break-word;
+ }
+ h1 code { background: transparent }
+
+ pre {
+ background: #111;
+ border: 0;
+ border-top: 1px solid #ccc;
+ border-bottom: 1px solid #ccc;
+ margin: 1em 0;
+ padding: 1ex;
+ }
+
+ #http-server-module-list {
+ display: flex;
+ flex-flow: column;
+ }
+ #http-server-module-list div {
+ display: flex;
+ }
+ #http-server-module-list dt {
+ min-width: 10%;
+ }
+ #http-server-module-list p {
+ margin-top: 0;
+ }
+
+ .toc ul,
+ #index {
+ list-style-type: none;
+ margin: 0;
+ padding: 0;
+ }
+ #index code {
+ background: transparent;
+ }
+ #index h3 {
+ border-bottom: 1px solid #ddd;
+ }
+ #index ul {
+ padding: 0;
+ }
+ #index h4 {
+ font-weight: bold;
+ }
+ #index h4 + ul {
+ margin-bottom:.6em;
+ }
+ /* Make TOC lists have 2+ columns when viewport is wide enough.
+ Assuming ~20-character identifiers and ~30% wide sidebar. */
+ @media (min-width: 200ex) { #index .two-column { column-count: 2 } }
+ @media (min-width: 300ex) { #index .two-column { column-count: 3 } }
+
+ dl {
+ margin-bottom: 2em;
+ }
+ dl dl:last-child {
+ margin-bottom: 4em;
+ }
+ dd {
+ margin: 0 0 1em 3em;
+ }
+ #header-classes + dl > dd {
+ margin-bottom: 3em;
+ }
+ dd dd {
+ margin-left: 2em;
+ }
+ dd p {
+ margin: 10px 0;
+ }
+ .name {
+ background: #111;
+ font-weight: bold;
+ font-size: .85em;
+ padding: 5px 10px;
+ display: inline-block;
+ min-width: 40%;
+ }
+ .name:hover {
+ background: #101010;
+ }
+ .name > span:first-child {
+ white-space: nowrap;
+ }
+ .name.class > span:nth-child(2) {
+ margin-left: .4em;
+ }
+ .inherited {
+ color: #777;
+ border-left: 5px solid #eee;
+ padding-left: 1em;
+ }
+ .inheritance em {
+ font-style: normal;
+ font-weight: bold;
+ }
+
+ /* Docstrings titles, e.g. in numpydoc format */
+ .desc h2 {
+ font-weight: 400;
+ font-size: 1.25em;
+ }
+ .desc h3 {
+ font-size: 1em;
+ }
+ .desc dt code {
+ background: inherit; /* Don't grey-back parameters */
+ }
+
+ .source summary,
+ .git-link-div {
+ color: #aaa;
+ text-align: right;
+ font-weight: 400;
+ font-size: .8em;
+ text-transform: uppercase;
+ }
+ .source summary > * {
+ white-space: nowrap;
+ cursor: pointer;
+ }
+ .git-link {
+ color: inherit;
+ margin-left: 1em;
+ }
+ .source pre {
+ max-height: 500px;
+ overflow: auto;
+ margin: 0;
+ }
+ .source pre code {
+ font-size: 12px;
+ overflow: visible;
+ }
+ .hlist {
+ list-style: none;
+ }
+ .hlist li {
+ display: inline;
+ }
+ .hlist li:after {
+ content: ',\2002';
+ }
+ .hlist li:last-child:after {
+ content: none;
+ }
+ .hlist .hlist {
+ display: inline;
+ padding-left: 1em;
+ }
+
+ img {
+ max-width: 100%;
+ }
+
+ .admonition {
+ padding: .1em .5em;
+ margin-bottom: 1em;
+ }
+ .admonition-title {
+ font-weight: bold;
+ }
+ .admonition.note,
+ .admonition.info,
+ .admonition.important {
+ background: #610;
+ }
+ .admonition.todo,
+ .admonition.versionadded,
+ .admonition.tip,
+ .admonition.hint {
+ background: #202;
+ }
+ .admonition.warning,
+ .admonition.versionchanged,
+ .admonition.deprecated {
+ background: #02b;
+ }
+ .admonition.error,
+ .admonition.danger,
+ .admonition.caution {
+ background: darkpink;
+ }
+
+
+<%def name="desktop()" filter="minify_css">
+ @media screen and (min-width: 700px) {
+ #sidebar {
+ width: 30%;
+ }
+ #content {
+ width: 70%;
+ max-width: 100ch;
+ padding: 3em 4em;
+ border-left: 1px solid #ddd;
+ }
+ pre code {
+ font-size: 1em;
+ }
+ .item .name {
+ font-size: 1em;
+ }
+ main {
+ display: flex;
+ flex-direction: row-reverse;
+ justify-content: flex-end;
+ }
+ .toc ul ul,
+ #index ul {
+ padding-left: 1.5em;
+ }
+ .toc > ul > li {
+ margin-top: .5em;
+ }
+ }
+
+
+<%def name="print()" filter="minify_css">
+@media print {
+ #sidebar h1 {
+ page-break-before: always;
+ }
+ .source {
+ display: none;
+ }
+}
+@media print {
+ * {
+ background: transparent !important;
+ color: #000 !important; /* Black prints faster: h5bp.com/s */
+ box-shadow: none !important;
+ text-shadow: none !important;
+ }
+
+ a[href]:after {
+ content: " (" attr(href) ")";
+ font-size: 90%;
+ }
+ /* Internal, documentation links, recognized by having a title,
+ don't need the URL explicity stated. */
+ a[href][title]:after {
+ content: none;
+ }
+
+ abbr[title]:after {
+ content: " (" attr(title) ")";
+ }
+
+ /*
+ * Don't show links for images, or javascript/internal links
+ */
+
+ .ir a:after,
+ a[href^="javascript:"]:after,
+ a[href^="#"]:after {
+ content: "";
+ }
+
+ pre,
+ blockquote {
+ border: 1px solid #999;
+ page-break-inside: avoid;
+ }
+
+ thead {
+ display: table-header-group; /* h5bp.com/t */
+ }
+
+ tr,
+ img {
+ page-break-inside: avoid;
+ }
+
+ img {
+ max-width: 100% !important;
+ }
+
+ @page {
+ margin: 0.5cm;
+ }
+
+ p,
+ h2,
+ h3 {
+ orphans: 3;
+ widows: 3;
+ }
+
+ h1,
+ h2,
+ h3,
+ h4,
+ h5,
+ h6 {
+ page-break-after: avoid;
+ }
+}
+
diff --git a/pdoc_templates/html.mako b/pdoc_templates/html.mako
new file mode 100644
index 0000000..9cf1137
--- /dev/null
+++ b/pdoc_templates/html.mako
@@ -0,0 +1,421 @@
+<%
+ import os
+
+ import pdoc
+ from pdoc.html_helpers import extract_toc, glimpse, to_html as _to_html, format_git_link
+
+
+ def link(d, name=None, fmt='{}'):
+ name = fmt.format(name or d.qualname + ('()' if isinstance(d, pdoc.Function) else ''))
+ if not isinstance(d, pdoc.Doc) or isinstance(d, pdoc.External) and not external_links:
+ return name
+ url = d.url(relative_to=module, link_prefix=link_prefix,
+ top_ancestor=not show_inherited_members)
+ return '{}'.format(d.refname, url, name)
+
+
+ def to_html(text):
+ return _to_html(text, module=module, link=link, latex_math=latex_math)
+%>
+
+<%def name="ident(name)">${name}
+
+<%def name="show_source(d)">
+ % if (show_source_code or git_link_template) and d.source and d.obj is not getattr(d.inherits, 'obj', None):
+ <% git_link = format_git_link(git_link_template, d) %>
+ % if show_source_code:
+
+
+ Expand source code
+ % if git_link:
+ Browse git
+ %endif
+
+ ${d.source | h}
+
+ % elif git_link:
+
+ %endif
+ %endif
+
+
+<%def name="show_desc(d, short=False)">
+ <%
+ inherits = ' inherited' if d.inherits else ''
+ docstring = glimpse(d.docstring) if short or inherits else d.docstring
+ %>
+ % if d.inherits:
+
+ Inherited from:
+ % if hasattr(d.inherits, 'cls'):
+ ${link(d.inherits.cls)}.${link(d.inherits, d.name)}
+ % else:
+ ${link(d.inherits)}
+ % endif
+
+ % endif
+
+ % if not isinstance(d, pdoc.Module):
+ ${show_source(d)}
+ % endif
+
+
+<%def name="show_module_list(modules)">
+Python module list
+
+% if not modules:
+ No modules found.
+% else:
+
+ % for name, desc in modules:
+
+
- ${name}
+
- ${desc | glimpse, to_html}
+
+% endif
+
+
+<%def name="show_column_list(items)">
+ <%
+ two_column = len(items) >= 6 and all(len(i.name) < 20 for i in items)
+ %>
+
+ % for item in items:
+ ${link(item, item.name)}
+
+
+<%def name="show_module(module)">
+ <%
+ variables = module.variables(sort=sort_identifiers)
+ classes = module.classes(sort=sort_identifiers)
+ functions = module.functions(sort=sort_identifiers)
+ submodules = module.submodules()
+ %>
+
+ <%def name="show_func(f)">
+
+ <%
+ params = ', '.join(f.params(annotate=show_type_annotations, link=link))
+ returns = show_type_annotations and f.return_annotation(link=link) or ''
+ if returns:
+ returns = ' ->\N{NBSP}' + returns
+ %>
+ ${f.funcdef()} ${ident(f.name)}(${params})${returns}
+
+ ${show_desc(f)}
+
+
+
+ % if http_server:
+
+ % endif
+ ${'Namespace' if module.is_namespace else 'Module'} ${module.name}
+
+
+
+ ${module.docstring | to_html}
+ ${show_source(module)}
+
+
+
+ % if submodules:
+
+
+ % for m in submodules:
+ ${link(m)}- ${show_desc(m, short=True)}
+ % endfor
+
+ % endif
+
+
+
+ % if variables:
+
+
+ % for v in variables:
+ var ${ident(v.name)}- ${show_desc(v)}
+ % endfor
+
+ % endif
+
+
+
+ % if functions:
+
+
+ % for f in functions:
+ ${show_func(f)}
+ % endfor
+
+ % endif
+
+
+
+ % if classes:
+
+
+ % for c in classes:
+ <%
+ class_vars = c.class_variables(show_inherited_members, sort=sort_identifiers)
+ smethods = c.functions(show_inherited_members, sort=sort_identifiers)
+ inst_vars = c.instance_variables(show_inherited_members, sort=sort_identifiers)
+ methods = c.methods(show_inherited_members, sort=sort_identifiers)
+ mro = c.mro()
+ subclasses = c.subclasses()
+ params = ', '.join(c.params(annotate=show_type_annotations, link=link))
+ %>
+
+ class ${ident(c.name)}
+ % if params:
+ (${params})
+ % endif
+ - ${show_desc(c)}
+
+ % if mro:
+
Ancestors
+
+ % for cls in mro:
+ - ${link(cls)}
+ % endfor
+
+ %endif
+
+ % if subclasses:
+ Subclasses
+
+ % for sub in subclasses:
+ - ${link(sub)}
+ % endfor
+
+ % endif
+ % if class_vars:
+ Class variables
+
+ % for v in class_vars:
+ var ${ident(v.name)}- ${show_desc(v)}
+ % endfor
+
+ % endif
+ % if smethods:
+ Static methods
+
+ % for f in smethods:
+ ${show_func(f)}
+ % endfor
+
+ % endif
+ % if inst_vars:
+ Instance variables
+
+ % for v in inst_vars:
+ var ${ident(v.name)}- ${show_desc(v)}
+ % endfor
+
+ % endif
+ % if methods:
+ Methods
+
+ % for f in methods:
+ ${show_func(f)}
+ % endfor
+
+ % endif
+
+ % if not show_inherited_members:
+ <%
+ members = c.inherited_members()
+ %>
+ % if members:
+ Inherited members
+
+ % for cls, mems in members:
+ ${link(cls)}:
+
+ % for m in mems:
+ ${link(m, name=m.name)}
+
+
+ % endif
+ % endif
+
+
+ % endfor
+
+ % endif
+
+
+
+<%def name="module_index(module)">
+ <%
+ variables = module.variables(sort=sort_identifiers)
+ classes = module.classes(sort=sort_identifiers)
+ functions = module.functions(sort=sort_identifiers)
+ submodules = module.submodules()
+ supermodule = module.supermodule
+ %>
+
+
+
+
+
+
+
+
+
+
+<%
+ module_list = 'modules' in context.keys() # Whether we're showing module list in server mode
+%>
+
+ % if module_list:
+ Python module list
+
+ % else:
+ ${module.name} API documentation
+
+ % endif
+
+
+
+ % if syntax_highlighting:
+
+ %endif
+
+ <%namespace name="css" file="css.mako" />
+
+
+
+
+ % if google_analytics:
+
+ % endif
+
+ % if latex_math:
+
+ % endif
+
+ <%include file="head.mako"/>
+
+
+
+ % if module_list:
+
+ ${show_module_list(modules)}
+
+ % else:
+
+ ${show_module(module)}
+
+ ${module_index(module)}
+ % endif
+
+
+
+
+% if syntax_highlighting:
+
+
+% endif
+
+% if http_server and module: ## Auto-reload on file change in dev mode
+
+% endif
+
+