[meanas.fdtd.misc] add basic pulse and beam shapes

- Properly account for rmin vs r0
- Change return values to match waveguide_2d
- Change operator definition to look more like waveguide_2d

remaining TODO:
- Fix docs
- Further consolidate operators vs waveguide_2d
- Figure out E/H field conversions

examples/fdfd.py

@@ -46,20 +46,24 @@ def test0(solver=generic_solver):
     # #### Create the grid, mask, and draw the device ####
     grid = gridlock.Grid(edge_coords)
     epsilon = grid.allocate(n_air**2, dtype=numpy.float32)
-    grid.draw_cylinder(epsilon,
+    grid.draw_cylinder(
+        epsilon,
         surface_normal=2,
         center=center,
         radius=max(radii),
         thickness=th,
-                       eps=n_ring**2,
-                       num_points=24)
-    grid.draw_cylinder(epsilon,
+        foreground=n_ring**2,
+        num_points=24,
+        )
+    grid.draw_cylinder(
+        epsilon,
         surface_normal=2,
         center=center,
         radius=min(radii),
         thickness=th*1.1,
-                       eps=n_air ** 2,
-                       num_points=24)
+        foreground=n_air ** 2,
+        num_points=24,
+        )
 
     dxes = [grid.dxyz, grid.autoshifted_dxyz()]
     for a in (0, 1, 2):
@@ -71,9 +75,9 @@ def test0(solver=generic_solver):
     J[1][15, grid.shape[1]//2, grid.shape[2]//2] = 1
 
 
-    '''
-    Solve!
-    '''
+    #
+    # Solve!
+    #
     sim_args = {
         'omega': omega,
         'dxes': dxes,
@@ -87,9 +91,9 @@ def test0(solver=generic_solver):
 
     E = unvec(x, grid.shape)
 
-    '''
-    Plot results
-    '''
+    #
+    # Plot results
+    #
     pyplot.figure()
     pyplot.pcolor(numpy.real(E[1][:, :, grid.shape[2]//2]), cmap='seismic')
     pyplot.axis('equal')
@@ -122,7 +126,7 @@ def test1(solver=generic_solver):
     # #### Create the grid and draw the device ####
     grid = gridlock.Grid(edge_coords)
     epsilon = grid.allocate(n_air**2, dtype=numpy.float32)
-    grid.draw_cuboid(epsilon, center=center, dimensions=[8e3, w, th], eps=n_wg**2)
+    grid.draw_cuboid(epsilon, center=center, dimensions=[8e3, w, th], foreground=n_wg**2)
 
     dxes = [grid.dxyz, grid.autoshifted_dxyz()]
     for a in (0, 1, 2):
@@ -169,9 +173,9 @@ def test1(solver=generic_solver):
 #    pcolor((numpy.abs(J3).sum(axis=2).sum(axis=0) > 0).astype(float).T)
     pyplot.show(block=True)
 
-    '''
-    Solve!
-    '''
+    #
+    # Solve!
+    #
     sim_args = {
         'omega': omega,
         'dxes': dxes,
@@ -188,9 +192,9 @@ def test1(solver=generic_solver):
 
     E = unvec(x, grid.shape)
 
-    '''
-    Plot results
-    '''
+    #
+    # Plot results
+    #
     center = grid.pos2ind([0, 0, 0], None).astype(int)
     pyplot.figure()
     pyplot.subplot(2, 2, 1)
@@ -232,7 +236,7 @@ def test1(solver=generic_solver):
     pyplot.grid(alpha=0.6)
     pyplot.title('Overlap with mode')
     pyplot.show()
-    print('Average overlap with mode:', sum(q)/len(q))
+    print('Average overlap with mode:', sum(q[8:32])/len(q[8:32]))
 
 
 def module_available(name):

examples/fdtd.py

@@ -157,7 +157,8 @@ def main():
         e[1][tuple(grid.shape//2)] += field_source(t)
         update_H(e, h)
 
-        print('iteration {}: average {} iterations per sec'.format(t, (t+1)/(time.perf_counter()-start)))
+        avg_rate = (t + 1)/(time.perf_counter() - start))
+        print(f'iteration {t}: average {avg_rate} iterations per sec')
         sys.stdout.flush()
 
         if t % 20 == 0:

examples/tcyl.py

@@ -3,7 +3,7 @@ import numpy
 from numpy.linalg import norm
 
 from meanas.fdmath import vec, unvec
-from meanas.fdfd import waveguide_mode, functional, scpml
+from meanas.fdfd import waveguide_cyl, functional, scpml
 from meanas.fdfd.solvers import generic as generic_solver
 
 import gridlock
@@ -37,29 +37,34 @@ def test1(solver=generic_solver):
     xyz_max = numpy.array([800, y_max, z_max]) + (pml_thickness + 2) * dx
 
     # Coordinates of the edges of the cells.
-    half_edge_coords = [numpy.arange(dx/2, m + dx/2, step=dx) for m in xyz_max]
+    half_edge_coords = [numpy.arange(dx / 2, m + dx / 2, step=dx) for m in xyz_max]
     edge_coords = [numpy.hstack((-h[::-1], h)) for h in half_edge_coords]
     edge_coords[0] = numpy.array([-dx, dx])
 
     # #### Create the grid and draw the device ####
     grid = gridlock.Grid(edge_coords)
     epsilon = grid.allocate(n_air**2, dtype=numpy.float32)
-    grid.draw_cuboid(epsilon, center=center, dimensions=[8e3, w, th], eps=n_wg**2)
+    grid.draw_cuboid(epsilon, center=center, dimensions=[8e3, w, th], foreground=n_wg**2)
 
     dxes = [grid.dxyz, grid.autoshifted_dxyz()]
     for a in (1, 2):
         for p in (-1, 1):
-            dxes = scmpl.stretch_with_scpml(dxes, omega=omega, axis=a, polarity=p,
-                                            thickness=pml_thickness)
+            dxes = scpml.stretch_with_scpml(
+                dxes,
+                omega=omega,
+                axis=a,
+                polarity=p,
+                thickness=pml_thickness,
+                )
 
     wg_args = {
         'omega': omega,
         'dxes': [(d[1], d[2]) for d in dxes],
-        'epsilon': vec(g.transpose([1, 2, 0]) for g in epsilon),
+        'epsilon': vec(epsilon.transpose([0, 2, 3, 1])),
         'r0': r0,
     }
 
-    wg_results = waveguide_mode.solve_waveguide_mode_cylindrical(mode_number=0, **wg_args)
+    wg_results = waveguide_cyl.solve_mode(mode_number=0, **wg_args)
 
     E = wg_results['E']
 
@@ -70,20 +75,17 @@ def test1(solver=generic_solver):
     '''
     Plot results
     '''
-    def pcolor(v):
+    def pcolor(fig, ax, v, title):
         vmax = numpy.max(numpy.abs(v))
-        pyplot.pcolor(v.T, cmap='seismic', vmin=-vmax, vmax=vmax)
-        pyplot.axis('equal')
-        pyplot.colorbar()
+        mappable = ax.pcolormesh(v.T, cmap='seismic', vmin=-vmax, vmax=vmax)
+        ax.set_aspect('equal', adjustable='box')
+        ax.set_title(title)
+        ax.figure.colorbar(mappable)
 
-    pyplot.figure()
-    pyplot.subplot(2, 2, 1)
-    pcolor(numpy.real(E[0][:, :]))
-    pyplot.subplot(2, 2, 2)
-    pcolor(numpy.real(E[1][:, :]))
-    pyplot.subplot(2, 2, 3)
-    pcolor(numpy.real(E[2][:, :]))
-    pyplot.subplot(2, 2, 4)
+    fig, axes = pyplot.subplots(2, 2)
+    pcolor(fig, axes[0][0], numpy.real(E[0]), 'Ex')
+    pcolor(fig, axes[0][1], numpy.real(E[1]), 'Ey')
+    pcolor(fig, axes[1][0], numpy.real(E[2]), 'Ez')
     pyplot.show()
 
 

make_docs.sh

@@ -12,7 +12,7 @@ cd ~/projects/meanas
 rm -rf _doc_mathimg
 pdoc --pdf --force --template-dir pdoc_templates -o doc . > doc.md
 pandoc --metadata=title:"meanas" --from=markdown+abbreviations --to=html --output=doc.htex --gladtex -s --css pdoc_templates/pdoc.css doc.md
-gladtex -a -n -d _doc_mathimg -c white doc.htex
+gladtex -a -n -d _doc_mathimg -c white -b black doc.htex
 
 # Approach 3: html with gladtex
 #pdoc3 --html --force --template-dir pdoc_templates -o doc .

meanas/__init__.py

@@ -11,7 +11,8 @@ __author__ = 'Jan Petykiewicz'
 
 
 try:
-    with open(pathlib.Path(__file__).parent / 'README.md', 'r') as f:
+    readme_path = pathlib.Path(__file__).parent / 'README.md'
+    with readme_path.open('r') as f:
         __doc__ = f.read()
 except Exception:
     pass

meanas/eigensolvers.py

@@ -1,12 +1,12 @@
 """
 Solvers for eigenvalue / eigenvector problems
 """
-from typing import Callable
+from collections.abc import Callable
 import numpy
 from numpy.typing import NDArray, ArrayLike
 from numpy.linalg import norm
-from scipy import sparse              # type: ignore
-import scipy.sparse.linalg as spalg   # type: ignore
+from scipy import sparse
+import scipy.sparse.linalg as spalg
 
 
 def power_iteration(
@@ -25,8 +25,9 @@ def power_iteration(
     Returns:
         (Largest-magnitude eigenvalue, Corresponding eigenvector estimate)
     """
+    rng = numpy.random.default_rng()
     if guess_vector is None:
-        v = numpy.random.rand(operator.shape[0]) + 1j * numpy.random.rand(operator.shape[0])
+        v = rng.random(operator.shape[0]) + 1j * rng.random(operator.shape[0])
     else:
         v = guess_vector
 

meanas/fdfd/__init__.py

@@ -1,4 +1,4 @@
-"""
+r"""
 Tools for finite difference frequency-domain (FDFD) simulations and calculations.
 
 These mostly involve picking a single frequency, then setting up and solving a
@@ -19,71 +19,71 @@ Submodules:
 From the "Frequency domain" section of `meanas.fdmath`, we have
 
 $$
- \\begin{aligned}
- \\tilde{E}_{l, \\vec{r}} &= \\tilde{E}_{\\vec{r}} e^{-\\imath \\omega l \\Delta_t} \\\\
- \\tilde{H}_{l - \\frac{1}{2}, \\vec{r} + \\frac{1}{2}} &= \\tilde{H}_{\\vec{r} + \\frac{1}{2}} e^{-\\imath \\omega (l - \\frac{1}{2}) \\Delta_t} \\\\
- \\tilde{J}_{l, \\vec{r}} &= \\tilde{J}_{\\vec{r}} e^{-\\imath \\omega (l - \\frac{1}{2}) \\Delta_t} \\\\
- \\tilde{M}_{l - \\frac{1}{2}, \\vec{r} + \\frac{1}{2}} &= \\tilde{M}_{\\vec{r} + \\frac{1}{2}} e^{-\\imath \\omega l \\Delta_t} \\\\
- \\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{\\vec{r}})
-    -\\Omega^2 \\epsilon_{\\vec{r}} \\cdot \\tilde{E}_{\\vec{r}} &= -\\imath \\Omega \\tilde{J}_{\\vec{r}} e^{\\imath \\omega \\Delta_t / 2} \\\\
- \\Omega &= 2 \\sin(\\omega \\Delta_t / 2) / \\Delta_t
- \\end{aligned}
+ \begin{aligned}
+ \tilde{E}_{l, \vec{r}} &= \tilde{E}_{\vec{r}} e^{-\imath \omega l \Delta_t} \\
+ \tilde{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &= \tilde{H}_{\vec{r} + \frac{1}{2}} e^{-\imath \omega (l - \frac{1}{2}) \Delta_t} \\
+ \tilde{J}_{l, \vec{r}} &= \tilde{J}_{\vec{r}} e^{-\imath \omega (l - \frac{1}{2}) \Delta_t} \\
+ \tilde{M}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &= \tilde{M}_{\vec{r} + \frac{1}{2}} e^{-\imath \omega l \Delta_t} \\
+ \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{\vec{r}})
+    -\Omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} &= -\imath \Omega \tilde{J}_{\vec{r}} e^{\imath \omega \Delta_t / 2} \\
+ \Omega &= 2 \sin(\omega \Delta_t / 2) / \Delta_t
+ \end{aligned}
 $$
 
 resulting in
 
 $$
- \\begin{aligned}
- \\tilde{\\partial}_t &\\Rightarrow -\\imath \\Omega e^{-\\imath \\omega \\Delta_t / 2}\\\\
-   \\hat{\\partial}_t &\\Rightarrow -\\imath \\Omega e^{ \\imath \\omega \\Delta_t / 2}\\\\
- \\end{aligned}
+ \begin{aligned}
+ \tilde{\partial}_t &\Rightarrow -\imath \Omega e^{-\imath \omega \Delta_t / 2}\\
+   \hat{\partial}_t &\Rightarrow -\imath \Omega e^{ \imath \omega \Delta_t / 2}\\
+ \end{aligned}
 $$
 
 Maxwell's equations are then
 
 $$
-  \\begin{aligned}
-  \\tilde{\\nabla} \\times \\tilde{E}_{\\vec{r}} &=
-         \\imath \\Omega e^{-\\imath \\omega \\Delta_t / 2} \\hat{B}_{\\vec{r} + \\frac{1}{2}}
-                                                          - \\hat{M}_{\\vec{r} + \\frac{1}{2}}  \\\\
-  \\hat{\\nabla} \\times \\hat{H}_{\\vec{r} + \\frac{1}{2}} &=
-        -\\imath \\Omega e^{ \\imath \\omega \\Delta_t / 2} \\tilde{D}_{\\vec{r}}
-                                                          + \\tilde{J}_{\\vec{r}} \\\\
-  \\tilde{\\nabla} \\cdot \\hat{B}_{\\vec{r} + \\frac{1}{2}} &= 0 \\\\
-  \\hat{\\nabla} \\cdot \\tilde{D}_{\\vec{r}} &= \\rho_{\\vec{r}}
- \\end{aligned}
+  \begin{aligned}
+  \tilde{\nabla} \times \tilde{E}_{\vec{r}} &=
+         \imath \Omega e^{-\imath \omega \Delta_t / 2} \hat{B}_{\vec{r} + \frac{1}{2}}
+                                                     - \hat{M}_{\vec{r} + \frac{1}{2}}  \\
+  \hat{\nabla} \times \hat{H}_{\vec{r} + \frac{1}{2}} &=
+        -\imath \Omega e^{ \imath \omega \Delta_t / 2} \tilde{D}_{\vec{r}}
+                                                     + \tilde{J}_{\vec{r}} \\
+  \tilde{\nabla} \cdot \hat{B}_{\vec{r} + \frac{1}{2}} &= 0 \\
+  \hat{\nabla} \cdot \tilde{D}_{\vec{r}} &= \rho_{\vec{r}}
+ \end{aligned}
 $$
 
-With $\\Delta_t \\to 0$, this simplifies to
+With $\Delta_t \to 0$, this simplifies to
 
 $$
- \\begin{aligned}
- \\tilde{E}_{l, \\vec{r}} &\\to \\tilde{E}_{\\vec{r}} \\\\
- \\tilde{H}_{l - \\frac{1}{2}, \\vec{r} + \\frac{1}{2}} &\\to \\tilde{H}_{\\vec{r} + \\frac{1}{2}} \\\\
- \\tilde{J}_{l, \\vec{r}} &\\to \\tilde{J}_{\\vec{r}} \\\\
- \\tilde{M}_{l - \\frac{1}{2}, \\vec{r} + \\frac{1}{2}} &\\to \\tilde{M}_{\\vec{r} + \\frac{1}{2}} \\\\
- \\Omega &\\to \\omega \\\\
- \\tilde{\\partial}_t &\\to -\\imath \\omega \\\\
-   \\hat{\\partial}_t &\\to -\\imath \\omega \\\\
- \\end{aligned}
+ \begin{aligned}
+ \tilde{E}_{l, \vec{r}} &\to \tilde{E}_{\vec{r}} \\
+ \tilde{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &\to \tilde{H}_{\vec{r} + \frac{1}{2}} \\
+ \tilde{J}_{l, \vec{r}} &\to \tilde{J}_{\vec{r}} \\
+ \tilde{M}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} &\to \tilde{M}_{\vec{r} + \frac{1}{2}} \\
+ \Omega &\to \omega \\
+ \tilde{\partial}_t &\to -\imath \omega \\
+   \hat{\partial}_t &\to -\imath \omega \\
+ \end{aligned}
 $$
 
 and then
 
 $$
-  \\begin{aligned}
-  \\tilde{\\nabla} \\times \\tilde{E}_{\\vec{r}} &=
-         \\imath \\omega \\hat{B}_{\\vec{r} + \\frac{1}{2}}
-                       - \\hat{M}_{\\vec{r} + \\frac{1}{2}}  \\\\
-  \\hat{\\nabla} \\times \\hat{H}_{\\vec{r} + \\frac{1}{2}} &=
-        -\\imath \\omega \\tilde{D}_{\\vec{r}}
-                       + \\tilde{J}_{\\vec{r}} \\\\
- \\end{aligned}
+  \begin{aligned}
+  \tilde{\nabla} \times \tilde{E}_{\vec{r}} &=
+         \imath \omega \hat{B}_{\vec{r} + \frac{1}{2}}
+                     - \hat{M}_{\vec{r} + \frac{1}{2}}  \\
+  \hat{\nabla} \times \hat{H}_{\vec{r} + \frac{1}{2}} &=
+        -\imath \omega \tilde{D}_{\vec{r}}
+                     + \tilde{J}_{\vec{r}} \\
+ \end{aligned}
 $$
 
 $$
- \\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{\\vec{r}})
-    -\\omega^2 \\epsilon_{\\vec{r}} \\cdot \\tilde{E}_{\\vec{r}} = -\\imath \\omega \\tilde{J}_{\\vec{r}} \\\\
+ \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{\vec{r}})
+    -\omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} = -\imath \omega \tilde{J}_{\vec{r}} \\
 $$
 
 # TODO FDFD?
@@ -91,5 +91,12 @@ $$
 
 
 """
-from . import solvers, operators, functional, scpml, waveguide_2d, waveguide_3d
+from . import (
+    solvers as solvers,
+    operators as operators,
+    functional as functional,
+    scpml as scpml,
+    waveguide_2d as waveguide_2d,
+    waveguide_3d as waveguide_3d,
+    )
 # from . import farfield, bloch TODO

meanas/fdfd/bloch.py

@@ -94,16 +94,17 @@ This module contains functions for generating and solving the
 
 """
 
-from typing import Callable, Any, cast, Sequence
+from typing import Any, cast
+from collections.abc import Callable, Sequence
 import logging
 import numpy
 from numpy import pi, real, trace
 from numpy.fft import fftfreq
 from numpy.typing import NDArray, ArrayLike
-import scipy                            # type: ignore
-import scipy.optimize                   # type: ignore
-from scipy.linalg import norm           # type: ignore
-import scipy.sparse.linalg as spalg     # type: ignore
+import scipy
+import scipy.optimize
+from scipy.linalg import norm
+import scipy.sparse.linalg as spalg
 
 from ..fdmath import fdfield_t, cfdfield_t
 
@@ -114,7 +115,6 @@ logger = logging.getLogger(__name__)
 try:
     import pyfftw.interfaces.numpy_fft  # type: ignore
     import pyfftw.interfaces            # type: ignore
-    import multiprocessing
     logger.info('Using pyfftw')
 
     pyfftw.interfaces.cache.enable()
@@ -155,7 +155,7 @@ def generate_kmn(
             All are given in the xyz basis (e.g. `|k|[0,0,0] = norm(G_matrix @ k0)`).
     """
     k0 = numpy.array(k0)
-    G_matrix = numpy.array(G_matrix, copy=False)
+    G_matrix = numpy.asarray(G_matrix)
 
     Gi_grids = numpy.array(numpy.meshgrid(*(fftfreq(n, 1 / n) for n in shape[:3]), indexing='ij'))
     Gi = numpy.moveaxis(Gi_grids, 0, -1)
@@ -232,7 +232,7 @@ def maxwell_operator(
             Raveled conv(1/mu_k, ik x conv(1/eps_k, ik x h_mn)), returned
             and overwritten in-place of `h`.
         """
-        hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
+        hin_m, hin_n = (hi.reshape(shape) for hi in numpy.split(h, 2))
 
         #{d,e,h}_xyz fields are complex 3-fields in (1/x, 1/y, 1/z) basis
 
@@ -303,12 +303,12 @@ def hmn_2_exyz(
     k_mag, m, n = generate_kmn(k0, G_matrix, shape)
 
     def operator(h: NDArray[numpy.complex128]) -> cfdfield_t:
-        hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
+        hin_m, hin_n = (hi.reshape(shape) for hi in numpy.split(h, 2))
         d_xyz = (n * hin_m
                - m * hin_n) * k_mag         # noqa: E128
 
         # divide by epsilon
-        return numpy.array([ei for ei in numpy.moveaxis(ifftn(d_xyz, axes=range(3)) / epsilon, 3, 0)])         # TODO avoid copy
+        return numpy.moveaxis(ifftn(d_xyz, axes=range(3)) / epsilon, 3, 0)
 
     return operator
 
@@ -341,7 +341,7 @@ def hmn_2_hxyz(
     _k_mag, m, n = generate_kmn(k0, G_matrix, shape)
 
     def operator(h: NDArray[numpy.complex128]) -> cfdfield_t:
-        hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
+        hin_m, hin_n = (hi.reshape(shape) for hi in numpy.split(h, 2))
         h_xyz = (m * hin_m
                + n * hin_n)     # noqa: E128
         return numpy.array([ifftn(hi) for hi in numpy.moveaxis(h_xyz, 3, 0)])
@@ -394,7 +394,7 @@ def inverse_maxwell_operator_approx(
         Returns:
             Raveled ik x conv(eps_k, ik x conv(mu_k, h_mn))
         """
-        hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)]
+        hin_m, hin_n = (hi.reshape(shape) for hi in numpy.split(h, 2))
 
         #{d,e,h}_xyz fields are complex 3-fields in (1/x, 1/y, 1/z) basis
 
@@ -538,7 +538,7 @@ def eigsolve(
         `(eigenvalues, eigenvectors)` where `eigenvalues[i]` corresponds to the
         vector `eigenvectors[i, :]`
     """
-    k0 = numpy.array(k0, copy=False)
+    k0 = numpy.asarray(k0)
 
     h_size = 2 * epsilon[0].size
 
@@ -561,11 +561,12 @@ def eigsolve(
     prev_theta = 0.5
     D = numpy.zeros(shape=y_shape, dtype=complex)
 
+    rng = numpy.random.default_rng()
     Z: NDArray[numpy.complex128]
     if y0 is None:
-        Z = numpy.random.rand(*y_shape) + 1j * numpy.random.rand(*y_shape)
+        Z = rng.random(y_shape) + 1j * rng.random(y_shape)
     else:
-        Z = numpy.array(y0, copy=False).T
+        Z = numpy.asarray(y0).T
 
     while True:
         Z *= num_modes / norm(Z)
@@ -573,7 +574,7 @@ def eigsolve(
         try:
             U = numpy.linalg.inv(ZtZ)
         except numpy.linalg.LinAlgError:
-            Z = numpy.random.rand(*y_shape) + 1j * numpy.random.rand(*y_shape)
+            Z = rng.random(y_shape) + 1j * rng.random(y_shape)
             continue
 
         trace_U = real(trace(U))
@@ -646,8 +647,7 @@ def eigsolve(
 
         Qi_memo: list[float | None] = [None, None]
 
-        def Qi_func(theta: float) -> float:
-            nonlocal Qi_memo
+        def Qi_func(theta: float, Qi_memo=Qi_memo, ZtZ=ZtZ, DtD=DtD, symZtD=symZtD) -> float:   # noqa: ANN001
             if Qi_memo[0] == theta:
                 return cast(float, Qi_memo[1])
 
@@ -656,7 +656,7 @@ def eigsolve(
             Q = c * c * ZtZ + s * s * DtD + 2 * s * c * symZtD
             try:
                 Qi = numpy.linalg.inv(Q)
-            except numpy.linalg.LinAlgError:
+            except numpy.linalg.LinAlgError as err:
                 logger.info('taylor Qi')
                 # if c or s small, taylor expand
                 if c < 1e-4 * s and c != 0:
@@ -666,12 +666,12 @@ def eigsolve(
                     ZtZi = numpy.linalg.inv(ZtZ)
                     Qi = ZtZi / (c * c) - 2 * s / (c * c * c) * (ZtZi @ (ZtZi @ symZtD).conj().T)
                 else:
-                    raise Exception('Inexplicable singularity in trace_func')
+                    raise Exception('Inexplicable singularity in trace_func') from err
             Qi_memo[0] = theta
             Qi_memo[1] = cast(float, Qi)
             return cast(float, Qi)
 
-        def trace_func(theta: float) -> float:
+        def trace_func(theta: float, ZtAZ=ZtAZ, DtAD=DtAD, symZtAD=symZtAD) -> float:           # noqa: ANN001
             c = numpy.cos(theta)
             s = numpy.sin(theta)
             Qi = Qi_func(theta)
@@ -680,15 +680,15 @@ def eigsolve(
             return numpy.abs(trace)
 
         if False:
-            def trace_deriv(theta):
+            def trace_deriv(theta, sgn: int = sgn, ZtAZ=ZtAZ, DtAD=DtAD, symZtD=symZtD, symZtAD=symZtAD, ZtZ=ZtZ, DtD=DtD):     # noqa: ANN001
                 Qi = Qi_func(theta)
                 c2 = numpy.cos(2 * theta)
                 s2 = numpy.sin(2 * theta)
-                F = -0.5*s2 *  (ZtAZ - DtAD) + c2 * symZtAD
+                F = -0.5 * s2 * (ZtAZ - DtAD) + c2 * symZtAD
                 trace_deriv = _rtrace_AtB(Qi, F)
 
                 G = Qi @ F.conj().T @ Qi.conj().T
-                H = -0.5*s2 * (ZtZ - DtD) + c2 * symZtD
+                H = -0.5 * s2 * (ZtZ - DtD) + c2 * symZtD
                 trace_deriv -= _rtrace_AtB(G, H)
 
                 trace_deriv *= 2
@@ -696,12 +696,12 @@ def eigsolve(
 
             U_sZtD = U @ symZtD
 
-            dE = 2.0 * (_rtrace_AtB(U, symZtAD) -
-                        _rtrace_AtB(ZtAZU, U_sZtD))
+            dE = 2.0 * (_rtrace_AtB(U, symZtAD)
+                        - _rtrace_AtB(ZtAZU, U_sZtD))
 
-            d2E = 2 * (_rtrace_AtB(U, DtAD) -
-                       _rtrace_AtB(ZtAZU, U @ (DtD - 4 * symZtD @ U_sZtD)) -
-                   4 * _rtrace_AtB(U, symZtAD @ U_sZtD))
+            d2E = 2 * (_rtrace_AtB(U, DtAD)
+                       - _rtrace_AtB(ZtAZU, U @ (DtD - 4 * symZtD @ U_sZtD))
+                       - 4 * _rtrace_AtB(U, symZtAD @ U_sZtD))
 
             # Newton-Raphson to find a root of the first derivative:
             theta = -dE / d2E
@@ -781,7 +781,7 @@ def linmin(x_guess, f0, df0, x_max, f_tol=0.1, df_tol=min(tolerance, 1e-6), x_to
                                           x_min, x_max, isave, dsave)
         for i in range(int(1e6)):
             if task != 'F':
-                logging.info('search converged in {} iterations'.format(i))
+                logging.info(f'search converged in {i} iterations')
                 break
             fx = f(x, dfx)
             x, fx, dfx, task = minpack2.dsrch(x, fx, dfx, f_tol, df_tol, x_tol, task,
@@ -799,3 +799,52 @@ def _rtrace_AtB(
 def _symmetrize(A: NDArray[numpy.complex128]) -> NDArray[numpy.complex128]:
     return (A + A.conj().T) * 0.5
 
+
+
+def inner_product(eL, hL, eR, hR) -> complex:
+    # assumes x-axis propagation
+
+    assert numpy.array_equal(eR.shape, hR.shape)
+    assert numpy.array_equal(eL.shape, hL.shape)
+    assert numpy.array_equal(eR.shape, eL.shape)
+
+    # Cross product, times 2 since it's <p | n>, then divide by 4. # TODO might want to abs() this?
+    norm2R = (eR[1] * hR[2] - eR[2] * hR[1]).sum() / 2
+    norm2L = (eL[1] * hL[2] - eL[2] * hL[1]).sum() / 2
+
+    # eRxhR_x = numpy.cross(eR.reshape(3, -1), hR.reshape(3, -1), axis=0).reshape(eR.shape)[0] / normR
+    # logger.info(f'power {eRxhR_x.sum() / 2})
+
+    eR /= numpy.sqrt(norm2R)
+    hR /= numpy.sqrt(norm2R)
+    eL /= numpy.sqrt(norm2L)
+    hL /= numpy.sqrt(norm2L)
+
+    # (eR x hL)[0] and (eL x hR)[0]
+    eRxhL_x = eR[1] * hL[2] - eR[2] - hL[1]
+    eLxhR_x = eL[1] * hR[2] - eL[2] - hR[1]
+
+    #return 1j *  (eRxhL_x - eLxhR_x).sum() / numpy.sqrt(norm2R * norm2L)
+    #return (eRxhL_x.sum() - eLxhR_x.sum()) / numpy.sqrt(norm2R * norm2L)
+    return eRxhL_x.sum() - eLxhR_x.sum()
+
+
+def trq(eI, hI, eO, hO) -> tuple[complex, complex]:
+    pp = inner_product(eO,  hO, eI,  hI)
+    pn = inner_product(eO,  hO, eI, -hI)
+    np = inner_product(eO, -hO, eI,  hI)
+    nn = inner_product(eO, -hO, eI, -hI)
+
+    assert pp == -nn
+    assert pn == -np
+
+    logger.info(f'''
+        {pp=:4g} {pn=:4g}
+        {nn=:4g} {np=:4g}
+            {nn * pp / pn=:4g}    {-np=:4g}
+        ''')
+
+    r = -pp / pn    # -<Pp|Bp>/<Pn/Bp> = -(-pp) / (-pn)
+    t = (np - nn * pp / pn) / 4
+
+    return t, r

meanas/fdfd/eme.py

@@ -0,0 +1,68 @@
+import numpy
+
+from ..fdmath import vec, unvec, dx_lists_t, vfdfield_t, vcfdfield_t
+from .waveguide_2d import inner_product
+
+
+def get_tr(ehL, wavenumbers_L, ehR, wavenumbers_R, dxes: dx_lists_t):
+    nL = len(wavenumbers_L)
+    nR = len(wavenumbers_R)
+    A12 = numpy.zeros((nL, nR), dtype=complex)
+    A21 = numpy.zeros((nL, nR), dtype=complex)
+    B11 = numpy.zeros((nL,), dtype=complex)
+    for ll in range(nL):
+        eL, hL = ehL[ll]
+        B11[ll] = inner_product(eL, hL, dxes=dxes, conj_h=False)
+        for rr in range(nR):
+            eR, hR = ehR[rr]
+            A12[ll, rr] = inner_product(eL, hR, dxes=dxes, conj_h=False)    # TODO optimize loop?
+            A21[ll, rr] = inner_product(eR, hL, dxes=dxes, conj_h=False)
+
+    # tt0 = 2 * numpy.linalg.pinv(A21 + numpy.conj(A12))
+    tt0, _resid, _rank, _sing = numpy.linalg.lstsq(A21 + A12, numpy.diag(2 * B11), rcond=None)
+
+    U, st, V = numpy.linalg.svd(tt0)
+    gain = st > 1
+    st[gain] = 1 / st[gain]
+    tt = U @ numpy.diag(st) @ V
+
+    # rr = 0.5 * (A21 - numpy.conj(A12)) @ tt
+    rr = numpy.diag(0.5 / B11) @ (A21 - A12) @ tt
+
+    return tt, rr
+
+
+def get_abcd(eL_xys, wavenumbers_L, eR_xys, wavenumbers_R, **kwargs):
+    t12, r12 = get_tr(eL_xys, wavenumbers_L, eR_xys, wavenumbers_R, **kwargs)
+    t21, r21 = get_tr(eR_xys, wavenumbers_R, eL_xys, wavenumbers_L, **kwargs)
+    t21i = numpy.linalg.pinv(t21)
+    A = t12 - r21 @ t21i @ r12
+    B = r21 @ t21i
+    C = -t21i @ r12
+    D = t21i
+    return sparse.block_array(((A, B), (C, D)))
+
+
+def get_s(
+        eL_xys,
+        wavenumbers_L,
+        eR_xys,
+        wavenumbers_R,
+        force_nogain: bool = False,
+        force_reciprocal: bool = False,
+        **kwargs):
+    t12, r12 = get_tr(eL_xys, wavenumbers_L, eR_xys, wavenumbers_R, **kwargs)
+    t21, r21 = get_tr(eR_xys, wavenumbers_R, eL_xys, wavenumbers_L, **kwargs)
+
+    ss = numpy.block([[r12, t12],
+                      [t21, r21]])
+
+    if force_nogain:
+        # force S @ S.H diagonal
+        U, sing, V = numpy.linalg.svd(ss)
+        ss = numpy.diag(sing) @ U @ V
+
+    if force_reciprocal:
+        ss = 0.5 * (ss + ss.T)
+
+    return ss

meanas/fdfd/farfield.py

@@ -1,7 +1,8 @@
 """
 Functions for performing near-to-farfield transformation (and the reverse).
 """
-from typing import Any, Sequence, cast
+from typing import Any, cast
+from collections.abc import Sequence
 import numpy
 from numpy.fft import fft2, fftshift, fftfreq, ifft2, ifftshift
 from numpy import pi

meanas/fdfd/functional.py

@@ -5,7 +5,7 @@ Functional versions of many FDFD operators. These can be useful for performing
 The functions generated here expect `cfdfield_t` inputs with shape (3, X, Y, Z),
 e.g. E = [E_x, E_y, E_z] where each (complex) component has shape (X, Y, Z)
 """
-from typing import Callable
+from collections.abc import Callable
 import numpy
 
 from ..fdmath import dx_lists_t, fdfield_t, cfdfield_t, cfdfield_updater_t
@@ -47,7 +47,6 @@ def e_full(
 
     if mu is None:
         return op_1
-    else:
     return op_mu
 
 
@@ -84,7 +83,6 @@ def eh_full(
 
     if mu is None:
         return op_1
-    else:
     return op_mu
 
 
@@ -116,7 +114,6 @@ def e2h(
 
     if mu is None:
         return e2h_1_1
-    else:
     return e2h_mu
 
 
@@ -151,7 +148,6 @@ def m2j(
 
     if mu is None:
         return m2j_1
-    else:
     return m2j_mu
 
 
@@ -189,10 +185,10 @@ def e_tfsf_source(
 
 
 def poynting_e_cross_h(dxes: dx_lists_t) -> Callable[[cfdfield_t, cfdfield_t], cfdfield_t]:
-    """
+    r"""
     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$
+    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

meanas/fdfd/operators.py

@@ -28,7 +28,7 @@ The following operators are included:
 """
 
 import numpy
-import scipy.sparse as sparse       # type: ignore
+from scipy import sparse
 
 from ..fdmath import vec, dx_lists_t, vfdfield_t, vcfdfield_t
 from ..fdmath.operators import shift_with_mirror, shift_circ, curl_forward, curl_back
@@ -40,19 +40,19 @@ __author__ = 'Jan Petykiewicz'
 def e_full(
         omega: complex,
         dxes: dx_lists_t,
-        epsilon: vfdfield_t,
+        epsilon: vfdfield_t | vcfdfield_t,
         mu: vfdfield_t | None = None,
         pec: vfdfield_t | None = None,
         pmc: vfdfield_t | None = None,
         ) -> sparse.spmatrix:
-    """
+    r"""
     Wave operator
-     $$ \\nabla \\times (\\frac{1}{\\mu} \\nabla \\times) - \\Omega^2 \\epsilon $$
+     $$ \nabla \times (\frac{1}{\mu} \nabla \times) - \Omega^2 \epsilon $$
 
         del x (1/mu * del x) - omega**2 * epsilon
 
      for use with the E-field, with wave equation
-     $$ (\\nabla \\times (\\frac{1}{\\mu} \\nabla \\times) - \\Omega^2 \\epsilon) E = -\\imath \\omega J $$
+     $$ (\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
 
@@ -131,14 +131,14 @@ def h_full(
         pec: vfdfield_t | None = None,
         pmc: vfdfield_t | None = None,
         ) -> sparse.spmatrix:
-    """
+    r"""
     Wave operator
-     $$ \\nabla \\times (\\frac{1}{\\epsilon} \\nabla \\times) - \\omega^2 \\mu $$
+     $$ \nabla \times (\frac{1}{\epsilon} \nabla \times) - \omega^2 \mu $$
 
         del x (1/epsilon * del x) - omega**2 * mu
 
      for use with the H-field, with wave equation
-     $$ (\\nabla \\times (\\frac{1}{\\epsilon} \\nabla \\times) - \\omega^2 \\mu) E = \\imath \\omega M $$
+     $$ (\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
 
@@ -188,28 +188,28 @@ def eh_full(
         pec: vfdfield_t | None = None,
         pmc: vfdfield_t | None = None,
         ) -> sparse.spmatrix:
-    """
+    r"""
     Wave operator for `[E, H]` field representation. This operator implements Maxwell's
      equations without cancelling out either E or H. The operator is
-    $$  \\begin{bmatrix}
-        -\\imath \\omega \\epsilon  &  \\nabla \\times      \\\\
-        \\nabla \\times             &  \\imath \\omega \\mu
-        \\end{bmatrix} $$
+    $$  \begin{bmatrix}
+        -\imath \omega \epsilon  &  \nabla \times      \\
+        \nabla \times            &  \imath \omega \mu
+        \end{bmatrix} $$
 
         [[-i * omega * epsilon,  del x         ],
          [del x,                 i * omega * mu]]
 
     for use with a field vector of the form `cat(vec(E), vec(H))`:
-    $$  \\begin{bmatrix}
-        -\\imath \\omega \\epsilon  &  \\nabla \\times      \\\\
-        \\nabla \\times             &  \\imath \\omega \\mu
-        \\end{bmatrix}
-        \\begin{bmatrix} E \\\\
+    $$  \begin{bmatrix}
+        -\imath \omega \epsilon  &  \nabla \times      \\
+        \nabla \times            &  \imath \omega \mu
+        \end{bmatrix}
+        \begin{bmatrix} E \\
                         H
-        \\end{bmatrix}
-        = \\begin{bmatrix} J \\\\
+        \end{bmatrix}
+        = \begin{bmatrix} J \\
                          -M
-          \\end{bmatrix} $$
+          \end{bmatrix} $$
 
     Args:
         omega: Angular frequency of the simulation
@@ -321,11 +321,11 @@ def poynting_e_cross(e: vcfdfield_t, dxes: dx_lists_t) -> sparse.spmatrix:
     """
     shape = [len(dx) for dx in dxes[0]]
 
-    fx, fy, fz = [shift_circ(i, shape, 1) for i in range(3)]
+    fx, fy, fz = (shift_circ(i, shape, 1) for i in range(3))
 
     dxag = [dx.ravel(order='C') for dx in numpy.meshgrid(*dxes[0], indexing='ij')]
     dxbg = [dx.ravel(order='C') for dx in numpy.meshgrid(*dxes[1], indexing='ij')]
-    Ex, Ey, Ez = [ei * da for ei, da in zip(numpy.split(e, 3), dxag)]
+    Ex, Ey, Ez = (ei * da for ei, da in zip(numpy.split(e, 3), dxag, strict=True))
 
     block_diags = [[ None,     fx @ -Ez, fx @  Ey],
                    [ fy @  Ez, None,     fy @ -Ex],
@@ -349,11 +349,11 @@ def poynting_h_cross(h: vcfdfield_t, dxes: dx_lists_t) -> sparse.spmatrix:
     """
     shape = [len(dx) for dx in dxes[0]]
 
-    fx, fy, fz = [shift_circ(i, shape, 1) for i in range(3)]
+    fx, fy, fz = (shift_circ(i, shape, 1) for i in range(3))
 
     dxag = [dx.ravel(order='C') for dx in numpy.meshgrid(*dxes[0], indexing='ij')]
     dxbg = [dx.ravel(order='C') for dx in numpy.meshgrid(*dxes[1], indexing='ij')]
-    Hx, Hy, Hz = [sparse.diags(hi * db) for hi, db in zip(numpy.split(h, 3), dxbg)]
+    Hx, Hy, Hz = (sparse.diags(hi * db) for hi, db in zip(numpy.split(h, 3), dxbg, strict=True))
 
     P = (sparse.bmat(
         [[ None,    -Hz @ fx,   Hy @ fx],

meanas/fdfd/scpml.py

@@ -2,7 +2,7 @@
 Functions for creating stretched coordinate perfectly matched layer (PML) absorbers.
 """
 
-from typing import Sequence, Callable
+from collections.abc import Sequence, Callable
 
 import numpy
 from numpy.typing import NDArray

meanas/fdfd/solvers.py

@@ -2,13 +2,14 @@
 Solvers and solver interface for FDFD problems.
 """
 
-from typing import Callable, Dict, Any, Optional
+from typing import Any
+from collections.abc import Callable
 import logging
 
 import numpy
 from numpy.typing import ArrayLike, NDArray
 from numpy.linalg import norm
-import scipy.sparse.linalg          # type: ignore
+import scipy.sparse.linalg
 
 from ..fdmath import dx_lists_t, vfdfield_t, vcfdfield_t
 from . import operators
@@ -34,16 +35,17 @@ def _scipy_qmr(
         Guess for solution (returned even if didn't converge)
     """
 
-    '''
-    Report on our progress
-    '''
+    #
+    #Report on our progress
+    #
     ii = 0
 
     def log_residual(xk: ArrayLike) -> None:
         nonlocal ii
         ii += 1
         if ii % 100 == 0:
-            logger.info('Solver residual at iteration {} : {}'.format(ii, norm(A @ xk - b)))
+            cur_norm = norm(A @ xk - b) / norm(b)
+            logger.info(f'Solver residual at iteration {ii} : {cur_norm}')
 
     if 'callback' in kwargs:
         def augmented_callback(xk: ArrayLike) -> None:
@@ -54,10 +56,9 @@ def _scipy_qmr(
     else:
         kwargs['callback'] = log_residual
 
-    '''
-    Run the actual solve
-    '''
-
+    #
+    # Run the actual solve
+    #
     x, _ = scipy.sparse.linalg.qmr(A, b, **kwargs)
     return x
 
@@ -67,12 +68,14 @@ def generic(
         dxes: dx_lists_t,
         J: vcfdfield_t,
         epsilon: vfdfield_t,
-        mu: Optional[vfdfield_t] = None,
-        pec: Optional[vfdfield_t] = None,
-        pmc: Optional[vfdfield_t] = None,
+        mu: vfdfield_t | None = None,
+        *,
+        pec: vfdfield_t | None = None,
+        pmc: vfdfield_t | None = None,
         adjoint: bool = False,
         matrix_solver: Callable[..., ArrayLike] = _scipy_qmr,
-        matrix_solver_opts: Optional[Dict[str, Any]] = None,
+        matrix_solver_opts: dict[str, Any] | None = None,
+        E_guess: vcfdfield_t | None = None,
         ) -> vcfdfield_t:
     """
     Conjugate gradient FDFD solver using CSR sparse matrices.
@@ -99,6 +102,8 @@ def generic(
                  which doesn't return convergence info and logs the residual
                  every 100 iterations.
         matrix_solver_opts: Passed as kwargs to `matrix_solver(...)`
+        E_guess: Guess at the solution E-field. `matrix_solver` must accept an
+                `x0` argument with the same purpose.
 
     Returns:
         E-field which solves the system.
@@ -119,6 +124,13 @@ def generic(
         A = Pl @ A0 @ Pr
         b = Pl @ b0
 
+    if E_guess is not None:
+        if adjoint:
+            x0 = Pr.H @ E_guess
+        else:
+            x0 = Pl @ E_guess
+        matrix_solver_opts['x0'] = x0
+
     x = matrix_solver(A.tocsr(), b, **matrix_solver_opts)
 
     if adjoint:

meanas/fdfd/waveguide_2d.py

@@ -1,4 +1,4 @@
-"""
+r"""
 Operators and helper functions for waveguides with unchanging cross-section.
 
 The propagation direction is chosen to be along the z axis, and all fields
@@ -12,180 +12,180 @@ As the z-dependence is known, all the functions in this file assume a 2D grid
 
 Consider Maxwell's equations in continuous space, in the frequency domain. Assuming
 a structure with some (x, y) cross-section extending uniformly into the z dimension,
-with a diagonal $\\epsilon$ tensor, we have
+with a diagonal $\epsilon$ tensor, we have
 
 $$
-\\begin{aligned}
-\\nabla \\times \\vec{E}(x, y, z) &= -\\imath \\omega \\mu \\vec{H} \\\\
-\\nabla \\times \\vec{H}(x, y, z) &=  \\imath \\omega \\epsilon \\vec{E} \\\\
-\\vec{E}(x,y,z) = (\\vec{E}_t(x, y) + E_z(x, y)\\vec{z}) e^{-\\gamma z} \\\\
-\\vec{H}(x,y,z) = (\\vec{H}_t(x, y) + H_z(x, y)\\vec{z}) e^{-\\gamma z} \\\\
-\\end{aligned}
+\begin{aligned}
+\nabla \times \vec{E}(x, y, z) &= -\imath \omega \mu \vec{H} \\
+\nabla \times \vec{H}(x, y, z) &=  \imath \omega \epsilon \vec{E} \\
+\vec{E}(x,y,z) &= (\vec{E}_t(x, y) + E_z(x, y)\vec{z}) e^{-\imath \beta z} \\
+\vec{H}(x,y,z) &= (\vec{H}_t(x, y) + H_z(x, y)\vec{z}) e^{-\imath \beta z} \\
+\end{aligned}
 $$
 
 Expanding the first two equations into vector components, we get
 
 $$
-\\begin{aligned}
--\\imath \\omega \\mu_{xx} H_x &= \\partial_y E_z - \\partial_z E_y \\\\
--\\imath \\omega \\mu_{yy} H_y &= \\partial_z E_x - \\partial_x E_z \\\\
--\\imath \\omega \\mu_{zz} H_z &= \\partial_x E_y - \\partial_y E_x \\\\
-\\imath \\omega \\epsilon_{xx} E_x &= \\partial_y H_z - \\partial_z H_y \\\\
-\\imath \\omega \\epsilon_{yy} E_y &= \\partial_z H_x - \\partial_x H_z \\\\
-\\imath \\omega \\epsilon_{zz} E_z &= \\partial_x H_y - \\partial_y H_x \\\\
-\\end{aligned}
+\begin{aligned}
+-\imath \omega \mu_{xx} H_x &= \partial_y E_z - \partial_z E_y \\
+-\imath \omega \mu_{yy} H_y &= \partial_z E_x - \partial_x E_z \\
+-\imath \omega \mu_{zz} H_z &= \partial_x E_y - \partial_y E_x \\
+\imath \omega \epsilon_{xx} E_x &= \partial_y H_z - \partial_z H_y \\
+\imath \omega \epsilon_{yy} E_y &= \partial_z H_x - \partial_x H_z \\
+\imath \omega \epsilon_{zz} E_z &= \partial_x H_y - \partial_y H_x \\
+\end{aligned}
 $$
 
-Substituting in our expressions for $\\vec{E}$, $\\vec{H}$ and discretizing:
+Substituting in our expressions for $\vec{E}$, $\vec{H}$ and discretizing:
 
 $$
-\\begin{aligned}
--\\imath \\omega \\mu_{xx} H_x &= \\tilde{\\partial}_y E_z + \\gamma E_y \\\\
--\\imath \\omega \\mu_{yy} H_y &= -\\gamma E_x - \\tilde{\\partial}_x E_z \\\\
--\\imath \\omega \\mu_{zz} H_z &= \\tilde{\\partial}_x E_y - \\tilde{\\partial}_y E_x \\\\
-\\imath \\omega \\epsilon_{xx} E_x &= \\hat{\\partial}_y H_z + \\gamma H_y \\\\
-\\imath \\omega \\epsilon_{yy} E_y &= -\\gamma H_x - \\hat{\\partial}_x H_z \\\\
-\\imath \\omega \\epsilon_{zz} E_z &= \\hat{\\partial}_x H_y - \\hat{\\partial}_y H_x \\\\
-\\end{aligned}
+\begin{aligned}
+-\imath \omega \mu_{xx} H_x &= \tilde{\partial}_y E_z + \imath \beta E_y \\
+-\imath \omega \mu_{yy} H_y &= -\imath \beta E_x - \tilde{\partial}_x E_z \\
+-\imath \omega \mu_{zz} H_z &= \tilde{\partial}_x E_y - \tilde{\partial}_y E_x \\
+\imath \omega \epsilon_{xx} E_x &= \hat{\partial}_y H_z + \imath \beta H_y \\
+\imath \omega \epsilon_{yy} E_y &= -\imath \beta H_x - \hat{\partial}_x H_z \\
+\imath \omega \epsilon_{zz} E_z &= \hat{\partial}_x H_y - \hat{\partial}_y H_x \\
+\end{aligned}
 $$
 
 Rewrite the last three equations as
-$$
-\\begin{aligned}
-\\gamma H_y &=  \\imath \\omega \\epsilon_{xx} E_x - \\hat{\\partial}_y H_z \\\\
-\\gamma H_x &= -\\imath \\omega \\epsilon_{yy} E_y - \\hat{\\partial}_x H_z \\\\
-\\imath \\omega E_z &= \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x H_y - \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y H_x \\\\
-\\end{aligned}
-$$
-
-Now apply $\\gamma \\tilde{\\partial}_x$ to the last equation,
-then substitute in for $\\gamma H_x$ and $\\gamma H_y$:
 
 $$
-\\begin{aligned}
-\\gamma \\tilde{\\partial}_x \\imath \\omega E_z &= \\gamma \\tilde{\\partial}_x \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x H_y
-                                                - \\gamma \\tilde{\\partial}_x \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y H_x \\\\
-        &= \\tilde{\\partial}_x \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x ( \\imath \\omega \\epsilon_{xx} E_x - \\hat{\\partial}_y H_z)
-         - \\tilde{\\partial}_x \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y (-\\imath \\omega \\epsilon_{yy} E_y - \\hat{\\partial}_x H_z)  \\\\
-        &= \\tilde{\\partial}_x \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x ( \\imath \\omega \\epsilon_{xx} E_x)
-         - \\tilde{\\partial}_x \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y (-\\imath \\omega \\epsilon_{yy} E_y)  \\\\
-\\gamma \\tilde{\\partial}_x E_z &= \\tilde{\\partial}_x \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x (\\epsilon_{xx} E_x)
-                                  + \\tilde{\\partial}_x \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y (\\epsilon_{yy} E_y) \\\\
-\\end{aligned}
+\begin{aligned}
+\imath \beta H_y &=  \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z \\
+\imath \beta H_x &= -\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z \\
+\imath \omega E_z &= \frac{1}{\epsilon_{zz}} \hat{\partial}_x H_y - \frac{1}{\epsilon_{zz}} \hat{\partial}_y H_x \\
+\end{aligned}
 $$
 
-With a similar approach (but using $\\gamma \\tilde{\\partial}_y$ instead), we can get
+Now apply $\imath \beta \tilde{\partial}_x$ to the last equation,
+then substitute in for $\imath \beta H_x$ and $\imath \beta H_y$:
 
 $$
-\\begin{aligned}
-\\gamma \\tilde{\\partial}_y E_z &= \\tilde{\\partial}_y \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x (\\epsilon_{xx} E_x)
-                                  + \\tilde{\\partial}_y \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y (\\epsilon_{yy} E_y) \\\\
-\\end{aligned}
+\begin{aligned}
+\imath \beta \tilde{\partial}_x \imath \omega E_z &= \imath \beta \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x H_y
+                                                   - \imath \beta \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y H_x \\
+        &= \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x ( \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z)
+         - \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y (-\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z)  \\
+        &= \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x ( \imath \omega \epsilon_{xx} E_x)
+         - \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y (-\imath \omega \epsilon_{yy} E_y)  \\
+\imath \beta \tilde{\partial}_x E_z &= \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x)
+                                     + \tilde{\partial}_x \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) \\
+\end{aligned}
 $$
 
-We can combine this equation for $\\gamma \\tilde{\\partial}_y E_z$ with
-the unused $\\imath \\omega \\mu_{xx} H_x$ and $\\imath \\omega \\mu_{yy} H_y$ equations to get
+With a similar approach (but using $\imath \beta \tilde{\partial}_y$ instead), we can get
 
 $$
-\\begin{aligned}
--\\imath \\omega \\mu_{xx} \\gamma H_x &=  \\gamma^2 E_y + \\gamma \\tilde{\\partial}_y E_z \\\\
--\\imath \\omega \\mu_{xx} \\gamma H_x &=  \\gamma^2 E_y + \\tilde{\\partial}_y (
-                                      \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x (\\epsilon_{xx} E_x)
-                                    + \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y (\\epsilon_{yy} E_y)
-                                    )\\\\
-\\end{aligned}
+\begin{aligned}
+\imath \beta \tilde{\partial}_y E_z &= \tilde{\partial}_y \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x)
+                                     + \tilde{\partial}_y \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y) \\
+\end{aligned}
+$$
+
+We can combine this equation for $\imath \beta \tilde{\partial}_y E_z$ with
+the unused $\imath \omega \mu_{xx} H_x$ and $\imath \omega \mu_{yy} H_y$ equations to get
+
+$$
+\begin{aligned}
+-\imath \omega \mu_{xx} \imath \beta H_x &=  -\beta^2 E_y + \imath \beta \tilde{\partial}_y E_z \\
+-\imath \omega \mu_{xx} \imath \beta H_x &=  -\beta^2 E_y + \tilde{\partial}_y (
+                                      \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x)
+                                    + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y)
+                                    )\\
+\end{aligned}
 $$
 
 and
 
 $$
-\\begin{aligned}
--\\imath \\omega \\mu_{yy} \\gamma H_y &= -\\gamma^2 E_x - \\gamma \\tilde{\\partial}_x E_z \\\\
--\\imath \\omega \\mu_{yy} \\gamma H_y &= -\\gamma^2 E_x - \\tilde{\\partial}_x (
-                                      \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x (\\epsilon_{xx} E_x)
-                                    + \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y (\\epsilon_{yy} E_y)
-                                    )\\\\
-\\end{aligned}
+\begin{aligned}
+-\imath \omega \mu_{yy} \imath \beta H_y &= \beta^2 E_x - \imath \beta \tilde{\partial}_x E_z \\
+-\imath \omega \mu_{yy} \imath \beta H_y &= \beta^2 E_x - \tilde{\partial}_x (
+                                      \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x)
+                                    + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y)
+                                    )\\
+\end{aligned}
 $$
 
-However, based on our rewritten equation for $\\gamma H_x$ and the so-far unused
-equation for $\\imath \\omega \\mu_{zz} H_z$ we can also write
+However, based on our rewritten equation for $\imath \beta H_x$ and the so-far unused
+equation for $\imath \omega \mu_{zz} H_z$ we can also write
 
 $$
-\\begin{aligned}
--\\imath \\omega \\mu_{xx} (\\gamma H_x) &= -\\imath \\omega \\mu_{xx} (-\\imath \\omega \\epsilon_{yy} E_y - \\hat{\\partial}_x H_z) \\\\
-                                         &= -\\omega^2 \\mu_{xx} \\epsilon_{yy} E_y
-                                            +\\imath \\omega \\mu_{xx} \\hat{\\partial}_x (
-                                                \\frac{1}{-\\imath \\omega \\mu_{zz}} (\\tilde{\\partial}_x E_y - \\tilde{\\partial}_y E_x)) \\\\
-                                         &= -\\omega^2 \\mu_{xx} \\epsilon_{yy} E_y
-                                            -\\mu_{xx} \\hat{\\partial}_x \\frac{1}{\\mu_{zz}} (\\tilde{\\partial}_x E_y - \\tilde{\\partial}_y E_x) \\\\
-\\end{aligned}
+\begin{aligned}
+-\imath \omega \mu_{xx} (\imath \beta H_x) &= -\imath \omega \mu_{xx} (-\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z) \\
+                 &= -\omega^2 \mu_{xx} \epsilon_{yy} E_y + \imath \omega \mu_{xx} \hat{\partial}_x (
+                         \frac{1}{-\imath \omega \mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x)) \\
+                 &= -\omega^2 \mu_{xx} \epsilon_{yy} E_y
+                         -\mu_{xx} \hat{\partial}_x \frac{1}{\mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x) \\
+\end{aligned}
 $$
 
 and, similarly,
 
 $$
-\\begin{aligned}
--\\imath \\omega \\mu_{yy} (\\gamma H_y) &= \\omega^2 \\mu_{yy} \\epsilon_{xx} E_x
-                                           +\\mu_{yy} \\hat{\\partial}_y \\frac{1}{\\mu_{zz}} (\\tilde{\\partial}_x E_y - \\tilde{\\partial}_y E_x) \\\\
-\\end{aligned}
+\begin{aligned}
+-\imath \omega \mu_{yy} (\imath \beta H_y) &= \omega^2 \mu_{yy} \epsilon_{xx} E_x
+                                           +\mu_{yy} \hat{\partial}_y \frac{1}{\mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x) \\
+\end{aligned}
 $$
 
 By combining both pairs of expressions, we get
 
 $$
-\\begin{aligned}
--\\gamma^2 E_x - \\tilde{\\partial}_x (
-    \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x (\\epsilon_{xx} E_x)
-  + \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y (\\epsilon_{yy} E_y)
-    ) &= \\omega^2 \\mu_{yy} \\epsilon_{xx} E_x
-        +\\mu_{yy} \\hat{\\partial}_y \\frac{1}{\\mu_{zz}} (\\tilde{\\partial}_x E_y - \\tilde{\\partial}_y E_x) \\\\
-\\gamma^2 E_y + \\tilde{\\partial}_y (
-    \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_x (\\epsilon_{xx} E_x)
-  + \\frac{1}{\\epsilon_{zz}} \\hat{\\partial}_y (\\epsilon_{yy} E_y)
-    ) &= -\\omega^2 \\mu_{xx} \\epsilon_{yy} E_y
-         -\\mu_{xx} \\hat{\\partial}_x \\frac{1}{\\mu_{zz}} (\\tilde{\\partial}_x E_y - \\tilde{\\partial}_y E_x) \\\\
-\\end{aligned}
+\begin{aligned}
+\beta^2 E_x - \tilde{\partial}_x (
+    \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x)
+  + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y)
+    ) &= \omega^2 \mu_{yy} \epsilon_{xx} E_x
+        +\mu_{yy} \hat{\partial}_y \frac{1}{\mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x) \\
+-\beta^2 E_y + \tilde{\partial}_y (
+    \frac{1}{\epsilon_{zz}} \hat{\partial}_x (\epsilon_{xx} E_x)
+  + \frac{1}{\epsilon_{zz}} \hat{\partial}_y (\epsilon_{yy} E_y)
+    ) &= -\omega^2 \mu_{xx} \epsilon_{yy} E_y
+         -\mu_{xx} \hat{\partial}_x \frac{1}{\mu_{zz}} (\tilde{\partial}_x E_y - \tilde{\partial}_y E_x) \\
+\end{aligned}
 $$
 
 Using these, we can construct the eigenvalue problem
 
 $$
-\\beta^2 \\begin{bmatrix} E_x \\\\
-                          E_y \\end{bmatrix} =
-    (\\omega^2 \\begin{bmatrix} \\mu_{yy} \\epsilon_{xx} & 0 \\\\
-                                                       0 & \\mu_{xx} \\epsilon_{yy} \\end{bmatrix} +
-                 \\begin{bmatrix} -\\mu_{yy} \\hat{\\partial}_y \\\\
-                                   \\mu_{xx} \\hat{\\partial}_x \\end{bmatrix} \\mu_{zz}^{-1}
-                 \\begin{bmatrix} -\\tilde{\\partial}_y & \\tilde{\\partial}_x \\end{bmatrix} +
-      \\begin{bmatrix} \\tilde{\\partial}_x \\\\
-                       \\tilde{\\partial}_y \\end{bmatrix} \\epsilon_{zz}^{-1}
-                 \\begin{bmatrix} \\hat{\\partial}_x \\epsilon_{xx} & \\hat{\\partial}_y \\epsilon_{yy} \\end{bmatrix})
-    \\begin{bmatrix} E_x \\\\
-                     E_y \\end{bmatrix}
+\beta^2 \begin{bmatrix} E_x \\
+                        E_y \end{bmatrix} =
+    (\omega^2 \begin{bmatrix} \mu_{yy} \epsilon_{xx} & 0 \\
+                                                   0 & \mu_{xx} \epsilon_{yy} \end{bmatrix} +
+              \begin{bmatrix} -\mu_{yy} \hat{\partial}_y \\
+                               \mu_{xx} \hat{\partial}_x \end{bmatrix} \mu_{zz}^{-1}
+              \begin{bmatrix} -\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} +
+      \begin{bmatrix} \tilde{\partial}_x \\
+                      \tilde{\partial}_y \end{bmatrix} \epsilon_{zz}^{-1}
+                 \begin{bmatrix} \hat{\partial}_x \epsilon_{xx} & \hat{\partial}_y \epsilon_{yy} \end{bmatrix})
+    \begin{bmatrix} E_x \\
+                    E_y \end{bmatrix}
 $$
 
-where $\\gamma = \\imath\\beta$. In the literature, $\\beta$ is usually used to denote
-the lossless/real part of the propagation constant, but in `meanas` it is allowed to
-be complex.
+In the literature, $\beta$ is usually used to denote the lossless/real part of the propagation constant,
+but in `meanas` it is allowed to be complex.
 
 An equivalent eigenvalue problem can be formed using the $H_x$ and $H_y$ fields, if those are more convenient.
 
-Note that $E_z$ was never discretized, so $\\gamma$ and $\\beta$ will need adjustment
-to account for numerical dispersion if the result is introduced into a space with a discretized z-axis.
+Note that $E_z$ was never discretized, so $\beta$ will need adjustment to account for numerical dispersion
+if the result is introduced into a space with a discretized z-axis.
 
 
 """
 # TODO update module docs
 
 from typing import Any
+from collections.abc import Sequence
 import numpy
 from numpy.typing import NDArray, ArrayLike
 from numpy.linalg import norm
-import scipy.sparse as sparse       # type: ignore
+from scipy import sparse
 
 from ..fdmath.operators import deriv_forward, deriv_back, cross
-from ..fdmath import unvec, dx_lists_t, vfdfield_t, vcfdfield_t
+from ..fdmath import vec, unvec, dx_lists_t, vfdfield_t, vcfdfield_t
 from ..eigensolvers import signed_eigensolve, rayleigh_quotient_iteration
 
 
@@ -198,7 +198,7 @@ def operator_e(
         epsilon: vfdfield_t,
         mu: vfdfield_t | None = None,
         ) -> sparse.spmatrix:
-    """
+    r"""
     Waveguide operator of the form
 
         omega**2 * mu * epsilon +
@@ -210,18 +210,18 @@ def operator_e(
     More precisely, the operator is
 
     $$
-    \\omega^2 \\begin{bmatrix} \\mu_{yy} \\epsilon_{xx} & 0 \\\\
-                                                       0 & \\mu_{xx} \\epsilon_{yy} \\end{bmatrix} +
-                 \\begin{bmatrix} -\\mu_{yy} \\hat{\\partial}_y \\\\
-                                   \\mu_{xx} \\hat{\\partial}_x \\end{bmatrix} \\mu_{zz}^{-1}
-                 \\begin{bmatrix} -\\tilde{\\partial}_y & \\tilde{\\partial}_x \\end{bmatrix} +
-      \\begin{bmatrix} \\tilde{\\partial}_x \\\\
-                       \\tilde{\\partial}_y \\end{bmatrix} \\epsilon_{zz}^{-1}
-                 \\begin{bmatrix} \\hat{\\partial}_x \\epsilon_{xx} & \\hat{\\partial}_y \\epsilon_{yy} \\end{bmatrix}
+    \omega^2 \begin{bmatrix} \mu_{yy} \epsilon_{xx} & 0 \\
+                                                      0 & \mu_{xx} \epsilon_{yy} \end{bmatrix} +
+                 \begin{bmatrix} -\mu_{yy} \hat{\partial}_y \\
+                                   \mu_{xx} \hat{\partial}_x \end{bmatrix} \mu_{zz}^{-1}
+                 \begin{bmatrix} -\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} +
+      \begin{bmatrix} \tilde{\partial}_x \\
+                       \tilde{\partial}_y \end{bmatrix} \epsilon_{zz}^{-1}
+                 \begin{bmatrix} \hat{\partial}_x \epsilon_{xx} & \hat{\partial}_y \epsilon_{yy} \end{bmatrix}
     $$
 
-    $\\tilde{\\partial}_x$ and $\\hat{\\partial}_x$ are the forward and backward derivatives along x,
-    and each $\\epsilon_{xx}$, $\\mu_{yy}$, etc. is a diagonal matrix containing the vectorized material
+    $\tilde{\partial}_x$ and $\hat{\partial}_x$ are the forward and backward derivatives along x,
+    and each $\epsilon_{xx}$, $\mu_{yy}$, etc. is a diagonal matrix containing the vectorized material
     property distribution.
 
     This operator can be used to form an eigenvalue problem of the form
@@ -253,9 +253,11 @@ def operator_e(
     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
+    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
 
 
@@ -265,7 +267,7 @@ def operator_h(
         epsilon: vfdfield_t,
         mu: vfdfield_t | None = None,
         ) -> sparse.spmatrix:
-    """
+    r"""
     Waveguide operator of the form
 
         omega**2 * epsilon * mu +
@@ -277,18 +279,18 @@ def operator_h(
     More precisely, the operator is
 
     $$
-    \\omega^2 \\begin{bmatrix} \\epsilon_{yy} \\mu_{xx} & 0 \\\\
-                                                      0 & \\epsilon_{xx} \\mu_{yy} \\end{bmatrix} +
-                 \\begin{bmatrix} -\\epsilon_{yy} \\tilde{\\partial}_y \\\\
-                                   \\epsilon_{xx} \\tilde{\\partial}_x \\end{bmatrix} \\epsilon_{zz}^{-1}
-                 \\begin{bmatrix} -\\hat{\\partial}_y & \\hat{\\partial}_x \\end{bmatrix} +
-      \\begin{bmatrix} \\hat{\\partial}_x \\\\
-                       \\hat{\\partial}_y \\end{bmatrix} \\mu_{zz}^{-1}
-                 \\begin{bmatrix} \\tilde{\\partial}_x \\mu_{xx} & \\tilde{\\partial}_y \\mu_{yy} \\end{bmatrix}
+    \omega^2 \begin{bmatrix} \epsilon_{yy} \mu_{xx} & 0 \\
+                                                  0 & \epsilon_{xx} \mu_{yy} \end{bmatrix} +
+                 \begin{bmatrix} -\epsilon_{yy} \tilde{\partial}_y \\
+                                  \epsilon_{xx} \tilde{\partial}_x \end{bmatrix} \epsilon_{zz}^{-1}
+                 \begin{bmatrix} -\hat{\partial}_y & \hat{\partial}_x \end{bmatrix} +
+      \begin{bmatrix} \hat{\partial}_x \\
+                      \hat{\partial}_y \end{bmatrix} \mu_{zz}^{-1}
+                 \begin{bmatrix} \tilde{\partial}_x \mu_{xx} & \tilde{\partial}_y \mu_{yy} \end{bmatrix}
     $$
 
-    $\\tilde{\\partial}_x$ and $\\hat{\\partial}_x$ are the forward and backward derivatives along x,
-    and each $\\epsilon_{xx}$, $\\mu_{yy}$, etc. is a diagonal matrix containing the vectorized material
+    $\tilde{\partial}_x$ and $\hat{\partial}_x$ are the forward and backward derivatives along x,
+    and each $\epsilon_{xx}$, $\mu_{yy}$, etc. is a diagonal matrix containing the vectorized material
     property distribution.
 
     This operator can be used to form an eigenvalue problem of the form
@@ -320,10 +322,11 @@ def operator_h(
     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
-
+    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
 
 
@@ -410,18 +413,13 @@ def _normalized_fields(
     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)
+    Sz_tavg = inner_product(e, h, dxes=dxes, prop_phase=prop_phase, conj_h=True).real
+    assert Sz_tavg > 0, f'Found a mode propagating in the wrong direction! {Sz_tavg=}'
 
-    energy = epsilon * e.conj() * e
+    energy = numpy.real(epsilon * e.conj() * e)
 
     norm_amplitude = 1 / numpy.sqrt(Sz_tavg)
     norm_angle = -numpy.angle(e[energy.argmax()])       # Will randomly add a negative sign when mode is symmetric
@@ -431,6 +429,7 @@ def _normalized_fields(
     sign = numpy.sign(E_weighted[:,
                                  :max(shape[0] // 2, 1),
                                  :max(shape[1] // 2, 1)].real.sum())
+    assert sign != 0
 
     norm_factor = sign * norm_amplitude * numpy.exp(1j * norm_angle)
 
@@ -531,10 +530,37 @@ def exy2e(
         dxes: dx_lists_t,
         epsilon: vfdfield_t,
         ) -> sparse.spmatrix:
-    """
+    r"""
     Operator which transforms the vector `e_xy` containing the vectorized E_x and E_y fields,
      into a vectorized E containing all three E components
 
+    From the operator derivation (see module docs), we have
+
+    $$
+    \imath \omega \epsilon_{zz} E_z = \hat{\partial}_x H_y - \hat{\partial}_y H_x \\
+    $$
+
+    as well as the intermediate equations
+
+    $$
+    \begin{aligned}
+    \imath \beta H_y &=  \imath \omega \epsilon_{xx} E_x - \hat{\partial}_y H_z \\
+    \imath \beta H_x &= -\imath \omega \epsilon_{yy} E_y - \hat{\partial}_x H_z \\
+    \end{aligned}
+    $$
+
+    Combining these, we get
+
+    $$
+    \begin{aligned}
+    E_z &= \frac{1}{- \omega \beta \epsilon_{zz}} ((
+             \hat{\partial}_y \hat{\partial}_x H_z
+            -\hat{\partial}_x \hat{\partial}_y H_z)
+          + \imath \omega (\hat{\partial}_x \epsilon_{xx} E_x + \hat{\partial}_y \epsilon{yy} E_y))
+        &= \frac{1}{\imath \beta \epsilon_{zz}} (\hat{\partial}_x \epsilon_{xx} E_x + \hat{\partial}_y \epsilon{yy} E_y)
+    \end{aligned}
+    $$
+
     Args:
         wavenumber: Wavenumber assuming fields have z-dependence of `exp(-i * wavenumber * z)`
                     It should satisfy `operator_e() @ e_xy == wavenumber**2 * e_xy`
@@ -717,8 +743,111 @@ def e_err(
     return float(norm(op) / norm(e))
 
 
+def sensitivity(
+        e_norm: vcfdfield_t,
+        h_norm: vcfdfield_t,
+        wavenumber: complex,
+        omega: complex,
+        dxes: dx_lists_t,
+        epsilon: vfdfield_t,
+        mu: vfdfield_t | None = None,
+        ) -> vcfdfield_t:
+    r"""
+      Given a waveguide structure (`dxes`, `epsilon`, `mu`) and mode fields
+    (`e_norm`, `h_norm`, `wavenumber`, `omega`), calculates the sensitivity of the wavenumber
+    $\beta$ to changes in the dielectric structure $\epsilon$.
+
+    The output is a vector of the same size as `vec(epsilon)`, with each element specifying the
+    sensitivity of `wavenumber` to changes in the corresponding element in `vec(epsilon)`, i.e.
+
+    $$sens_{i} = \frac{\partial\beta}{\partial\epsilon_i}$$
+
+    An adjoint approach is used to calculate the sensitivity; the derivation is provided here:
+
+    Starting with the eigenvalue equation
+
+    $$\beta^2 E_{xy} = A_E E_{xy}$$
+
+    where $A_E$ is the waveguide operator from `operator_e()`, and $E_{xy} = \begin{bmatrix} E_x \\
+                                                                                             E_y \end{bmatrix}$,
+    we can differentiate with respect to one of the $\epsilon$ elements (i.e. at one Yee grid point), $\epsilon_i$:
+
+    $$
+    (2 \beta) \partial_{\epsilon_i}(\beta) E_{xy} + \beta^2 \partial_{\epsilon_i} E_{xy}
+        = \partial_{\epsilon_i}(A_E) E_{xy} + A_E \partial_{\epsilon_i} E_{xy}
+    $$
+
+    We then multiply by $H_{yx}^\star = \begin{bmatrix}H_y^\star \\ -H_x^\star \end{bmatrix}$ from the left:
+
+    $$
+    (2 \beta) \partial_{\epsilon_i}(\beta) H_{yx}^\star E_{xy} + \beta^2 H_{yx}^\star \partial_{\epsilon_i} E_{xy}
+        = H_{yx}^\star \partial_{\epsilon_i}(A_E) E_{xy} + H_{yx}^\star A_E \partial_{\epsilon_i} E_{xy}
+    $$
+
+    However, $H_{yx}^\star$ is actually a left-eigenvector of $A_E$. This can be verified by inspecting
+    the form of `operator_h` ($A_H$) and comparing its conjugate transpose to `operator_e` ($A_E$). Also, note
+    $H_{yx}^\star \cdot E_{xy} = H^\star \times E$ recalls the mode orthogonality relation. See doi:10.5194/ars-9-85-201
+    for a similar approach. Therefore,
+
+    $$
+    H_{yx}^\star A_E \partial_{\epsilon_i} E_{xy} = \beta^2 H_{yx}^\star \partial_{\epsilon_i} E_{xy}
+    $$
+
+    and we can simplify to
+
+    $$
+    \partial_{\epsilon_i}(\beta)
+        = \frac{1}{2 \beta} \frac{H_{yx}^\star \partial_{\epsilon_i}(A_E) E_{xy} }{H_{yx}^\star E_{xy}}
+    $$
+
+    This expression can be quickly calculated for all $i$ by writing out the various terms of
+    $\partial_{\epsilon_i} A_E$ and recognizing that the vector-matrix-vector products (i.e. scalars)
+    $sens_i = \vec{v}_{left} \partial_{\epsilon_i} (\epsilon_{xyz}) \vec{v}_{right}$, indexed by $i$, can be expressed as
+    elementwise multiplications $\vec{sens} = \vec{v}_{left} \star \vec{v}_{right}$
+
+
+    Args:
+        e_norm: Normalized, vectorized E_xyz field for the mode. E.g. as returned by `normalized_fields_e`.
+        h_norm: Normalized, vectorized H_xyz field for the mode. E.g. as returned by `normalized_fields_e`.
+        wavenumber: Propagation constant for the mode. The z-axis is assumed to be continuous (i.e. without numerical dispersion).
+        omega: The angular frequency of the system.
+        dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
+        epsilon: Vectorized dielectric constant grid
+        mu: Vectorized magnetic permeability grid (default 1 everywhere)
+
+    Returns:
+        Sparse matrix representation of the operator.
+    """
+    if mu is None:
+        mu = numpy.ones_like(epsilon)
+
+    Dfx, Dfy = deriv_forward(dxes[0])
+    Dbx, Dby = deriv_back(dxes[1])
+
+    eps_x, eps_y, eps_z = numpy.split(epsilon, 3)
+    eps_xy = sparse.diags(numpy.hstack((eps_x, eps_y)))
+    eps_z_inv = sparse.diags(1 / eps_z)
+
+    mu_x, mu_y, _mu_z = numpy.split(mu, 3)
+    mu_yx = sparse.diags(numpy.hstack((mu_y, mu_x)))
+
+    da_exxhyy = vec(dxes[1][0][:, None] * dxes[0][1][None, :])
+    da_eyyhxx = vec(dxes[1][1][None, :] * dxes[0][0][:, None])
+    ev_xy = numpy.concatenate(numpy.split(e_norm, 3)[:2]) * numpy.concatenate([da_exxhyy, da_eyyhxx])
+    hx, hy, hz = numpy.split(h_norm, 3)
+    hv_yx_conj = numpy.conj(numpy.concatenate([hy, -hx]))
+
+    sens_xy1 = (hv_yx_conj @ (omega * omega * mu_yx)) * ev_xy
+    sens_xy2 = (hv_yx_conj @ sparse.vstack((Dfx, Dfy)) @ eps_z_inv @ sparse.hstack((Dbx, Dby))) * ev_xy
+    sens_z   = (hv_yx_conj @ sparse.vstack((Dfx, Dfy)) @ (-eps_z_inv * eps_z_inv)) * (sparse.hstack((Dbx, Dby)) @ eps_xy @ ev_xy)
+    norm = hv_yx_conj @ ev_xy
+
+    sens_tot = numpy.concatenate([sens_xy1 + sens_xy2, sens_z]) / (2 * wavenumber * norm)
+    return sens_tot
+
+
 def solve_modes(
-        mode_numbers: list[int],
+        mode_numbers: Sequence[int],
         omega: complex,
         dxes: dx_lists_t,
         epsilon: vfdfield_t,
@@ -739,32 +868,38 @@ def solve_modes(
             ability to find the correct mode. Default 2.
 
     Returns:
-        e_xys: list of vfdfield_t specifying fields
+        e_xys: NDArray of vfdfield_t specifying fields. First dimension is mode number.
         wavenumbers: list of wavenumbers
     """
 
-    '''
-    Solve for the largest-magnitude eigenvalue of the real operator
-    '''
+    #
+    # Solve for the largest-magnitude eigenvalue of the real operator
+    #
     dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes]
     mu_real = None if mu is None else numpy.real(mu)
     A_r = operator_e(numpy.real(omega), dxes_real, numpy.real(epsilon), mu_real)
 
     eigvals, eigvecs = signed_eigensolve(A_r, max(mode_numbers) + mode_margin)
-    e_xys = eigvecs[:, -(numpy.array(mode_numbers) + 1)]
+    keep_inds = -(numpy.array(mode_numbers) + 1)
+    e_xys = eigvecs[:, keep_inds].T
+    eigvals = eigvals[keep_inds]
 
-    '''
-    Now solve for the eigenvector of the full operator, using the real operator's
-     eigenvector as an initial guess for Rayleigh quotient iteration.
-    '''
+    #
+    # Now solve for the eigenvector of the full operator, using the real operator's
+    #  eigenvector as an initial guess for Rayleigh quotient iteration.
+    #
     A = operator_e(omega, dxes, epsilon, mu)
     for nn in range(len(mode_numbers)):
-        eigvals[nn], e_xys[:, nn] = rayleigh_quotient_iteration(A, e_xys[:, nn])
+        eigvals[nn], e_xys[nn, :] = rayleigh_quotient_iteration(A, e_xys[nn, :])
 
     # Calculate the wave-vector (force the real part to be positive)
     wavenumbers = numpy.sqrt(eigvals)
     wavenumbers *= numpy.sign(numpy.real(wavenumbers))
 
+    order = wavenumbers.argsort()[::-1]
+    e_xys = e_xys[order]
+    wavenumbers = wavenumbers[order]
+
     return e_xys, wavenumbers
 
 
@@ -786,4 +921,38 @@ def solve_mode(
     """
     kwargs['mode_numbers'] = [mode_number]
     e_xys, wavenumbers = solve_modes(*args, **kwargs)
-    return e_xys[:, 0], wavenumbers[0]
+    return e_xys[0], wavenumbers[0]
+
+
+def inner_product(    # TODO documentation
+        e1: vcfdfield_t,
+        h2: vcfdfield_t,
+        dxes: dx_lists_t,
+        prop_phase: float = 0,
+        conj_h: bool = False,
+        trapezoid: bool = False,
+        ) -> complex:
+
+    shape = [s.size for s in dxes[0]]
+
+    # 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)
+
+    E1 = unvec(e1, shape)
+    H2 = unvec(h2, shape) * phase
+
+    if conj_h:
+        H2 = numpy.conj(H2)
+
+    # Find time-averaged Sz and normalize to it
+    dxes_real = [[numpy.real(dxyz) for dxyz in dxeh] for dxeh in dxes]
+    if trapezoid:
+        Sz_a = numpy.trapezoid(numpy.trapezoid(E1[0] * H2[1], numpy.cumsum(dxes_real[0][1])), numpy.cumsum(dxes_real[1][0]))
+        Sz_b = numpy.trapezoid(numpy.trapezoid(E1[1] * H2[0], numpy.cumsum(dxes_real[0][0])), numpy.cumsum(dxes_real[1][1]))
+    else:
+        Sz_a = E1[0] * H2[1] * dxes_real[1][0][:, None] * dxes_real[0][1][None, :]
+        Sz_b = E1[1] * H2[0] * dxes_real[0][0][:, None] * dxes_real[1][1][None, :]
+    Sz = 0.5 * (Sz_a.sum() - Sz_b.sum())
+    return Sz
+
+

meanas/fdfd/waveguide_3d.py

@@ -4,9 +4,11 @@ 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 Sequence, Any
+from typing import Any
+from collections.abc import Sequence
 import numpy
 from numpy.typing import NDArray
+from numpy import complexfloating
 
 from ..fdmath import vec, unvec, dx_lists_t, fdfield_t, cfdfield_t
 from . import operators, waveguide_2d
@@ -21,7 +23,7 @@ def solve_mode(
         slices: Sequence[slice],
         epsilon: fdfield_t,
         mu: fdfield_t | None = None,
-        ) -> dict[str, complex | NDArray[numpy.float_]]:
+        ) -> dict[str, complex | NDArray[complexfloating]]:
     """
     Given a 3D grid, selects a slice from the grid and attempts to
      solve for an eigenmode propagating through that slice.
@@ -40,8 +42,8 @@ def solve_mode(
     Returns:
         ```
         {
-            'E': list[NDArray[numpy.float_]],
-            'H': list[NDArray[numpy.float_]],
+            'E': NDArray[complexfloating],
+            'H': NDArray[complexfloating],
             'wavenumber': complex,
         }
         ```
@@ -51,9 +53,9 @@ def solve_mode(
 
     slices = tuple(slices)
 
-    '''
-    Solve the 2D problem in the specified plane
-    '''
+    #
+    # 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)
@@ -71,9 +73,10 @@ def solve_mode(
     }
     e_xy, wavenumber_2d = waveguide_2d.solve_mode(mode_number, **args_2d)
 
-    '''
-    Apply corrections and expand to 3D
-    '''
+    #
+    # Apply corrections and expand to 3D
+    #
+
     # Correct wavenumber to account for numerical dispersion.
     wavenumber = 2 / dx_prop * numpy.arcsin(wavenumber_2d * dx_prop / 2)
 

meanas/fdfd/waveguide_cyl.py

@@ -1,31 +1,102 @@
-"""
+r"""
 Operators and helper functions for cylindrical waveguides with unchanging cross-section.
 
-WORK IN PROGRESS, CURRENTLY BROKEN
+Waveguide operator is derived according to 10.1364/OL.33.001848.
+The curl equations in the complex coordinate system become
 
-As the z-dependence is known, all the functions in this file assume a 2D grid
+$$
+\begin{aligned}
+-\imath \omega \mu_{xx} H_x &= \tilde{\partial}_y E_z + \imath \beta frac{E_y}{\tilde{t}_x} \\
+-\imath \omega \mu_{yy} H_y &= -\imath \beta E_x - \frac{1}{\hat{t}_x} \tilde{\partial}_x \tilde{t}_x E_z \\
+-\imath \omega \mu_{zz} H_z &= \tilde{\partial}_x E_y - \tilde{\partial}_y E_x \\
+\imath \omega \epsilon_{xx} E_x &= \hat{\partial}_y H_z + \imath \beta \frac{H_y}{\hat{T}} \\
+\imath \omega \epsilon_{yy} E_y &= -\imath \beta H_x - \{1}{\tilde{t}_x} \hat{\partial}_x \hat{t}_x} H_z \\
+\imath \omega \epsilon_{zz} E_z &= \hat{\partial}_x H_y - \hat{\partial}_y H_x \\
+\end{aligned}
+$$
+
+where $t_x = 1 + \frac{\Delta_{x, m}}{R_0}$ is the grid spacing adjusted by the nominal radius $R0$.
+
+Rewrite the last three equations as
+
+$$
+\begin{aligned}
+\imath \beta H_y &=  \imath \omega \hat{t}_x \epsilon_{xx} E_x - \hat{t}_x \hat{\partial}_y H_z \\
+\imath \beta H_x &= -\imath \omega \hat{t}_x \epsilon_{yy} E_y - \hat{t}_x \hat{\partial}_x H_z \\
+\imath \omega E_z &= \frac{1}{\epsilon_{zz}} \hat{\partial}_x H_y - \frac{1}{\epsilon_{zz}} \hat{\partial}_y H_x \\
+\end{aligned}
+$$
+
+The derivation then follows the same steps as the straight waveguide, leading to the eigenvalue problem
+
+$$
+\beta^2 \begin{bmatrix} E_x \\
+                        E_y \end{bmatrix} =
+    (\omega^2 \begin{bmatrix} T_b T_b \mu_{yy} \epsilon_{xx} & 0 \\
+                                                            0 & T_a T_a \mu_{xx} \epsilon_{yy} \end{bmatrix} +
+              \begin{bmatrix} -T_b \mu_{yy} \hat{\partial}_y \\
+                               T_a \mu_{xx} \hat{\partial}_x \end{bmatrix} T_b \mu_{zz}^{-1}
+              \begin{bmatrix} -\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} +
+      \begin{bmatrix} \tilde{\partial}_x \\
+                      \tilde{\partial}_y \end{bmatrix} T_a \epsilon_{zz}^{-1}
+                 \begin{bmatrix} \hat{\partial}_x T_b \epsilon_{xx} & \hat{\partial}_y T_a \epsilon_{yy} \end{bmatrix})
+    \begin{bmatrix} E_x \\
+                    E_y \end{bmatrix}
+$$
+
+which resembles the straight waveguide eigenproblem with additonal $T_a$ and $T_b$ terms. These
+are diagonal matrices containing the $t_x$ values:
+
+$$
+\begin{aligned}
+T_a &=  1 + \frac{\Delta_{x, m               }}{R_0}
+T_b &=  1 + \frac{\Delta_{x, m + \frac{1}{2} }}{R_0}
+\end{aligned}
+
+
+TODO: consider 10.1364/OE.20.021583 for an alternate approach
+$$
+
+As in the straight waveguide case, 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 Any, cast
+from collections.abc import Sequence
+import logging
 
 import numpy
-import scipy.sparse as sparse       # type: ignore
+from numpy.typing import NDArray, ArrayLike
+from scipy import sparse
 
-from ..fdmath import vec, unvec, dx_lists_t, fdfield_t, vfdfield_t, cfdfield_t
+from ..fdmath import vec, unvec, dx_lists_t, vfdfield_t, vcfdfield_t
 from ..fdmath.operators import deriv_forward, deriv_back
 from ..eigensolvers import signed_eigensolve, rayleigh_quotient_iteration
+from . import waveguide_2d
+
+logger = logging.getLogger(__name__)
 
 
 def cylindrical_operator(
-        omega: complex,
+        omega: float,
         dxes: dx_lists_t,
         epsilon: vfdfield_t,
-        r0: float,
+        rmin: float,
         ) -> sparse.spmatrix:
-    """
+    r"""
     Cylindrical coordinate waveguide operator of the form
 
-    TODO
+    $$
+        (\omega^2 \begin{bmatrix} T_b T_b \mu_{yy} \epsilon_{xx} & 0 \\
+                                                                0 & T_a T_a \mu_{xx} \epsilon_{yy} \end{bmatrix} +
+                  \begin{bmatrix} -T_b \mu_{yy} \hat{\partial}_y \\
+                                   T_a \mu_{xx} \hat{\partial}_x \end{bmatrix} T_b \mu_{zz}^{-1}
+                  \begin{bmatrix} -\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} +
+          \begin{bmatrix} \tilde{\partial}_x \\
+                          \tilde{\partial}_y \end{bmatrix} T_a \epsilon_{zz}^{-1}
+                     \begin{bmatrix} \hat{\partial}_x T_b \epsilon_{xx} & \hat{\partial}_y T_a \epsilon_{yy} \end{bmatrix})
+        \begin{bmatrix} E_x \\
+                        E_y \end{bmatrix}
+    $$
 
     for use with a field vector of the form `[E_r, E_y]`.
 
@@ -35,12 +106,13 @@ def cylindrical_operator(
     which can then be solved for the eigenmodes of the system
     (an `exp(-i * wavenumber * theta)` theta-dependence is assumed for the fields).
 
+    (NOTE: See module docs and 10.1364/OL.33.001848)
+
     Args:
         omega: The angular frequency of the system
         dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
         epsilon: Vectorized dielectric constant grid
-        r0: Radius of curvature for the simulation. This should be the minimum value of
-            r within the simulation domain.
+        rmin: Radius at the left edge of the simulation domain (at minimum 'x')
 
     Returns:
         Sparse matrix representation of the operator
@@ -49,44 +121,34 @@ def cylindrical_operator(
     Dfx, Dfy = deriv_forward(dxes[0])
     Dbx, Dby = 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)))
+    Ta, Tb = dxes2T(dxes=dxes, rmin=rmin)
 
     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
+    eps_x = sparse.diags_array(eps_parts[0])
+    eps_y = sparse.diags_array(eps_parts[1])
+    eps_z_inv = sparse.diags_array(1 / eps_parts[2])
 
     omega2 = omega * omega
+    diag = sparse.block_diag
 
-    op = (omega2 * diag((Tx, Ty)) + pa) @ diag((a0, a1)) + \
-        - (sparse.bmat(((None, Ty), (Tx, None))) + pb / omega2) @ diag((b0, b1))
+    sq0 = omega2 * diag((Tb @ Tb @ eps_x,
+                         Ta @ Ta @ eps_y))
+    lin0 = sparse.vstack((-Tb @ Dby, Ta @ Dbx)) @ Tb @ sparse.hstack((-Dfy, Dfx))
+    lin1 = sparse.vstack((Dfx, Dfy)) @ Ta @ eps_z_inv @ sparse.hstack((Dbx @ Tb @ eps_x,
+                                                                       Dby @ Ta @ eps_y))
+    op = sq0 + lin0 + lin1
     return op
 
 
-def solve_mode(
-        mode_number: int,
-        omega: complex,
+def solve_modes(
+        mode_numbers: Sequence[int],
+        omega: float,
         dxes: dx_lists_t,
         epsilon: vfdfield_t,
-        r0: float,
-        ) -> dict[str, complex | cfdfield_t]:
+        rmin: float,
+        mode_margin: int = 2,
+        ) -> tuple[vcfdfield_t, NDArray[numpy.complex128]]:
     """
-    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.
 
@@ -96,48 +158,345 @@ def solve_mode(
         dxes: Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types.
               The first coordinate is assumed to be r, the second is y.
         epsilon: Dielectric constant
-        r0: Radius of curvature for the simulation. This should be the minimum value of
+        rmin: Radius of curvature for the simulation. This should be the minimum value of
                r within the simulation domain.
 
     Returns:
-        ```
-        {
-            'E': list[NDArray[numpy.complex_]],
-            'H': list[NDArray[numpy.complex_]],
-            'wavenumber': complex,
-        }
-        ```
+        e_xys: NDArray of vfdfield_t specifying fields. First dimension is mode number.
+        angular_wavenumbers: list of wavenumbers in 1/rad units.
     """
 
-    '''
-    Solve for the largest-magnitude eigenvalue of the real operator
-    '''
+    #
+    # 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 = 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)]
+    A_r = cylindrical_operator(numpy.real(omega), dxes_real, numpy.real(epsilon), rmin=rmin)
+    eigvals, eigvecs = signed_eigensolve(A_r, max(mode_numbers) + mode_margin)
+    keep_inds = -(numpy.array(mode_numbers) + 1)
+    e_xys = eigvecs[:, keep_inds].T
+    eigvals = eigvals[keep_inds]
 
-    '''
-    Now solve for the eigenvector of the full operator, using the real operator's
-     eigenvector as an initial guess for Rayleigh quotient iteration.
-    '''
-    A = cylindrical_operator(omega, dxes, epsilon, r0)
-    eigval, e_xy = rayleigh_quotient_iteration(A, e_xy)
+    #
+    # Now solve for the eigenvector of the full operator, using the real operator's
+    #  eigenvector as an initial guess for Rayleigh quotient iteration.
+    #
+    A = cylindrical_operator(omega, dxes, epsilon, rmin=rmin)
+    for nn in range(len(mode_numbers)):
+        eigvals[nn], e_xys[nn, :] = rayleigh_quotient_iteration(A, e_xys[nn, :])
 
     # Calculate the wave-vector (force the real part to be positive)
-    wavenumber = numpy.sqrt(eigval)
-    wavenumber *= numpy.sign(numpy.real(wavenumber))
+    wavenumbers = numpy.sqrt(eigvals)
+    wavenumbers *= numpy.sign(numpy.real(wavenumbers))
 
-    # TODO: Perform correction on wavenumber to account for numerical dispersion.
+    # Wavenumbers assume the mode is at rmin, which is unlikely
+    # Instead, return the wavenumber in inverse radians
+    angular_wavenumbers = wavenumbers * cast(complex, rmin)
 
-    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),
-    }
+    order = angular_wavenumbers.argsort()[::-1]
+    e_xys = e_xys[order]
+    angular_wavenumbers = angular_wavenumbers[order]
 
-    return fields
+    return e_xys, angular_wavenumbers
+
+
+def solve_mode(
+        mode_number: int,
+        *args: Any,
+        **kwargs: Any,
+        ) -> tuple[vcfdfield_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, angular_wavenumber)
+    """
+    kwargs['mode_numbers'] = [mode_number]
+    e_xys, angular_wavenumbers = solve_modes(*args, **kwargs)
+    return e_xys[0], angular_wavenumbers[0]
+
+
+def linear_wavenumbers(
+        e_xys: vcfdfield_t,
+        angular_wavenumbers: ArrayLike,
+        epsilon: vfdfield_t,
+        dxes: dx_lists_t,
+        rmin: float,
+        ) -> NDArray[numpy.complex128]:
+    """
+    Calculate linear wavenumbers (1/distance) based on angular wavenumbers (1/rad)
+      and the mode's energy distribution.
+
+    Args:
+        e_xys: Vectorized mode fields with shape (num_modes, 2 * x *y)
+        angular_wavenumbers: Wavenumbers assuming fields have theta-dependence of
+            `exp(-i * angular_wavenumber * theta)`. They should satisfy
+            `operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy`
+        epsilon: Vectorized dielectric constant grid with shape (3, x, y)
+        dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
+        rmin: Radius at the left edge of the simulation domain (at minimum 'x')
+
+    Returns:
+        NDArray containing the calculated linear (1/distance) wavenumbers
+    """
+    angular_wavenumbers = numpy.asarray(angular_wavenumbers)
+    mode_radii = numpy.empty_like(angular_wavenumbers, dtype=float)
+
+    wavenumbers = numpy.empty_like(angular_wavenumbers)
+    shape2d = (len(dxes[0][0]), len(dxes[0][1]))
+    epsilon2d = unvec(epsilon, shape2d)[:2]
+    grid_radii = rmin + numpy.cumsum(dxes[0][0])
+    for ii in range(angular_wavenumbers.size):
+        efield = unvec(e_xys[ii], shape2d, 2)
+        energy = numpy.real((efield * efield.conj()) * epsilon2d)
+        energy_vs_x = energy.sum(axis=(0, 2))
+        mode_radii[ii] = (grid_radii * energy_vs_x).sum() / energy_vs_x.sum()
+
+    logger.info(f'{mode_radii=}')
+    lin_wavenumbers = angular_wavenumbers / mode_radii
+    return lin_wavenumbers
+
+
+def exy2h(
+        angular_wavenumber: complex,
+        omega: float,
+        dxes: dx_lists_t,
+        rmin: float,
+        epsilon: vfdfield_t,
+        mu: vfdfield_t | None = 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:
+        angular_wavenumber: Wavenumber assuming fields have theta-dependence of
+            `exp(-i * angular_wavenumber * theta)`. It should satisfy
+            `operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy`
+        omega: The angular frequency of the system
+        dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
+        rmin: Radius at the left edge of the simulation domain (at minimum 'x')
+        epsilon: Vectorized dielectric constant grid
+        mu: Vectorized magnetic permeability grid (default 1 everywhere)
+
+    Returns:
+        Sparse matrix representing the operator.
+    """
+    e2hop = e2h(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, mu=mu)
+    return e2hop @ exy2e(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon)
+
+
+def exy2e(
+        angular_wavenumber: complex,
+        omega: float,
+        dxes: dx_lists_t,
+        rmin: float,
+        epsilon: vfdfield_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
+
+    Unlike the straight waveguide case, the H_z components do not cancel and must be calculated
+    from E_x and E_y in order to then calculate E_z.
+
+    Args:
+        angular_wavenumber: Wavenumber assuming fields have theta-dependence of
+            `exp(-i * angular_wavenumber * theta)`. It should satisfy
+            `operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy`
+        omega: The angular frequency of the system
+        dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
+        rmin: Radius at the left edge of the simulation domain (at minimum 'x')
+        epsilon: Vectorized dielectric constant grid
+
+    Returns:
+        Sparse matrix representing the operator.
+    """
+    Dfx, Dfy = deriv_forward(dxes[0])
+    Dbx, Dby = deriv_back(dxes[1])
+    wavenumber = angular_wavenumber / rmin
+
+    Ta, Tb = dxes2T(dxes=dxes, rmin=rmin)
+    Tai = sparse.diags_array(1 / Ta.diagonal())
+    Tbi = sparse.diags_array(1 / Tb.diagonal())
+
+    epsilon_parts = numpy.split(epsilon, 3)
+    epsilon_x, epsilon_y = (sparse.diags_array(epsi) for epsi in epsilon_parts[:2])
+    epsilon_z_inv = sparse.diags_array(1 / epsilon_parts[2])
+
+    n_pts = dxes[0][0].size * dxes[0][1].size
+    zeros = sparse.coo_array((n_pts, n_pts))
+    keep_x = sparse.block_array([[sparse.eye_array(n_pts), None], [None, zeros]])
+    keep_y = sparse.block_array([[zeros, None], [None, sparse.eye_array(n_pts)]])
+
+    mu_z = numpy.ones(n_pts)
+    mu_z_inv = sparse.diags_array(1 / mu_z)
+    exy2hz = 1 / (-1j * omega) * mu_z_inv @ sparse.hstack((Dfy, -Dfx))
+    hxy2ez = 1 / (1j * omega) * epsilon_z_inv @ sparse.hstack((Dby, -Dbx))
+
+    exy2hy = Tb / (1j * wavenumber) @ (-1j * omega * sparse.hstack((epsilon_x, zeros)) - Dby @ exy2hz)
+    exy2hx = Tb / (1j * wavenumber) @ ( 1j * omega * sparse.hstack((zeros, epsilon_y)) - Tai @ Dbx @ Tb @ exy2hz)
+
+    exy2ez = hxy2ez @ sparse.vstack((exy2hx, exy2hy))
+
+    op = sparse.vstack((sparse.eye_array(2 * n_pts),
+                        exy2ez))
+    return op
+
+
+def e2h(
+        angular_wavenumber: complex,
+        omega: float,
+        dxes: dx_lists_t,
+        rmin: float,
+        mu: vfdfield_t | None = None
+        ) -> sparse.spmatrix:
+    r"""
+    Returns an operator which, when applied to a vectorized E eigenfield, produces
+     the vectorized H eigenfield.
+
+    This operator is created directly from the initial coordinate-transformed equations:
+    $$
+    \begin{aligned}
+    \imath \omega \epsilon_{xx} E_x &= \hat{\partial}_y H_z + \imath \beta \frac{H_y}{\hat{T}} \\
+    \imath \omega \epsilon_{yy} E_y &= -\imath \beta H_x - \{1}{\tilde{t}_x} \hat{\partial}_x \hat{t}_x} H_z \\
+    \imath \omega \epsilon_{zz} E_z &= \hat{\partial}_x H_y - \hat{\partial}_y H_x \\
+    \end{aligned}
+    $$
+
+    Args:
+        angular_wavenumber: Wavenumber assuming fields have theta-dependence of
+            `exp(-i * angular_wavenumber * theta)`. It should satisfy
+            `operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy`
+        omega: The angular frequency of the system
+        dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
+        rmin: Radius at the left edge of the simulation domain (at minimum 'x')
+        mu: Vectorized magnetic permeability grid (default 1 everywhere)
+
+    Returns:
+        Sparse matrix representation of the operator.
+    """
+    Dfx, Dfy = deriv_forward(dxes[0])
+    Ta, Tb = dxes2T(dxes=dxes, rmin=rmin)
+    Tai = sparse.diags_array(1 / Ta.diagonal())
+    Tbi = sparse.diags_array(1 / Tb.diagonal())
+
+    jB = 1j * angular_wavenumber / rmin
+    op = sparse.block_array([[    None, -jB * Tai,           -Dfy],
+                             [jB * Tbi,      None, Tbi @ Dfx @ Ta],
+                             [     Dfy,      -Dfx,           None]]) / (-1j * omega)
+    if mu is not None:
+        op = sparse.diags_array(1 / mu) @ op
+    return op
+
+
+def dxes2T(
+        dxes: dx_lists_t,
+        rmin: float,
+        ) -> tuple[NDArray[numpy.float64], NDArray[numpy.float64]]:
+    r"""
+    Returns the $T_a$ and $T_b$ diagonal matrices which are used to apply the cylindrical
+      coordinate transformation in various operators.
+
+    Args:
+        dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
+        rmin: Radius at the left edge of the simulation domain (at minimum 'x')
+
+    Returns:
+        Sparse matrix representations of the operators Ta and Tb
+    """
+    ra = rmin + numpy.cumsum(dxes[0][0])                      # Radius at Ey points
+    rb = rmin + dxes[0][0] / 2.0 + numpy.cumsum(dxes[1][0])   # Radius at Ex points
+    ta = ra / rmin
+    tb = rb / rmin
+
+    Ta = sparse.diags_array(vec(ta[:, None].repeat(dxes[0][1].size, axis=1)))
+    Tb = sparse.diags_array(vec(tb[:, None].repeat(dxes[1][1].size, axis=1)))
+    return Ta, Tb
+
+
+def normalized_fields_e(
+        e_xy: ArrayLike,
+        angular_wavenumber: complex,
+        omega: float,
+        dxes: dx_lists_t,
+        rmin: float,
+        epsilon: vfdfield_t,
+        mu: vfdfield_t | None = None,
+        prop_phase: float = 0,
+        ) -> tuple[vcfdfield_t, vcfdfield_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
+        angular_wavenumber: Wavenumber assuming fields have theta-dependence of
+            `exp(-i * angular_wavenumber * theta)`. It should satisfy
+            `operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy`
+        omega: The angular frequency of the system
+        dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
+        rmin: Radius at the left edge of the simulation domain (at minimum 'x')
+        epsilon: Vectorized dielectric constant grid
+        mu: Vectorized magnetic permeability grid (default 1 everywhere)
+        prop_phase: Phase shift `(dz * corrected_wavenumber)` over 1 cell in propagation direction.
+                    Default 0 (continuous propagation direction, i.e. dz->0).
+
+    Returns:
+        `(e, h)`, where each field is vectorized, normalized,
+        and contains all three vector components.
+    """
+    e = exy2e(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon) @ e_xy
+    h = exy2h(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon, mu=mu) @ e_xy
+    e_norm, h_norm = _normalized_fields(e=e, h=h, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon,
+                                        mu=mu, prop_phase=prop_phase)
+    return e_norm, h_norm
+
+
+def _normalized_fields(
+        e: vcfdfield_t,
+        h: vcfdfield_t,
+        omega: complex,
+        dxes: dx_lists_t,
+        rmin: float,
+        epsilon: vfdfield_t,
+        mu: vfdfield_t | None = None,
+        prop_phase: float = 0,
+        ) -> tuple[vcfdfield_t, vcfdfield_t]:
+    h *= -1
+    # TODO documentation for normalized_fields
+    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)]
+
+    # 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_tavg = waveguide_2d.inner_product(e, h, dxes=dxes, prop_phase=prop_phase, conj_h=True).real   # Note, using linear poynting vector
+    assert Sz_tavg > 0, f'Found a mode propagating in the wrong direction! {Sz_tavg=}'
+
+    energy = numpy.real(epsilon * e.conj() * e)
+
+    norm_amplitude = 1 / numpy.sqrt(Sz_tavg)
+    norm_angle = -numpy.angle(e[energy.argmax()])       # Will randomly add a negative sign when mode is symmetric
+
+    # Try to break symmetry to assign a consistent sign [experimental]
+    E_weighted = unvec(e * energy * numpy.exp(1j * norm_angle), shape)
+    sign = numpy.sign(E_weighted[:,
+                                 :max(shape[0] // 2, 1),
+                                 :max(shape[1] // 2, 1)].real.sum())
+    assert sign != 0
+
+    norm_factor = sign * norm_amplitude * numpy.exp(1j * norm_angle)
+
+    print('\nAAA\n', waveguide_2d.inner_product(e, h, dxes, prop_phase=prop_phase))
+    e *= norm_factor
+    h *= norm_factor
+    print(f'{sign=} {norm_amplitude=} {norm_angle=} {prop_phase=}')
+    print(waveguide_2d.inner_product(e, h, dxes, prop_phase=prop_phase))
+
+    return e, h

meanas/fdmath/__init__.py

@@ -1,4 +1,4 @@
-"""
+r"""
 
 Basic discrete calculus for finite difference (fd) simulations.
 
@@ -43,11 +43,11 @@ Scalar derivatives and cell shifts
 ----------------------------------
 
 Define the discrete forward derivative as
- $$ [\\tilde{\\partial}_x f]_{m + \\frac{1}{2}} = \\frac{1}{\\Delta_{x, m}} (f_{m + 1} - f_m) $$
+ $$ [\tilde{\partial}_x f]_{m + \frac{1}{2}} = \frac{1}{\Delta_{x, m}} (f_{m + 1} - f_m) $$
  where $f$ is a function defined at discrete locations on the x-axis (labeled using $m$).
- The value at $m$ occupies a length $\\Delta_{x, m}$ along the x-axis. Note that $m$
+ The value at $m$ occupies a length $\Delta_{x, m}$ along the x-axis. Note that $m$
  is an index along the x-axis, _not_ necessarily an x-coordinate, since each length
- $\\Delta_{x, m}, \\Delta_{x, m+1}, ...$ is independently chosen.
+ $\Delta_{x, m}, \Delta_{x, m+1}, ...$ is independently chosen.
 
 If we treat `f` as a 1D array of values, with the `i`-th value `f[i]` taking up a length `dx[i]`
 along the x-axis, the forward derivative is
@@ -56,13 +56,13 @@ along the x-axis, the forward derivative is
 
 
 Likewise, discrete reverse derivative is
- $$ [\\hat{\\partial}_x f ]_{m - \\frac{1}{2}} = \\frac{1}{\\Delta_{x, m}} (f_{m} - f_{m - 1}) $$
+ $$ [\hat{\partial}_x f ]_{m - \frac{1}{2}} = \frac{1}{\Delta_{x, m}} (f_{m} - f_{m - 1}) $$
  or
 
     deriv_back(f)[i] = (f[i] - f[i - 1]) / dx[i]
 
 The derivatives' values are shifted by a half-cell relative to the original function, and
-will have different cell widths if all the `dx[i]` ( $\\Delta_{x, m}$ ) are not
+will have different cell widths if all the `dx[i]` ( $\Delta_{x, m}$ ) are not
 identical:
 
     [figure: derivatives and cell sizes]
@@ -88,19 +88,19 @@ identical:
 
 In the above figure,
  `f0 =` $f_0$, `f1 =` $f_1$
- `Df0 =` $[\\tilde{\\partial}f]_{0 + \\frac{1}{2}}$
- `Df1 =` $[\\tilde{\\partial}f]_{1 + \\frac{1}{2}}$
- `df0 =` $[\\hat{\\partial}f]_{0 - \\frac{1}{2}}$
+ `Df0 =` $[\tilde{\partial}f]_{0 + \frac{1}{2}}$
+ `Df1 =` $[\tilde{\partial}f]_{1 + \frac{1}{2}}$
+ `df0 =` $[\hat{\partial}f]_{0 - \frac{1}{2}}$
  etc.
 
-The fractional subscript $m + \\frac{1}{2}$ is used to indicate values defined
+The fractional subscript $m + \frac{1}{2}$ is used to indicate values defined
  at shifted locations relative to the original $m$, with corresponding lengths
- $$ \\Delta_{x, m + \\frac{1}{2}} = \\frac{1}{2} * (\\Delta_{x, m} + \\Delta_{x, m + 1}) $$
+ $$ \Delta_{x, m + \frac{1}{2}} = \frac{1}{2} * (\Delta_{x, m} + \Delta_{x, m + 1}) $$
 
-Just as $m$ is not itself an x-coordinate, neither is $m + \\frac{1}{2}$;
+Just as $m$ is not itself an x-coordinate, neither is $m + \frac{1}{2}$;
 carefully note the positions of the various cells in the above figure vs their labels.
 If the positions labeled with $m$ are considered the "base" or "original" grid,
-the positions labeled with $m + \\frac{1}{2}$ are said to lie on a "dual" or
+the positions labeled with $m + \frac{1}{2}$ are said to lie on a "dual" or
 "derived" grid.
 
 For the remainder of the `Discrete calculus` section, all figures will show
@@ -113,12 +113,12 @@ Gradients and fore-vectors
 --------------------------
 
 Expanding to three dimensions, we can define two gradients
-  $$ [\\tilde{\\nabla} f]_{m,n,p} = \\vec{x} [\\tilde{\\partial}_x f]_{m + \\frac{1}{2},n,p} +
-                                    \\vec{y} [\\tilde{\\partial}_y f]_{m,n + \\frac{1}{2},p} +
-                                    \\vec{z} [\\tilde{\\partial}_z f]_{m,n,p + \\frac{1}{2}}  $$
-  $$ [\\hat{\\nabla} f]_{m,n,p} = \\vec{x} [\\hat{\\partial}_x f]_{m + \\frac{1}{2},n,p} +
-                                  \\vec{y} [\\hat{\\partial}_y f]_{m,n + \\frac{1}{2},p} +
-                                  \\vec{z} [\\hat{\\partial}_z f]_{m,n,p + \\frac{1}{2}}  $$
+  $$ [\tilde{\nabla} f]_{m,n,p} = \vec{x} [\tilde{\partial}_x f]_{m + \frac{1}{2},n,p} +
+                                  \vec{y} [\tilde{\partial}_y f]_{m,n + \frac{1}{2},p} +
+                                  \vec{z} [\tilde{\partial}_z f]_{m,n,p + \frac{1}{2}}  $$
+  $$ [\hat{\nabla} f]_{m,n,p} = \vec{x} [\hat{\partial}_x f]_{m + \frac{1}{2},n,p} +
+                                \vec{y} [\hat{\partial}_y f]_{m,n + \frac{1}{2},p} +
+                                \vec{z} [\hat{\partial}_z f]_{m,n,p + \frac{1}{2}}  $$
 
  or
 
@@ -144,12 +144,12 @@ y in y, and z in z.
 
 We call the resulting object a "fore-vector" or "back-vector", depending
 on the direction of the shift. We write it as
-  $$ \\tilde{g}_{m,n,p} = \\vec{x} g^x_{m + \\frac{1}{2},n,p} +
-                          \\vec{y} g^y_{m,n + \\frac{1}{2},p} +
-                          \\vec{z} g^z_{m,n,p + \\frac{1}{2}} $$
-  $$ \\hat{g}_{m,n,p} = \\vec{x} g^x_{m - \\frac{1}{2},n,p} +
-                        \\vec{y} g^y_{m,n - \\frac{1}{2},p} +
-                        \\vec{z} g^z_{m,n,p - \\frac{1}{2}} $$
+  $$ \tilde{g}_{m,n,p} = \vec{x} g^x_{m + \frac{1}{2},n,p} +
+                         \vec{y} g^y_{m,n + \frac{1}{2},p} +
+                         \vec{z} g^z_{m,n,p + \frac{1}{2}} $$
+  $$ \hat{g}_{m,n,p} = \vec{x} g^x_{m - \frac{1}{2},n,p} +
+                       \vec{y} g^y_{m,n - \frac{1}{2},p} +
+                       \vec{z} g^z_{m,n,p - \frac{1}{2}} $$
 
 
     [figure: gradient / fore-vector]
@@ -172,15 +172,15 @@ Divergences
 
 There are also two divergences,
 
-  $$ d_{n,m,p} = [\\tilde{\\nabla} \\cdot \\hat{g}]_{n,m,p}
-               = [\\tilde{\\partial}_x g^x]_{m,n,p} +
-                 [\\tilde{\\partial}_y g^y]_{m,n,p} +
-                 [\\tilde{\\partial}_z g^z]_{m,n,p}   $$
+  $$ d_{n,m,p} = [\tilde{\nabla} \cdot \hat{g}]_{n,m,p}
+               = [\tilde{\partial}_x g^x]_{m,n,p} +
+                 [\tilde{\partial}_y g^y]_{m,n,p} +
+                 [\tilde{\partial}_z g^z]_{m,n,p}   $$
 
-  $$ d_{n,m,p} = [\\hat{\\nabla} \\cdot \\tilde{g}]_{n,m,p}
-               = [\\hat{\\partial}_x g^x]_{m,n,p} +
-                 [\\hat{\\partial}_y g^y]_{m,n,p} +
-                 [\\hat{\\partial}_z g^z]_{m,n,p}  $$
+  $$ d_{n,m,p} = [\hat{\nabla} \cdot \tilde{g}]_{n,m,p}
+               = [\hat{\partial}_x g^x]_{m,n,p} +
+                 [\hat{\partial}_y g^y]_{m,n,p} +
+                 [\hat{\partial}_z g^z]_{m,n,p}  $$
 
  or
 
@@ -203,7 +203,7 @@ where `g = [gx, gy, gz]` is a fore- or back-vector field.
 
 Since we applied the forward divergence to the back-vector (and vice-versa), the resulting scalar value
 is defined at the back-vector's (fore-vector's) location $(m,n,p)$ and not at the locations of its components
-$(m \\pm \\frac{1}{2},n,p)$ etc.
+$(m \pm \frac{1}{2},n,p)$ etc.
 
     [figure: divergence]
                                     ^^
@@ -227,23 +227,23 @@ Curls
 
 The two curls are then
 
-  $$ \\begin{aligned}
-     \\hat{h}_{m + \\frac{1}{2}, n + \\frac{1}{2}, p + \\frac{1}{2}} &= \\\\
-     [\\tilde{\\nabla} \\times \\tilde{g}]_{m + \\frac{1}{2}, n + \\frac{1}{2}, p + \\frac{1}{2}} &=
-        \\vec{x} (\\tilde{\\partial}_y g^z_{m,n,p + \\frac{1}{2}} - \\tilde{\\partial}_z g^y_{m,n + \\frac{1}{2},p}) \\\\
-     &+ \\vec{y} (\\tilde{\\partial}_z g^x_{m + \\frac{1}{2},n,p} - \\tilde{\\partial}_x g^z_{m,n,p + \\frac{1}{2}}) \\\\
-     &+ \\vec{z} (\\tilde{\\partial}_x g^y_{m,n + \\frac{1}{2},p} - \\tilde{\\partial}_y g^z_{m + \\frac{1}{2},n,p})
-     \\end{aligned} $$
+  $$ \begin{aligned}
+     \hat{h}_{m + \frac{1}{2}, n + \frac{1}{2}, p + \frac{1}{2}} &= \\
+     [\tilde{\nabla} \times \tilde{g}]_{m + \frac{1}{2}, n + \frac{1}{2}, p + \frac{1}{2}} &=
+        \vec{x} (\tilde{\partial}_y g^z_{m,n,p + \frac{1}{2}} - \tilde{\partial}_z g^y_{m,n + \frac{1}{2},p}) \\
+     &+ \vec{y} (\tilde{\partial}_z g^x_{m + \frac{1}{2},n,p} - \tilde{\partial}_x g^z_{m,n,p + \frac{1}{2}}) \\
+     &+ \vec{z} (\tilde{\partial}_x g^y_{m,n + \frac{1}{2},p} - \tilde{\partial}_y g^z_{m + \frac{1}{2},n,p})
+     \end{aligned} $$
 
  and
 
-  $$ \\tilde{h}_{m - \\frac{1}{2}, n - \\frac{1}{2}, p - \\frac{1}{2}} =
-     [\\hat{\\nabla} \\times \\hat{g}]_{m - \\frac{1}{2}, n - \\frac{1}{2}, p - \\frac{1}{2}} $$
+  $$ \tilde{h}_{m - \frac{1}{2}, n - \frac{1}{2}, p - \frac{1}{2}} =
+     [\hat{\nabla} \times \hat{g}]_{m - \frac{1}{2}, n - \frac{1}{2}, p - \frac{1}{2}} $$
 
-  where $\\hat{g}$ and $\\tilde{g}$ are located at $(m,n,p)$
-  with components at $(m \\pm \\frac{1}{2},n,p)$ etc.,
-  while $\\hat{h}$ and $\\tilde{h}$ are located at $(m \\pm \\frac{1}{2}, n \\pm \\frac{1}{2}, p \\pm \\frac{1}{2})$
-  with components at $(m, n \\pm \\frac{1}{2}, p \\pm \\frac{1}{2})$ etc.
+  where $\hat{g}$ and $\tilde{g}$ are located at $(m,n,p)$
+  with components at $(m \pm \frac{1}{2},n,p)$ etc.,
+  while $\hat{h}$ and $\tilde{h}$ are located at $(m \pm \frac{1}{2}, n \pm \frac{1}{2}, p \pm \frac{1}{2})$
+  with components at $(m, n \pm \frac{1}{2}, p \pm \frac{1}{2})$ etc.
 
 
     [code: curls]
@@ -287,27 +287,27 @@ Maxwell's Equations
 
 If we discretize both space (m,n,p) and time (l), Maxwell's equations become
 
- $$ \\begin{aligned}
-  \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}} &= -\\tilde{\\partial}_t \\hat{B}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}}
-                                                                         - \\hat{M}_{l, \\vec{r} + \\frac{1}{2}}  \\\\
-  \\hat{\\nabla} \\times \\hat{H}_{l-\\frac{1}{2},\\vec{r} + \\frac{1}{2}} &= \\hat{\\partial}_t \\tilde{D}_{l, \\vec{r}}
-                                                                             + \\tilde{J}_{l-\\frac{1}{2},\\vec{r}} \\\\
-  \\tilde{\\nabla} \\cdot \\hat{B}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}} &= 0 \\\\
-  \\hat{\\nabla} \\cdot \\tilde{D}_{l,\\vec{r}} &= \\rho_{l,\\vec{r}}
- \\end{aligned} $$
+ $$ \begin{aligned}
+  \tilde{\nabla} \times \tilde{E}_{l,\vec{r}} &= -\tilde{\partial}_t \hat{B}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}}
+                                                                   - \hat{M}_{l, \vec{r} + \frac{1}{2}}  \\
+  \hat{\nabla} \times \hat{H}_{l-\frac{1}{2},\vec{r} + \frac{1}{2}} &= \hat{\partial}_t \tilde{D}_{l, \vec{r}}
+                                                                   + \tilde{J}_{l-\frac{1}{2},\vec{r}} \\
+  \tilde{\nabla} \cdot \hat{B}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}} &= 0 \\
+  \hat{\nabla} \cdot \tilde{D}_{l,\vec{r}} &= \rho_{l,\vec{r}}
+ \end{aligned} $$
 
  with
 
- $$ \\begin{aligned}
-  \\hat{B}_{\\vec{r}} &= \\mu_{\\vec{r} + \\frac{1}{2}} \\cdot \\hat{H}_{\\vec{r} + \\frac{1}{2}} \\\\
-  \\tilde{D}_{\\vec{r}} &= \\epsilon_{\\vec{r}} \\cdot \\tilde{E}_{\\vec{r}}
- \\end{aligned} $$
+ $$ \begin{aligned}
+  \hat{B}_{\vec{r}} &= \mu_{\vec{r} + \frac{1}{2}} \cdot \hat{H}_{\vec{r} + \frac{1}{2}} \\
+  \tilde{D}_{\vec{r}} &= \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}}
+ \end{aligned} $$
 
-where the spatial subscripts are abbreviated as $\\vec{r} = (m, n, p)$ and
-$\\vec{r} + \\frac{1}{2} = (m + \\frac{1}{2}, n + \\frac{1}{2}, p + \\frac{1}{2})$,
-$\\tilde{E}$ and $\\hat{H}$ are the electric and magnetic fields,
-$\\tilde{J}$ and $\\hat{M}$ are the electric and magnetic current distributions,
-and $\\epsilon$ and $\\mu$ are the dielectric permittivity and magnetic permeability.
+where the spatial subscripts are abbreviated as $\vec{r} = (m, n, p)$ and
+$\vec{r} + \frac{1}{2} = (m + \frac{1}{2}, n + \frac{1}{2}, p + \frac{1}{2})$,
+$\tilde{E}$ and $\hat{H}$ are the electric and magnetic fields,
+$\tilde{J}$ and $\hat{M}$ are the electric and magnetic current distributions,
+and $\epsilon$ and $\mu$ are the dielectric permittivity and magnetic permeability.
 
 The above is Yee's algorithm, written in a form analogous to Maxwell's equations.
 The time derivatives can be expanded to form the update equations:
@@ -369,34 +369,34 @@ Each component forms its own grid, offset from the others:
 
 The divergence equations can be derived by taking the divergence of the curl equations
 and combining them with charge continuity,
- $$ \\hat{\\nabla} \\cdot \\tilde{J} + \\hat{\\partial}_t \\rho = 0 $$
+ $$ \hat{\nabla} \cdot \tilde{J} + \hat{\partial}_t \rho = 0 $$
  implying that the discrete Maxwell's equations do not produce spurious charges.
 
 
 Wave equation
 -------------
 
-Taking the backward curl of the $\\tilde{\\nabla} \\times \\tilde{E}$ equation and
-replacing the resulting $\\hat{\\nabla} \\times \\hat{H}$ term using its respective equation,
-and setting $\\hat{M}$ to zero, we can form the discrete wave equation:
+Taking the backward curl of the $\tilde{\nabla} \times \tilde{E}$ equation and
+replacing the resulting $\hat{\nabla} \times \hat{H}$ term using its respective equation,
+and setting $\hat{M}$ to zero, we can form the discrete wave equation:
 
 $$
-  \\begin{aligned}
-  \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}} &=
-      -\\tilde{\\partial}_t \\hat{B}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}}
-                          - \\hat{M}_{l-1, \\vec{r} + \\frac{1}{2}}  \\\\
-  \\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}} &=
-    -\\tilde{\\partial}_t \\hat{H}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}}  \\\\
-  \\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}}) &=
-    \\hat{\\nabla} \\times (-\\tilde{\\partial}_t \\hat{H}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}})  \\\\
-  \\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}}) &=
-    -\\tilde{\\partial}_t \\hat{\\nabla} \\times \\hat{H}_{l-\\frac{1}{2}, \\vec{r} + \\frac{1}{2}}  \\\\
-  \\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}}) &=
-    -\\tilde{\\partial}_t \\hat{\\partial}_t \\epsilon_{\\vec{r}} \\tilde{E}_{l, \\vec{r}} + \\hat{\\partial}_t \\tilde{J}_{l-\\frac{1}{2},\\vec{r}} \\\\
-  \\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l,\\vec{r}})
-            + \\tilde{\\partial}_t \\hat{\\partial}_t \\epsilon_{\\vec{r}} \\cdot \\tilde{E}_{l, \\vec{r}}
-            &= \\tilde{\\partial}_t \\tilde{J}_{l - \\frac{1}{2}, \\vec{r}}
-  \\end{aligned}
+  \begin{aligned}
+  \tilde{\nabla} \times \tilde{E}_{l,\vec{r}} &=
+     -\tilde{\partial}_t \hat{B}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}}
+                       - \hat{M}_{l-1, \vec{r} + \frac{1}{2}}  \\
+  \mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}} &=
+   -\tilde{\partial}_t \hat{H}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}}  \\
+  \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}}) &=
+   \hat{\nabla} \times (-\tilde{\partial}_t \hat{H}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}})  \\
+  \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}}) &=
+   -\tilde{\partial}_t \hat{\nabla} \times \hat{H}_{l-\frac{1}{2}, \vec{r} + \frac{1}{2}}  \\
+  \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}}) &=
+   -\tilde{\partial}_t \hat{\partial}_t \epsilon_{\vec{r}} \tilde{E}_{l, \vec{r}} + \hat{\partial}_t \tilde{J}_{l-\frac{1}{2},\vec{r}} \\
+  \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l,\vec{r}})
+           + \tilde{\partial}_t \hat{\partial}_t \epsilon_{\vec{r}} \cdot \tilde{E}_{l, \vec{r}}
+           &= \tilde{\partial}_t \tilde{J}_{l - \frac{1}{2}, \vec{r}}
+  \end{aligned}
 $$
 
 
@@ -406,27 +406,27 @@ Frequency domain
 We can substitute in a time-harmonic fields
 
 $$
- \\begin{aligned}
- \\tilde{E}_{l, \\vec{r}} &= \\tilde{E}_{\\vec{r}} e^{-\\imath \\omega l \\Delta_t} \\\\
- \\tilde{J}_{l, \\vec{r}} &= \\tilde{J}_{\\vec{r}} e^{-\\imath \\omega (l - \\frac{1}{2}) \\Delta_t}
- \\end{aligned}
+ \begin{aligned}
+ \tilde{E}_{l, \vec{r}} &= \tilde{E}_{\vec{r}} e^{-\imath \omega l \Delta_t} \\
+ \tilde{J}_{l, \vec{r}} &= \tilde{J}_{\vec{r}} e^{-\imath \omega (l - \frac{1}{2}) \Delta_t}
+ \end{aligned}
 $$
 
 resulting in
 
 $$
- \\begin{aligned}
- \\tilde{\\partial}_t &\\Rightarrow (e^{ \\imath \\omega \\Delta_t} - 1) / \\Delta_t = \\frac{-2 \\imath}{\\Delta_t} \\sin(\\omega \\Delta_t / 2) e^{-\\imath \\omega \\Delta_t / 2} = -\\imath \\Omega e^{-\\imath \\omega \\Delta_t / 2}\\\\
-   \\hat{\\partial}_t &\\Rightarrow (1 - e^{-\\imath \\omega \\Delta_t}) / \\Delta_t = \\frac{-2 \\imath}{\\Delta_t} \\sin(\\omega \\Delta_t / 2) e^{ \\imath \\omega \\Delta_t / 2} = -\\imath \\Omega e^{ \\imath \\omega \\Delta_t / 2}\\\\
- \\Omega &= 2 \\sin(\\omega \\Delta_t / 2) / \\Delta_t
- \\end{aligned}
+ \begin{aligned}
+ \tilde{\partial}_t &\Rightarrow (e^{ \imath \omega \Delta_t} - 1) / \Delta_t = \frac{-2 \imath}{\Delta_t} \sin(\omega \Delta_t / 2) e^{-\imath \omega \Delta_t / 2} = -\imath \Omega e^{-\imath \omega \Delta_t / 2}\\
+   \hat{\partial}_t &\Rightarrow (1 - e^{-\imath \omega \Delta_t}) / \Delta_t = \frac{-2 \imath}{\Delta_t} \sin(\omega \Delta_t / 2) e^{ \imath \omega \Delta_t / 2} = -\imath \Omega e^{ \imath \omega \Delta_t / 2}\\
+ \Omega &= 2 \sin(\omega \Delta_t / 2) / \Delta_t
+ \end{aligned}
 $$
 
 This gives the frequency-domain wave equation,
 
 $$
- \\hat{\\nabla} \\times (\\mu^{-1}_{\\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{\\vec{r}})
-    -\\Omega^2 \\epsilon_{\\vec{r}} \\cdot \\tilde{E}_{\\vec{r}} = -\\imath \\Omega \\tilde{J}_{\\vec{r}} e^{\\imath \\omega \\Delta_t / 2} \\\\
+ \hat{\nabla} \times (\mu^{-1}_{\vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{\vec{r}})
+    -\Omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} = -\imath \Omega \tilde{J}_{\vec{r}} e^{\imath \omega \Delta_t / 2} \\
 $$
 
 
@@ -436,48 +436,48 @@ Plane waves and Dispersion relation
 With uniform material distribution and no sources
 
 $$
- \\begin{aligned}
- \\mu_{\\vec{r} + \\frac{1}{2}} &= \\mu \\\\
- \\epsilon_{\\vec{r}} &= \\epsilon \\\\
- \\tilde{J}_{\\vec{r}} &= 0 \\\\
- \\end{aligned}
+ \begin{aligned}
+ \mu_{\vec{r} + \frac{1}{2}} &= \mu \\
+ \epsilon_{\vec{r}} &= \epsilon \\
+ \tilde{J}_{\vec{r}} &= 0 \\
+ \end{aligned}
 $$
 
 the frequency domain wave equation simplifies to
 
-$$ \\hat{\\nabla} \\times \\tilde{\\nabla} \\times \\tilde{E}_{\\vec{r}} - \\Omega^2 \\epsilon \\mu \\tilde{E}_{\\vec{r}} = 0 $$
+$$ \hat{\nabla} \times \tilde{\nabla} \times \tilde{E}_{\vec{r}} - \Omega^2 \epsilon \mu \tilde{E}_{\vec{r}} = 0 $$
 
-Since $\\hat{\\nabla} \\cdot \\tilde{E}_{\\vec{r}} = 0$, we can simplify
+Since $\hat{\nabla} \cdot \tilde{E}_{\vec{r}} = 0$, we can simplify
 
 $$
- \\begin{aligned}
- \\hat{\\nabla} \\times \\tilde{\\nabla} \\times \\tilde{E}_{\\vec{r}}
-   &= \\tilde{\\nabla}(\\hat{\\nabla} \\cdot \\tilde{E}_{\\vec{r}}) - \\hat{\\nabla} \\cdot \\tilde{\\nabla} \\tilde{E}_{\\vec{r}} \\\\
-   &= - \\hat{\\nabla} \\cdot \\tilde{\\nabla} \\tilde{E}_{\\vec{r}} \\\\
-   &= - \\tilde{\\nabla}^2 \\tilde{E}_{\\vec{r}}
- \\end{aligned}
+ \begin{aligned}
+ \hat{\nabla} \times \tilde{\nabla} \times \tilde{E}_{\vec{r}}
+  &= \tilde{\nabla}(\hat{\nabla} \cdot \tilde{E}_{\vec{r}}) - \hat{\nabla} \cdot \tilde{\nabla} \tilde{E}_{\vec{r}} \\
+  &= - \hat{\nabla} \cdot \tilde{\nabla} \tilde{E}_{\vec{r}} \\
+  &= - \tilde{\nabla}^2 \tilde{E}_{\vec{r}}
+ \end{aligned}
 $$
 
 and we get
 
-$$  \\tilde{\\nabla}^2 \\tilde{E}_{\\vec{r}} + \\Omega^2 \\epsilon \\mu \\tilde{E}_{\\vec{r}} = 0 $$
+$$  \tilde{\nabla}^2 \tilde{E}_{\vec{r}} + \Omega^2 \epsilon \mu \tilde{E}_{\vec{r}} = 0 $$
 
 We can convert this to three scalar-wave equations of the form
 
-$$ (\\tilde{\\nabla}^2 + K^2) \\phi_{\\vec{r}} = 0 $$
+$$ (\tilde{\nabla}^2 + K^2) \phi_{\vec{r}} = 0 $$
 
-with $K^2 = \\Omega^2 \\mu \\epsilon$. Now we let
+with $K^2 = \Omega^2 \mu \epsilon$. Now we let
 
-$$  \\phi_{\\vec{r}} = A e^{\\imath (k_x m \\Delta_x + k_y n \\Delta_y + k_z p \\Delta_z)}  $$
+$$  \phi_{\vec{r}} = A e^{\imath (k_x m \Delta_x + k_y n \Delta_y + k_z p \Delta_z)}  $$
 
 resulting in
 
 $$
- \\begin{aligned}
- \\tilde{\\partial}_x &\\Rightarrow (e^{ \\imath k_x \\Delta_x} - 1) / \\Delta_t = \\frac{-2 \\imath}{\\Delta_x} \\sin(k_x \\Delta_x / 2) e^{ \\imath k_x \\Delta_x / 2} = \\imath K_x e^{ \\imath k_x \\Delta_x / 2}\\\\
-   \\hat{\\partial}_x &\\Rightarrow (1 - e^{-\\imath k_x \\Delta_x}) / \\Delta_t = \\frac{-2 \\imath}{\\Delta_x} \\sin(k_x \\Delta_x / 2) e^{-\\imath k_x \\Delta_x / 2} = \\imath K_x e^{-\\imath k_x \\Delta_x / 2}\\\\
- K_x &= 2 \\sin(k_x \\Delta_x / 2) / \\Delta_x \\\\
- \\end{aligned}
+ \begin{aligned}
+ \tilde{\partial}_x &\Rightarrow (e^{ \imath k_x \Delta_x} - 1) / \Delta_t = \frac{-2 \imath}{\Delta_x} \sin(k_x \Delta_x / 2) e^{ \imath k_x \Delta_x / 2} = \imath K_x e^{ \imath k_x \Delta_x / 2}\\
+   \hat{\partial}_x &\Rightarrow (1 - e^{-\imath k_x \Delta_x}) / \Delta_t = \frac{-2 \imath}{\Delta_x} \sin(k_x \Delta_x / 2) e^{-\imath k_x \Delta_x / 2} = \imath K_x e^{-\imath k_x \Delta_x / 2}\\
+ K_x &= 2 \sin(k_x \Delta_x / 2) / \Delta_x \\
+ \end{aligned}
 $$
 
 with similar expressions for the y and z dimnsions (and $K_y, K_z$).
@@ -485,20 +485,20 @@ with similar expressions for the y and z dimnsions (and $K_y, K_z$).
 This implies
 
 $$
-  \\tilde{\\nabla}^2 = -(K_x^2 + K_y^2 + K_z^2) \\phi_{\\vec{r}} \\\\
-  K_x^2 + K_y^2 + K_z^2 = \\Omega^2 \\mu \\epsilon = \\Omega^2 / c^2
+  \tilde{\nabla}^2 = -(K_x^2 + K_y^2 + K_z^2) \phi_{\vec{r}} \\
+  K_x^2 + K_y^2 + K_z^2 = \Omega^2 \mu \epsilon = \Omega^2 / c^2
 $$
 
-where $c = \\sqrt{\\mu \\epsilon}$.
+where $c = \sqrt{\mu \epsilon}$.
 
-Assuming real $(k_x, k_y, k_z), \\omega$ will be real only if
+Assuming real $(k_x, k_y, k_z), \omega$ will be real only if
 
-$$ c^2 \\Delta_t^2 = \\frac{\\Delta_t^2}{\\mu \\epsilon} < 1/(\\frac{1}{\\Delta_x^2} + \\frac{1}{\\Delta_y^2} + \\frac{1}{\\Delta_z^2}) $$
+$$ c^2 \Delta_t^2 = \frac{\Delta_t^2}{\mu \epsilon} < 1/(\frac{1}{\Delta_x^2} + \frac{1}{\Delta_y^2} + \frac{1}{\Delta_z^2}) $$
 
-If $\\Delta_x = \\Delta_y = \\Delta_z$, this simplifies to $c \\Delta_t < \\Delta_x / \\sqrt{3}$.
+If $\Delta_x = \Delta_y = \Delta_z$, this simplifies to $c \Delta_t < \Delta_x / \sqrt{3}$.
 This last form can be interpreted as enforcing causality; the distance that light
-travels in one timestep (i.e., $c \\Delta_t$) must be less than the diagonal
-of the smallest cell ( $\\Delta_x / \\sqrt{3}$ when on a uniform cubic grid).
+travels in one timestep (i.e., $c \Delta_t$) must be less than the diagonal
+of the smallest cell ( $\Delta_x / \sqrt{3}$ when on a uniform cubic grid).
 
 
 Grid description
@@ -515,8 +515,8 @@ to make the illustration simpler; we need at least two cells in the x dimension
 demonstrate how nonuniform `dx` affects the various components.
 
 Place the E fore-vectors at integer indices $r = (m, n, p)$ and the H back-vectors
-at fractional indices $r + \\frac{1}{2} = (m + \\frac{1}{2}, n + \\frac{1}{2},
-p + \\frac{1}{2})$. Remember that these are indices and not coordinates; they can
+at fractional indices $r + \frac{1}{2} = (m + \frac{1}{2}, n + \frac{1}{2},
+p + \frac{1}{2})$. Remember that these are indices and not coordinates; they can
 correspond to arbitrary (monotonically increasing) coordinates depending on the cell widths.
 
 Draw lines to denote the planes on which the H components and back-vectors are defined.
@@ -718,14 +718,14 @@ composed of the three diagonal tensor components:
 or
 
 $$
- \\epsilon = \\begin{bmatrix} \\epsilon_{xx} & 0 & 0 \\\\
-                              0 & \\epsilon_{yy} & 0 \\\\
-                              0 & 0 & \\epsilon_{zz} \\end{bmatrix}
+ \epsilon = \begin{bmatrix} \epsilon_{xx} & 0 & 0 \\
+                            0 & \epsilon_{yy} & 0 \\
+                            0 & 0 & \epsilon_{zz} \end{bmatrix}
 $$
 $$
- \\mu = \\begin{bmatrix} \\mu_{xx} & 0 & 0 \\\\
-                         0 & \\mu_{yy} & 0 \\\\
-                         0 & 0 & \\mu_{zz} \\end{bmatrix}
+ \mu = \begin{bmatrix} \mu_{xx} & 0 & 0 \\
+                         0 & \mu_{yy} & 0 \\
+                         0 & 0 & \mu_{zz} \end{bmatrix}
 $$
 
 where the off-diagonal terms (e.g. `epsilon_xy`) are assumed to be zero.
@@ -741,8 +741,24 @@ the true values can be multiplied back in after the simulation is complete if no
 normalized results are needed.
 """
 
-from .types import fdfield_t, vfdfield_t, cfdfield_t, vcfdfield_t, dx_lists_t, dx_lists_mut
-from .types import fdfield_updater_t, cfdfield_updater_t
-from .vectorization import vec, unvec
-from . import operators, functional, types, vectorization
+from .types import (
+    fdfield_t as fdfield_t,
+    vfdfield_t as vfdfield_t,
+    cfdfield_t as cfdfield_t,
+    vcfdfield_t as vcfdfield_t,
+    dx_lists_t as dx_lists_t,
+    dx_lists_mut as dx_lists_mut,
+    fdfield_updater_t as fdfield_updater_t,
+    cfdfield_updater_t as cfdfield_updater_t,
+    )
+from .vectorization import (
+    vec as vec,
+    unvec as unvec,
+    )
+from . import (
+    operators as operators,
+    functional as functional,
+    types as types,
+    vectorization as vectorization,
+    )
 

meanas/fdmath/functional.py

@@ -3,16 +3,18 @@ Math functions for finite difference simulations
 
 Basic discrete calculus etc.
 """
-from typing import Sequence, Callable
+from typing import TypeVar
+from collections.abc import Sequence, Callable
 
 import numpy
 from numpy.typing import NDArray
+from numpy import floating, complexfloating
 
 from .types import fdfield_t, fdfield_updater_t
 
 
 def deriv_forward(
-        dx_e: Sequence[NDArray[numpy.float_]] | None = None,
+        dx_e: Sequence[NDArray[floating | complexfloating]] | None = None,
         ) -> tuple[fdfield_updater_t, fdfield_updater_t, fdfield_updater_t]:
     """
     Utility operators for taking discretized derivatives (backward variant).
@@ -36,7 +38,7 @@ def deriv_forward(
 
 
 def deriv_back(
-        dx_h: Sequence[NDArray[numpy.float_]] | None = None,
+        dx_h: Sequence[NDArray[floating | complexfloating]] | None = None,
         ) -> tuple[fdfield_updater_t, fdfield_updater_t, fdfield_updater_t]:
     """
     Utility operators for taking discretized derivatives (forward variant).
@@ -59,10 +61,13 @@ def deriv_back(
     return derivs
 
 
+TT = TypeVar('TT', bound='NDArray[floating | complexfloating]')
+
+
 def curl_forward(
-        dx_e: Sequence[NDArray[numpy.float_]] | None = None,
-        ) -> fdfield_updater_t:
-    """
+        dx_e: Sequence[NDArray[floating | complexfloating]] | None = None,
+        ) -> Callable[[TT], TT]:
+    r"""
     Curl operator for use with the E field.
 
     Args:
@@ -71,11 +76,11 @@ def curl_forward(
 
     Returns:
         Function `f` for taking the discrete forward curl of a field,
-        `f(E)` -> curlE $= \\nabla_f \\times E$
+        `f(E)` -> curlE $= \nabla_f \times E$
     """
     Dx, Dy, Dz = deriv_forward(dx_e)
 
-    def ce_fun(e: fdfield_t) -> fdfield_t:
+    def ce_fun(e: TT) -> TT:
         output = numpy.empty_like(e)
         output[0] = Dy(e[2])
         output[1] = Dz(e[0])
@@ -89,9 +94,9 @@ def curl_forward(
 
 
 def curl_back(
-        dx_h: Sequence[NDArray[numpy.float_]] | None = None,
-        ) -> fdfield_updater_t:
-    """
+        dx_h: Sequence[NDArray[floating | complexfloating]] | None = None,
+        ) -> Callable[[TT], TT]:
+    r"""
     Create a function which takes the backward curl of a field.
 
     Args:
@@ -100,11 +105,11 @@ def curl_back(
 
     Returns:
         Function `f` for taking the discrete backward curl of a field,
-        `f(H)` -> curlH $= \\nabla_b \\times H$
+        `f(H)` -> curlH $= \nabla_b \times H$
     """
     Dx, Dy, Dz = deriv_back(dx_h)
 
-    def ch_fun(h: fdfield_t) -> fdfield_t:
+    def ch_fun(h: TT) -> TT:
         output = numpy.empty_like(h)
         output[0] = Dy(h[2])
         output[1] = Dz(h[0])
@@ -118,7 +123,7 @@ def curl_back(
 
 
 def curl_forward_parts(
-        dx_e: Sequence[NDArray[numpy.float_]] | None = None,
+        dx_e: Sequence[NDArray[floating | complexfloating]] | None = None,
         ) -> Callable:
     Dx, Dy, Dz = deriv_forward(dx_e)
 
@@ -131,7 +136,7 @@ def curl_forward_parts(
 
 
 def curl_back_parts(
-        dx_h: Sequence[NDArray[numpy.float_]] | None = None,
+        dx_h: Sequence[NDArray[floating | complexfloating]] | None = None,
         ) -> Callable:
     Dx, Dy, Dz = deriv_back(dx_h)
 

meanas/fdmath/operators.py

@@ -3,10 +3,11 @@ Matrix operators for finite difference simulations
 
 Basic discrete calculus etc.
 """
-from typing import Sequence
+from collections.abc import Sequence
 import numpy
 from numpy.typing import NDArray
-import scipy.sparse as sparse   # type: ignore
+from numpy import floating, complexfloating
+from scipy import sparse
 
 from .types import vfdfield_t
 
@@ -29,12 +30,12 @@ def shift_circ(
         Sparse matrix for performing the circular shift.
     """
     if len(shape) not in (2, 3):
-        raise Exception('Invalid shape: {}'.format(shape))
+        raise Exception(f'Invalid shape: {shape}')
     if axis not in range(len(shape)):
-        raise Exception('Invalid direction: {}, shape is {}'.format(axis, shape))
+        raise Exception(f'Invalid direction: {axis}, shape is {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)]
+    shifts = [abs(shift_distance) if a == axis else 0 for a in range(len(shape))]
+    shifted_diags = [(numpy.arange(n) + s) % n for n, s in zip(shape, shifts, strict=True)]
     ijk = numpy.meshgrid(*shifted_diags, indexing='ij')
 
     n = numpy.prod(shape)
@@ -69,12 +70,11 @@ def shift_with_mirror(
         Sparse matrix for performing the shift-with-mirror.
     """
     if len(shape) not in (2, 3):
-        raise Exception('Invalid shape: {}'.format(shape))
+        raise Exception(f'Invalid shape: {shape}')
     if axis not in range(len(shape)):
-        raise Exception('Invalid direction: {}, shape is {}'.format(axis, shape))
+        raise Exception(f'Invalid direction: {axis}, shape is {shape}')
     if shift_distance >= shape[axis]:
-        raise Exception('Shift ({}) is too large for axis {} of size {}'.format(
-                        shift_distance, axis, shape[axis]))
+        raise Exception(f'Shift ({shift_distance}) is too large for axis {axis} of size {shape[axis]}')
 
     def mirrored_range(n: int, s: int) -> NDArray[numpy.int_]:
         v = numpy.arange(n) + s
@@ -82,8 +82,8 @@ def shift_with_mirror(
         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)]
+    shifts = [shift_distance if a == axis else 0 for a in range(len(shape))]
+    shifted_diags = [mirrored_range(n, s) for n, s in zip(shape, shifts, strict=True)]
     ijk = numpy.meshgrid(*shifted_diags, indexing='ij')
 
     n = numpy.prod(shape)
@@ -97,7 +97,7 @@ def shift_with_mirror(
 
 
 def deriv_forward(
-        dx_e: Sequence[NDArray[numpy.float_]],
+        dx_e: Sequence[NDArray[floating | complexfloating]],
         ) -> list[sparse.spmatrix]:
     """
     Utility operators for taking discretized derivatives (forward variant).
@@ -124,7 +124,7 @@ def deriv_forward(
 
 
 def deriv_back(
-        dx_h: Sequence[NDArray[numpy.float_]],
+        dx_h: Sequence[NDArray[floating | complexfloating]],
         ) -> list[sparse.spmatrix]:
     """
     Utility operators for taking discretized derivatives (backward variant).
@@ -198,7 +198,7 @@ def avg_forward(axis: int, shape: Sequence[int]) -> sparse.spmatrix:
         Sparse matrix for forward average operation.
     """
     if len(shape) not in (2, 3):
-        raise Exception('Invalid shape: {}'.format(shape))
+        raise Exception(f'Invalid shape: {shape}')
 
     n = numpy.prod(shape)
     return 0.5 * (sparse.eye(n) + shift_circ(axis, shape))
@@ -219,7 +219,7 @@ def avg_back(axis: int, shape: Sequence[int]) -> sparse.spmatrix:
 
 
 def curl_forward(
-        dx_e: Sequence[NDArray[numpy.float_]],
+        dx_e: Sequence[NDArray[floating | complexfloating]],
         ) -> sparse.spmatrix:
     """
     Curl operator for use with the E field.
@@ -235,7 +235,7 @@ def curl_forward(
 
 
 def curl_back(
-        dx_h: Sequence[NDArray[numpy.float_]],
+        dx_h: Sequence[NDArray[floating | complexfloating]],
         ) -> sparse.spmatrix:
     """
     Curl operator for use with the H field.

meanas/fdmath/types.py

@@ -1,26 +1,26 @@
 """
 Types shared across multiple submodules
 """
-from typing import Sequence, Callable, MutableSequence
-import numpy
+from collections.abc import Sequence, Callable, MutableSequence
 from numpy.typing import NDArray
+from numpy import floating, complexfloating
 
 
 # Field types
-fdfield_t = NDArray[numpy.float_]
+fdfield_t = NDArray[floating]
 """Vector field with shape (3, X, Y, Z) (e.g. `[E_x, E_y, E_z]`)"""
 
-vfdfield_t = NDArray[numpy.float_]
+vfdfield_t = NDArray[floating]
 """Linearized vector field (single vector of length 3*X*Y*Z)"""
 
-cfdfield_t = NDArray[numpy.complex_]
+cfdfield_t = NDArray[complexfloating]
 """Complex vector field with shape (3, X, Y, Z) (e.g. `[E_x, E_y, E_z]`)"""
 
-vcfdfield_t = NDArray[numpy.complex_]
+vcfdfield_t = NDArray[complexfloating]
 """Linearized complex vector field (single vector of length 3*X*Y*Z)"""
 
 
-dx_lists_t = Sequence[Sequence[NDArray[numpy.float_]]]
+dx_lists_t = Sequence[Sequence[NDArray[floating | complexfloating]]]
 """
  'dxes' datastructure which contains grid cell width information in the following format:
 
@@ -31,7 +31,7 @@ dx_lists_t = Sequence[Sequence[NDArray[numpy.float_]]]
    and `dy_h[0]` is the y-width of the `y=0` cells, as used when calculating dH/dy, etc.
 """
 
-dx_lists_mut = MutableSequence[MutableSequence[NDArray[numpy.float_]]]
+dx_lists_mut = MutableSequence[MutableSequence[NDArray[floating | complexfloating]]]
 """Mutable version of `dx_lists_t`"""
 
 

meanas/fdmath/vectorization.py

@@ -4,7 +4,8 @@ and a 1D array representation of that field `[f_x0, f_x1, f_x2,... f_y0,... f_z0
 Vectorized versions of the field use row-major (ie., C-style) ordering.
 """
 
-from typing import overload, Sequence
+from typing import overload
+from collections.abc import Sequence
 import numpy
 from numpy.typing import ArrayLike
 
@@ -27,14 +28,16 @@ def vec(f: cfdfield_t) -> vcfdfield_t:
 def vec(f: ArrayLike) -> vfdfield_t | vcfdfield_t:
     pass
 
-def vec(f: fdfield_t | cfdfield_t | ArrayLike | None) -> vfdfield_t | vcfdfield_t | None:
+def vec(
+        f: fdfield_t | cfdfield_t | ArrayLike | None,
+        ) -> vfdfield_t | vcfdfield_t | None:
     """
-    Create a 1D ndarray from a 3D vector field which spans a 1-3D region.
+    Create a 1D ndarray from a vector field which spans a 1-3D region.
 
     Returns `None` if called with `f=None`.
 
     Args:
-        f: A vector field, `[f_x, f_y, f_z]` where each `f_` component is a 1- to
+        f: A vector field, e.g. `[f_x, f_y, f_z]` where each `f_` component is a 1- to
            3-D ndarray (`f_*` should all be the same size). Doesn't fail with `f=None`.
 
     Returns:
@@ -46,33 +49,38 @@ def vec(f: fdfield_t | cfdfield_t | ArrayLike | None) -> vfdfield_t | vcfdfield_
 
 
 @overload
-def unvec(v: None, shape: Sequence[int]) -> None:
+def unvec(v: None, shape: Sequence[int], nvdim: int = 3) -> None:
     pass
 
 @overload
-def unvec(v: vfdfield_t, shape: Sequence[int]) -> fdfield_t:
+def unvec(v: vfdfield_t, shape: Sequence[int], nvdim: int = 3) -> fdfield_t:
     pass
 
 @overload
-def unvec(v: vcfdfield_t, shape: Sequence[int]) -> cfdfield_t:
+def unvec(v: vcfdfield_t, shape: Sequence[int], nvdim: int = 3) -> cfdfield_t:
     pass
 
-def unvec(v: vfdfield_t | vcfdfield_t | None, shape: Sequence[int]) -> fdfield_t | cfdfield_t | None:
+def unvec(
+        v: vfdfield_t | vcfdfield_t | None,
+        shape: Sequence[int],
+        nvdim: int = 3,
+        ) -> fdfield_t | cfdfield_t | None:
     """
-    Perform the inverse of vec(): take a 1D ndarray and output a 3D field
-     of form `[f_x, f_y, f_z]` where each of `f_*` is a len(shape)-dimensional
+    Perform the inverse of vec(): take a 1D ndarray and output an `nvdim`-component field
+     of form e.g. `[f_x, f_y, f_z]` (`nvdim=3`) where each of `f_*` is a len(shape)-dimensional
      ndarray.
 
     Returns `None` if called with `v=None`.
 
     Args:
-        v: 1D ndarray representing a 3D vector field of shape shape (or None)
+        v: 1D ndarray representing a vector field of shape shape (or None)
         shape: shape of the vector field
+        nvdim: Number of components in each vector
 
     Returns:
         `[f_x, f_y, f_z]` where each `f_` is a `len(shape)` dimensional ndarray (or `None`)
     """
     if v is None:
         return None
-    return v.reshape((3, *shape), order='C')
+    return v.reshape((nvdim, *shape), order='C')
 

meanas/fdtd/__init__.py

@@ -1,4 +1,4 @@
-"""
+r"""
 Utilities for running finite-difference time-domain (FDTD) simulations
 
 See the discussion of `Maxwell's Equations` in `meanas.fdmath` for basic
@@ -11,9 +11,9 @@ Timestep
 From the discussion of "Plane waves and the Dispersion relation" in `meanas.fdmath`,
 we have
 
-$$ c^2 \\Delta_t^2 = \\frac{\\Delta_t^2}{\\mu \\epsilon} < 1/(\\frac{1}{\\Delta_x^2} + \\frac{1}{\\Delta_y^2} + \\frac{1}{\\Delta_z^2}) $$
+$$ c^2 \Delta_t^2 = \frac{\Delta_t^2}{\mu \epsilon} < 1/(\frac{1}{\Delta_x^2} + \frac{1}{\Delta_y^2} + \frac{1}{\Delta_z^2}) $$
 
-or, if $\\Delta_x = \\Delta_y = \\Delta_z$, then $c \\Delta_t < \\frac{\\Delta_x}{\\sqrt{3}}$.
+or, if $\Delta_x = \Delta_y = \Delta_z$, then $c \Delta_t < \frac{\Delta_x}{\sqrt{3}}$.
 
 Based on this, we can set
 
@@ -27,81 +27,81 @@ Poynting Vector and Energy Conservation
 
 Let
 
-$$ \\begin{aligned}
-  \\tilde{S}_{l, l', \\vec{r}} &=& &\\tilde{E}_{l, \\vec{r}} \\otimes \\hat{H}_{l', \\vec{r} + \\frac{1}{2}}  \\\\
-  &=&  &\\vec{x} (\\tilde{E}^y_{l,m+1,n,p} \\hat{H}^z_{l',\\vec{r} + \\frac{1}{2}} - \\tilde{E}^z_{l,m+1,n,p} \\hat{H}^y_{l', \\vec{r} + \\frac{1}{2}}) \\\\
-  & &+ &\\vec{y} (\\tilde{E}^z_{l,m,n+1,p} \\hat{H}^x_{l',\\vec{r} + \\frac{1}{2}} - \\tilde{E}^x_{l,m,n+1,p} \\hat{H}^z_{l', \\vec{r} + \\frac{1}{2}}) \\\\
-  & &+ &\\vec{z} (\\tilde{E}^x_{l,m,n,p+1} \\hat{H}^y_{l',\\vec{r} + \\frac{1}{2}} - \\tilde{E}^y_{l,m,n,p+1} \\hat{H}^z_{l', \\vec{r} + \\frac{1}{2}})
-   \\end{aligned}
+$$ \begin{aligned}
+  \tilde{S}_{l, l', \vec{r}} &=& &\tilde{E}_{l, \vec{r}} \otimes \hat{H}_{l', \vec{r} + \frac{1}{2}}  \\
+  &=&  &\vec{x} (\tilde{E}^y_{l,m+1,n,p} \hat{H}^z_{l',\vec{r} + \frac{1}{2}} - \tilde{E}^z_{l,m+1,n,p} \hat{H}^y_{l', \vec{r} + \frac{1}{2}}) \\
+  & &+ &\vec{y} (\tilde{E}^z_{l,m,n+1,p} \hat{H}^x_{l',\vec{r} + \frac{1}{2}} - \tilde{E}^x_{l,m,n+1,p} \hat{H}^z_{l', \vec{r} + \frac{1}{2}}) \\
+  & &+ &\vec{z} (\tilde{E}^x_{l,m,n,p+1} \hat{H}^y_{l',\vec{r} + \frac{1}{2}} - \tilde{E}^y_{l,m,n,p+1} \hat{H}^z_{l', \vec{r} + \frac{1}{2}})
+   \end{aligned}
 $$
 
-where $\\vec{r} = (m, n, p)$ and $\\otimes$ is a modified cross product
-in which the $\\tilde{E}$ terms are shifted as indicated.
+where $\vec{r} = (m, n, p)$ and $\otimes$ is a modified cross product
+in which the $\tilde{E}$ terms are shifted as indicated.
 
 By taking the divergence and rearranging terms, we can show that
 
 $$
-  \\begin{aligned}
-  \\hat{\\nabla} \\cdot \\tilde{S}_{l, l', \\vec{r}}
-   &= \\hat{\\nabla} \\cdot (\\tilde{E}_{l, \\vec{r}} \\otimes \\hat{H}_{l', \\vec{r} + \\frac{1}{2}})  \\\\
-   &= \\hat{H}_{l', \\vec{r} + \\frac{1}{2}} \\cdot \\tilde{\\nabla} \\times \\tilde{E}_{l, \\vec{r}} -
-      \\tilde{E}_{l, \\vec{r}} \\cdot \\hat{\\nabla} \\times \\hat{H}_{l', \\vec{r} + \\frac{1}{2}} \\\\
-   &= \\hat{H}_{l', \\vec{r} + \\frac{1}{2}} \\cdot
-            (-\\tilde{\\partial}_t \\mu_{\\vec{r} + \\frac{1}{2}} \\hat{H}_{l - \\frac{1}{2}, \\vec{r} + \\frac{1}{2}} -
-                \\hat{M}_{l, \\vec{r} + \\frac{1}{2}}) -
-      \\tilde{E}_{l, \\vec{r}} \\cdot (\\hat{\\partial}_t \\tilde{\\epsilon}_{\\vec{r}} \\tilde{E}_{l'+\\frac{1}{2}, \\vec{r}} +
-                \\tilde{J}_{l', \\vec{r}}) \\\\
-   &= \\hat{H}_{l'} \\cdot (-\\mu / \\Delta_t)(\\hat{H}_{l + \\frac{1}{2}} - \\hat{H}_{l - \\frac{1}{2}}) -
-      \\tilde{E}_l \\cdot (\\epsilon / \\Delta_t )(\\tilde{E}_{l'+\\frac{1}{2}} - \\tilde{E}_{l'-\\frac{1}{2}})
-      - \\hat{H}_{l'} \\cdot \\hat{M}_{l} - \\tilde{E}_l \\cdot \\tilde{J}_{l'} \\\\
-  \\end{aligned}
+  \begin{aligned}
+  \hat{\nabla} \cdot \tilde{S}_{l, l', \vec{r}}
+   &= \hat{\nabla} \cdot (\tilde{E}_{l, \vec{r}} \otimes \hat{H}_{l', \vec{r} + \frac{1}{2}})  \\
+   &= \hat{H}_{l', \vec{r} + \frac{1}{2}} \cdot \tilde{\nabla} \times \tilde{E}_{l, \vec{r}} -
+      \tilde{E}_{l, \vec{r}} \cdot \hat{\nabla} \times \hat{H}_{l', \vec{r} + \frac{1}{2}} \\
+   &= \hat{H}_{l', \vec{r} + \frac{1}{2}} \cdot
+          (-\tilde{\partial}_t \mu_{\vec{r} + \frac{1}{2}} \hat{H}_{l - \frac{1}{2}, \vec{r} + \frac{1}{2}} -
+              \hat{M}_{l, \vec{r} + \frac{1}{2}}) -
+      \tilde{E}_{l, \vec{r}} \cdot (\hat{\partial}_t \tilde{\epsilon}_{\vec{r}} \tilde{E}_{l'+\frac{1}{2}, \vec{r}} +
+              \tilde{J}_{l', \vec{r}}) \\
+   &= \hat{H}_{l'} \cdot (-\mu / \Delta_t)(\hat{H}_{l + \frac{1}{2}} - \hat{H}_{l - \frac{1}{2}}) -
+      \tilde{E}_l \cdot (\epsilon / \Delta_t )(\tilde{E}_{l'+\frac{1}{2}} - \tilde{E}_{l'-\frac{1}{2}})
+      - \hat{H}_{l'} \cdot \hat{M}_{l} - \tilde{E}_l \cdot \tilde{J}_{l'} \\
+  \end{aligned}
 $$
 
 where in the last line the spatial subscripts have been dropped to emphasize
 the time subscripts $l, l'$, i.e.
 
 $$
-  \\begin{aligned}
-  \\tilde{E}_l &= \\tilde{E}_{l, \\vec{r}} \\\\
-  \\hat{H}_l &= \\tilde{H}_{l, \\vec{r} + \\frac{1}{2}}  \\\\
-  \\tilde{\\epsilon} &= \\tilde{\\epsilon}_{\\vec{r}}  \\\\
-  \\end{aligned}
+  \begin{aligned}
+  \tilde{E}_l &= \tilde{E}_{l, \vec{r}} \\
+  \hat{H}_l &= \tilde{H}_{l, \vec{r} + \frac{1}{2}}  \\
+  \tilde{\epsilon} &= \tilde{\epsilon}_{\vec{r}}  \\
+  \end{aligned}
 $$
 
 etc.
-For $l' = l + \\frac{1}{2}$ we get
+For $l' = l + \frac{1}{2}$ we get
 
 $$
-  \\begin{aligned}
-  \\hat{\\nabla} \\cdot \\tilde{S}_{l, l + \\frac{1}{2}}
-   &= \\hat{H}_{l + \\frac{1}{2}} \\cdot
-            (-\\mu / \\Delta_t)(\\hat{H}_{l + \\frac{1}{2}} - \\hat{H}_{l - \\frac{1}{2}}) -
-      \\tilde{E}_l \\cdot (\\epsilon / \\Delta_t)(\\tilde{E}_{l+1} - \\tilde{E}_l)
-      - \\hat{H}_{l'} \\cdot \\hat{M}_l - \\tilde{E}_l \\cdot \\tilde{J}_{l + \\frac{1}{2}} \\\\
-   &= (-\\mu / \\Delta_t)(\\hat{H}^2_{l + \\frac{1}{2}} - \\hat{H}_{l + \\frac{1}{2}} \\cdot \\hat{H}_{l - \\frac{1}{2}}) -
-      (\\epsilon / \\Delta_t)(\\tilde{E}_{l+1} \\cdot \\tilde{E}_l - \\tilde{E}^2_l)
-      - \\hat{H}_{l'} \\cdot \\hat{M}_l - \\tilde{E}_l \\cdot \\tilde{J}_{l + \\frac{1}{2}} \\\\
-   &= -(\\mu \\hat{H}^2_{l + \\frac{1}{2}}
-       +\\epsilon \\tilde{E}_{l+1} \\cdot \\tilde{E}_l) / \\Delta_t \\ \\
-      +(\\mu \\hat{H}_{l + \\frac{1}{2}} \\cdot \\hat{H}_{l - \\frac{1}{2}}
-       +\\epsilon \\tilde{E}^2_l) / \\Delta_t \\ \\
-      - \\hat{H}_{l+\\frac{1}{2}} \\cdot \\hat{M}_l \\ \\
-      - \\tilde{E}_l \\cdot \\tilde{J}_{l+\\frac{1}{2}} \\\\
-  \\end{aligned}
+  \begin{aligned}
+  \hat{\nabla} \cdot \tilde{S}_{l, l + \frac{1}{2}}
+   &= \hat{H}_{l + \frac{1}{2}} \cdot
+            (-\mu / \Delta_t)(\hat{H}_{l + \frac{1}{2}} - \hat{H}_{l - \frac{1}{2}}) -
+      \tilde{E}_l \cdot (\epsilon / \Delta_t)(\tilde{E}_{l+1} - \tilde{E}_l)
+      - \hat{H}_{l'} \cdot \hat{M}_l - \tilde{E}_l \cdot \tilde{J}_{l + \frac{1}{2}} \\
+   &= (-\mu / \Delta_t)(\hat{H}^2_{l + \frac{1}{2}} - \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}}) -
+      (\epsilon / \Delta_t)(\tilde{E}_{l+1} \cdot \tilde{E}_l - \tilde{E}^2_l)
+      - \hat{H}_{l'} \cdot \hat{M}_l - \tilde{E}_l \cdot \tilde{J}_{l + \frac{1}{2}} \\
+   &= -(\mu \hat{H}^2_{l + \frac{1}{2}}
+       +\epsilon \tilde{E}_{l+1} \cdot \tilde{E}_l) / \Delta_t \\
+      +(\mu \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}}
+       +\epsilon \tilde{E}^2_l) / \Delta_t \\
+      - \hat{H}_{l+\frac{1}{2}} \cdot \hat{M}_l \\
+      - \tilde{E}_l \cdot \tilde{J}_{l+\frac{1}{2}} \\
+  \end{aligned}
 $$
 
-and for $l' = l - \\frac{1}{2}$,
+and for $l' = l - \frac{1}{2}$,
 
 $$
-  \\begin{aligned}
-  \\hat{\\nabla} \\cdot \\tilde{S}_{l, l - \\frac{1}{2}}
-   &=  (\\mu \\hat{H}^2_{l - \\frac{1}{2}}
-       +\\epsilon \\tilde{E}_{l-1} \\cdot \\tilde{E}_l) / \\Delta_t \\ \\
-      -(\\mu \\hat{H}_{l + \\frac{1}{2}} \\cdot \\hat{H}_{l - \\frac{1}{2}}
-       +\\epsilon \\tilde{E}^2_l) / \\Delta_t \\ \\
-      - \\hat{H}_{l-\\frac{1}{2}} \\cdot \\hat{M}_l \\ \\
-      - \\tilde{E}_l \\cdot \\tilde{J}_{l-\\frac{1}{2}} \\\\
-  \\end{aligned}
+  \begin{aligned}
+  \hat{\nabla} \cdot \tilde{S}_{l, l - \frac{1}{2}}
+   &=  (\mu \hat{H}^2_{l - \frac{1}{2}}
+       +\epsilon \tilde{E}_{l-1} \cdot \tilde{E}_l) / \Delta_t \\
+      -(\mu \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}}
+       +\epsilon \tilde{E}^2_l) / \Delta_t \\
+      - \hat{H}_{l-\frac{1}{2}} \cdot \hat{M}_l \\
+      - \tilde{E}_l \cdot \tilde{J}_{l-\frac{1}{2}} \\
+  \end{aligned}
 $$
 
 These two results form the discrete time-domain analogue to Poynting's theorem.
@@ -109,25 +109,25 @@ They hint at the expressions for the energy, which can be calculated at the same
 time-index as either the E or H field:
 
 $$
- \\begin{aligned}
- U_l &= \\epsilon \\tilde{E}^2_l + \\mu \\hat{H}_{l + \\frac{1}{2}} \\cdot \\hat{H}_{l - \\frac{1}{2}} \\\\
- U_{l + \\frac{1}{2}} &= \\epsilon \\tilde{E}_l \\cdot \\tilde{E}_{l + 1} + \\mu \\hat{H}^2_{l + \\frac{1}{2}} \\\\
- \\end{aligned}
+ \begin{aligned}
+ U_l &= \epsilon \tilde{E}^2_l + \mu \hat{H}_{l + \frac{1}{2}} \cdot \hat{H}_{l - \frac{1}{2}} \\
+ U_{l + \frac{1}{2}} &= \epsilon \tilde{E}_l \cdot \tilde{E}_{l + 1} + \mu \hat{H}^2_{l + \frac{1}{2}} \\
+ \end{aligned}
 $$
 
 Rewriting the Poynting theorem in terms of the energy expressions,
 
 $$
-  \\begin{aligned}
-  (U_{l+\\frac{1}{2}} - U_l) / \\Delta_t
-   &= -\\hat{\\nabla} \\cdot \\tilde{S}_{l, l + \\frac{1}{2}} \\ \\
-      - \\hat{H}_{l+\\frac{1}{2}} \\cdot \\hat{M}_l \\ \\
-      - \\tilde{E}_l \\cdot \\tilde{J}_{l+\\frac{1}{2}} \\\\
-  (U_l - U_{l-\\frac{1}{2}}) / \\Delta_t
-   &= -\\hat{\\nabla} \\cdot \\tilde{S}_{l, l - \\frac{1}{2}} \\ \\
-      - \\hat{H}_{l-\\frac{1}{2}} \\cdot \\hat{M}_l \\ \\
-      - \\tilde{E}_l \\cdot \\tilde{J}_{l-\\frac{1}{2}} \\\\
- \\end{aligned}
+  \begin{aligned}
+  (U_{l+\frac{1}{2}} - U_l) / \Delta_t
+   &= -\hat{\nabla} \cdot \tilde{S}_{l, l + \frac{1}{2}} \\
+      - \hat{H}_{l+\frac{1}{2}} \cdot \hat{M}_l \\
+      - \tilde{E}_l \cdot \tilde{J}_{l+\frac{1}{2}} \\
+  (U_l - U_{l-\frac{1}{2}}) / \Delta_t
+   &= -\hat{\nabla} \cdot \tilde{S}_{l, l - \frac{1}{2}} \\
+      - \hat{H}_{l-\frac{1}{2}} \cdot \hat{M}_l \\
+      - \tilde{E}_l \cdot \tilde{J}_{l-\frac{1}{2}} \\
+ \end{aligned}
 $$
 
 This result is exact and should practically hold to within numerical precision. No time-
@@ -147,10 +147,10 @@ of the time-domain fields.
 The Ricker wavelet (normalized second derivative of a Gaussian) is commonly used for the pulse
 shape. It can be written
 
-$$ f_r(t) = (1 - \\frac{1}{2} (\\omega (t - \\tau))^2) e^{-(\\frac{\\omega (t - \\tau)}{2})^2} $$
+$$ f_r(t) = (1 - \frac{1}{2} (\omega (t - \tau))^2) e^{-(\frac{\omega (t - \tau)}{2})^2} $$
 
-with $\\tau > \\frac{2 * \\pi}{\\omega}$ as a minimum delay to avoid a discontinuity at
-t=0 (assuming the source is off for t<0 this gives $\\sim 10^{-3}$ error at t=0).
+with $\tau > \frac{2 * \pi}{\omega}$ as a minimum delay to avoid a discontinuity at
+t=0 (assuming the source is off for t<0 this gives $\sim 10^{-3}$ error at t=0).
 
 
 
@@ -159,8 +159,22 @@ Boundary conditions
 # TODO notes about boundaries / PMLs
 """
 
-from .base import maxwell_e, maxwell_h
-from .pml import cpml_params, updates_with_cpml
-from .energy import (poynting, poynting_divergence, energy_hstep, energy_estep,
-                     delta_energy_h2e, delta_energy_j)
-from .boundaries import conducting_boundary
+from .base import (
+    maxwell_e as maxwell_e,
+    maxwell_h as maxwell_h,
+    )
+from .pml import (
+    cpml_params as cpml_params,
+    updates_with_cpml as updates_with_cpml,
+    )
+from .energy import (
+    poynting as poynting,
+    poynting_divergence as poynting_divergence,
+    energy_hstep as energy_hstep,
+    energy_estep as energy_estep,
+    delta_energy_h2e as delta_energy_h2e,
+    delta_energy_j as delta_energy_j,
+    )
+from .boundaries import (
+    conducting_boundary as conducting_boundary,
+    )

meanas/fdtd/boundaries.py

@@ -15,13 +15,17 @@ def conducting_boundary(
         ) -> tuple[fdfield_updater_t, fdfield_updater_t]:
     dirs = [0, 1, 2]
     if direction not in dirs:
-        raise Exception('Invalid direction: {}'.format(direction))
+        raise Exception(f'Invalid direction: {direction}')
     dirs.remove(direction)
     u, v = dirs
 
+    boundary_slice: list[Any]
+    shifted1_slice: list[Any]
+    shifted2_slice: list[Any]
+
     if polarity < 0:
-        boundary_slice = [slice(None)] * 3      # type: list[Any]
-        shifted1_slice = [slice(None)] * 3      # type: list[Any]
+        boundary_slice = [slice(None)] * 3
+        shifted1_slice = [slice(None)] * 3
         boundary_slice[direction] = 0
         shifted1_slice[direction] = 1
 
@@ -42,7 +46,7 @@ def conducting_boundary(
     if polarity > 0:
         boundary_slice = [slice(None)] * 3
         shifted1_slice = [slice(None)] * 3
-        shifted2_slice = [slice(None)] * 3      # type: list[Any]
+        shifted2_slice = [slice(None)] * 3
         boundary_slice[direction] = -1
         shifted1_slice[direction] = -2
         shifted2_slice[direction] = -3
@@ -64,4 +68,4 @@ def conducting_boundary(
 
         return ep, hp
 
-    raise Exception('Bad polarity: {}'.format(polarity))
+    raise Exception(f'Bad polarity: {polarity}')

meanas/fdtd/energy.py

@@ -12,7 +12,7 @@ def poynting(
         h: fdfield_t,
         dxes: dx_lists_t | None = None,
         ) -> fdfield_t:
-    """
+    r"""
     Calculate the poynting vector `S` ($S$).
 
     This is the energy transfer rate (amount of energy `U` per `dt` transferred
@@ -44,16 +44,16 @@ def poynting(
 
     The full relationship is
     $$
-      \\begin{aligned}
-      (U_{l+\\frac{1}{2}} - U_l) / \\Delta_t
-       &= -\\hat{\\nabla} \\cdot \\tilde{S}_{l, l + \\frac{1}{2}} \\ \\
-          - \\hat{H}_{l+\\frac{1}{2}} \\cdot \\hat{M}_l \\ \\
-          - \\tilde{E}_l \\cdot \\tilde{J}_{l+\\frac{1}{2}} \\\\
-      (U_l - U_{l-\\frac{1}{2}}) / \\Delta_t
-       &= -\\hat{\\nabla} \\cdot \\tilde{S}_{l, l - \\frac{1}{2}} \\ \\
-          - \\hat{H}_{l-\\frac{1}{2}} \\cdot \\hat{M}_l \\ \\
-          - \\tilde{E}_l \\cdot \\tilde{J}_{l-\\frac{1}{2}} \\\\
-     \\end{aligned}
+      \begin{aligned}
+      (U_{l+\frac{1}{2}} - U_l) / \Delta_t
+       &= -\hat{\nabla} \cdot \tilde{S}_{l, l + \frac{1}{2}} \\
+          - \hat{H}_{l+\frac{1}{2}} \cdot \hat{M}_l \\
+          - \tilde{E}_l \cdot \tilde{J}_{l+\frac{1}{2}} \\
+      (U_l - U_{l-\frac{1}{2}}) / \Delta_t
+       &= -\hat{\nabla} \cdot \tilde{S}_{l, l - \frac{1}{2}} \\
+          - \hat{H}_{l-\frac{1}{2}} \cdot \hat{M}_l \\
+          - \tilde{E}_l \cdot \tilde{J}_{l-\frac{1}{2}} \\
+     \end{aligned}
     $$
 
     These equalities are exact and should practically hold to within numerical precision.

meanas/fdtd/misc.py

@@ -0,0 +1,167 @@
+from typing import Callable
+from collections.abc import Sequence
+import logging
+
+import numpy
+from numpy.typing import NDArray, ArrayLike
+from numpy import pi
+
+
+logger = logging.getLogger(__name__)
+
+
+pulse_fn_t = Callable[[int | NDArray], tuple[float, float, float]]
+
+
+def gaussian_packet(
+        wl: float,
+        dwl: float,
+        dt: float,
+        turn_on: float = 1e-10,
+        one_sided: bool = False,
+        ) -> tuple[pulse_fn_t, float]:
+    """
+    Gaussian pulse (or gaussian ramp) for FDTD excitation
+
+    exp(-a*t*t) ==> exp(-omega * omega / (4 * a))   [fourier, ignoring leading const.]
+
+    FWHM_time is 2 * sqrt(2 * log(2)) * sqrt(2 / a)
+    FWHM_omega is 2 * sqrt(2 * log(2)) * sqrt(2 * a)  = 4 * sqrt(log(2) * a)
+
+    Args:
+        wl: wavelength
+        dwl: Gaussian's FWHM in wavelength space
+        dt: Timestep
+        turn_on: Max allowable amplitude at t=0
+        one_sided: If `True`, source amplitude never decreases after reaching max
+
+    Returns:
+        Source function: src(timestep) -> (envelope[tt], cos[... * tt], sin[... * tt])
+        Delay: number of initial timesteps for which envelope[tt] will be 0
+    """
+    logger.warning('meanas.fdtd.misc functions are still very WIP!')    # TODO
+    # dt * dw = 4 * ln(2)
+
+    omega = 2 * pi / wl
+    freq = 1 / wl
+    fwhm_omega = dwl * omega * omega / (2 * pi)         # dwl -> d_omega (approx)
+    alpha = (fwhm_omega * fwhm_omega) * numpy.log(2) / 8
+    delay = numpy.sqrt(-numpy.log(turn_on) / alpha)
+    delay = numpy.ceil(delay * freq) / freq     # force delay to integer number of periods to maintain phase
+    logger.info(f'src_time {2 * delay / dt}')
+
+    def source_phasor(ii: int | NDArray) -> tuple[float, float, float]:
+        t0 = ii * dt - delay
+        envelope = numpy.sqrt(numpy.sqrt(2 * alpha / pi)) * numpy.exp(-alpha * t0 * t0)
+
+        if one_sided and t0 > 0:
+            envelope = 1
+
+        cc = numpy.cos(omega * t0)
+        ss = numpy.sin(omega * t0)
+
+        return envelope, cc, ss
+
+    # nrm = numpy.exp(-omega * omega / alpha) / 2
+
+    return source_phasor, delay
+
+
+def ricker_pulse(
+        wl: float,
+        dt: float,
+        turn_on: float = 1e-10,
+        ) -> tuple[pulse_fn_t, float]:
+    """
+    Ricker wavelet (second derivative of a gaussian pulse)
+
+    t0 = ii * dt - delay
+    R = w_peak * t0 / 2
+    f(t) = (1 - 2 * (pi * f_peak * t0) ** 2) * exp(-(pi * f_peak * t0)**2
+         = (1 - (w_peak * t0)**2 / 2 exp(-(w_peak * t0 / 2) **2)
+         = (1 - 2 * R * R) * exp(-R * R)
+
+    # NOTE: don't use cosine/sine for J, just for phasor readout
+
+    Args:
+        wl: wavelength
+        dt: Timestep
+        turn_on: Max allowable amplitude at t=0
+
+    Returns:
+        Source function: src(timestep) -> (envelope[tt], cos[... * tt], sin[... * tt])
+        Delay: number of initial timesteps for which envelope[tt] will be 0
+    """
+    logger.warning('meanas.fdtd.misc functions are still very WIP!')    # TODO
+    omega = 2 * pi / wl
+    freq = 1 / wl
+    r0 = omega / 2
+
+    from scipy.optimize import root_scalar
+    delay_results = root_scalar(lambda xx: (1 - omega * omega * tt * tt / 2) * numpy.exp(-omega * omega / 4 * tt * tt) - turn_on, x0=0, x1=-2 / omega)
+    delay = delay_results.root
+    delay = numpy.ceil(delay * freq) / freq     # force delay to integer number of periods to maintain phase
+
+    def source_phasor(ii: int | NDArray) -> tuple[float, float, float]:
+        t0 = ii * dt - delay
+        rr = omega * t0 / 2
+        ff = (1 - 2 * rr * rr) * numpy.exp(-rr * rr)
+
+        cc = numpy.cos(omega * t0)
+        ss = numpy.sin(omega * t0)
+
+        return ff, cc, ss
+
+    return source_phasor, delay
+
+
+def gaussian_beam(
+        xyz: list[NDArray],
+        center: ArrayLike,
+        waist_radius: float,
+        wl: float,
+        tilt: float = 0,
+        ) -> NDArray[numpy.complex128]:
+    """
+    Gaussian beam
+    (solution to paraxial Helmholtz equation)
+
+    Default (no tilt) corresponds to a beam propagating in the -z direction.
+
+    Args:
+        xyz: List of [[x0, x1, ...], [y0, ...], [z0, ...]] positions specifying grid
+            locations at which the field will be sampled.
+        center: [x, y, z] location of beam waist
+        waist_radius: Beam radius at the waist
+        wl: Wavelength
+        tilt: Rotation around y axis. Default (0) has beam propagating in -z direction.
+    """
+    logger.warning('meanas.fdtd.misc functions are still very WIP!')    # TODO
+    w0 = waist_radius
+    grids = numpy.asarray(numpy.meshgrid(*xyz, indexing='ij'))
+    grids -= numpy.asarray(center)[:, None, None, None]
+
+    rot = numpy.array([
+        [ numpy.cos(tilt), 0, numpy.sin(tilt)],
+        [               0, 1,               0],
+        [-numpy.sin(tilt), 0, numpy.cos(tilt)],
+        ])
+
+    xx, yy, zz = numpy.einsum('ij,jxyz->ixyz', rot, grids)
+    r2 = xx * xx + yy * yy
+    z2 = zz * zz
+
+    zr = pi * w0 * w0 / wl
+    zr2 = zr * zr
+    wz2 = w0 * w0 * (1 + z2 / zr2)
+    wz = numpy.sqrt(wz2)        # == fwhm(z) / sqrt(2 * ln(2))
+
+    kk = 2 * pi / wl
+    Rz = zz * (1 + zr2 / z2)
+    gouy = numpy.arctan(zz / zr)
+
+    gaussian = w0 / wz * numpy.exp(-r2 / wz2) * numpy.exp(1j * (kk * zz + kk * r2 / 2 / Rz - gouy))
+
+    row = gaussian[:, :, gaussian.shape[2] // 2]
+    norm = numpy.sqrt((row * row.conj()).sum())
+    return gaussian / norm

meanas/fdtd/pml.py

@@ -7,7 +7,8 @@ PML implementations
 """
 # TODO retest pmls!
 
-from typing import Callable, Sequence, Any
+from typing import Any
+from collections.abc import Callable, Sequence
 from copy import deepcopy
 import numpy
 from numpy.typing import NDArray, DTypeLike
@@ -33,10 +34,10 @@ def cpml_params(
         ) -> dict[str, Any]:
 
     if axis not in range(3):
-        raise Exception('Invalid axis: {}'.format(axis))
+        raise Exception(f'Invalid axis: {axis}')
 
     if polarity not in (-1, 1):
-        raise Exception('Invalid polarity: {}'.format(polarity))
+        raise Exception(f'Invalid polarity: {polarity}')
 
     if thickness <= 2:
         raise Exception('It would be wise to have a pml with 4+ cells of thickness')
@@ -111,7 +112,7 @@ def updates_with_cpml(
     params_H: list[list[tuple[Any, Any, Any, Any]]] = deepcopy(params_E)
 
     for axis in range(3):
-        for pp, polarity in enumerate((-1, 1)):
+        for pp, _polarity in enumerate((-1, 1)):
             cpml_param = cpml_params[axis][pp]
             if cpml_param is None:
                 psi_E[axis][pp] = (None, None)
@@ -184,7 +185,7 @@ def updates_with_cpml(
     def update_H(
             e: fdfield_t,
             h: fdfield_t,
-            mu: fdfield_t = numpy.ones(3),
+            mu: fdfield_t | tuple[int, int, int] = (1, 1, 1),
             ) -> None:
         dyEx = Dfy(e[0])
         dzEx = Dfz(e[0])

meanas/test/conftest.py

@@ -3,7 +3,8 @@
 Test fixtures
 
 """
-from typing import Iterable, Any
+# ruff: noqa: ARG001
+from typing import Any
 import numpy
 from numpy.typing import NDArray
 import pytest       # type: ignore
@@ -20,18 +21,18 @@ FixtureRequest = Any
                         (5, 5, 5),
                         # (7, 7, 7),
                        ])
-def shape(request: FixtureRequest) -> Iterable[tuple[int, ...]]:
-    yield (3, *request.param)
+def shape(request: FixtureRequest) -> tuple[int, ...]:
+    return (3, *request.param)
 
 
 @pytest.fixture(scope='module', params=[1.0, 1.5])
-def epsilon_bg(request: FixtureRequest) -> Iterable[float]:
-    yield request.param
+def epsilon_bg(request: FixtureRequest) -> float:
+    return request.param
 
 
 @pytest.fixture(scope='module', params=[1.0, 2.5])
-def epsilon_fg(request: FixtureRequest) -> Iterable[float]:
-    yield request.param
+def epsilon_fg(request: FixtureRequest) -> float:
+    return request.param
 
 
 @pytest.fixture(scope='module', params=['center', '000', 'random'])
@@ -40,7 +41,7 @@ def epsilon(
         shape: tuple[int, ...],
         epsilon_bg: float,
         epsilon_fg: float,
-        ) -> Iterable[NDArray[numpy.float64]]:
+        ) -> NDArray[numpy.float64]:
     is3d = (numpy.array(shape) == 1).sum() == 0
     if is3d:
         if request.param == '000':
@@ -60,17 +61,17 @@ def epsilon(
                                   high=max(epsilon_bg, epsilon_fg),
                                   size=shape)
 
-    yield epsilon
+    return epsilon
 
 
 @pytest.fixture(scope='module', params=[1.0])  # 1.5
-def j_mag(request: FixtureRequest) -> Iterable[float]:
-    yield request.param
+def j_mag(request: FixtureRequest) -> float:
+    return request.param
 
 
 @pytest.fixture(scope='module', params=[1.0, 1.5])
-def dx(request: FixtureRequest) -> Iterable[float]:
-    yield request.param
+def dx(request: FixtureRequest) -> float:
+    return request.param
 
 
 @pytest.fixture(scope='module', params=['uniform', 'centerbig'])
@@ -78,7 +79,7 @@ def dxes(
         request: FixtureRequest,
         shape: tuple[int, ...],
         dx: float,
-        ) -> Iterable[list[list[NDArray[numpy.float64]]]]:
+        ) -> list[list[NDArray[numpy.float64]]]:
     if request.param == 'uniform':
         dxes = [[numpy.full(s, dx) for s in shape[1:]] for _ in range(2)]
     elif request.param == 'centerbig':
@@ -90,5 +91,5 @@ def dxes(
         dxe = [PRNG.uniform(low=1.0 * dx, high=1.1 * dx, size=s) for s in shape[1:]]
         dxh = [(d + numpy.roll(d, -1)) / 2 for d in dxe]
         dxes = [dxe, dxh]
-    yield dxes
+    return dxes
 

meanas/test/test_fdfd.py

@@ -1,4 +1,4 @@
-from typing import Iterable
+# ruff: noqa: ARG001
 import dataclasses
 import pytest       # type: ignore
 import numpy
@@ -61,24 +61,24 @@ def test_poynting_planes(sim: 'FDResult') -> None:
 # Also see conftest.py
 
 @pytest.fixture(params=[1 / 1500])
-def omega(request: FixtureRequest) -> Iterable[float]:
-    yield request.param
+def omega(request: FixtureRequest) -> float:
+    return request.param
 
 
 @pytest.fixture(params=[None])
-def pec(request: FixtureRequest) -> Iterable[NDArray[numpy.float64] | None]:
-    yield request.param
+def pec(request: FixtureRequest) -> NDArray[numpy.float64] | None:
+    return request.param
 
 
 @pytest.fixture(params=[None])
-def pmc(request: FixtureRequest) -> Iterable[NDArray[numpy.float64] | None]:
-    yield request.param
+def pmc(request: FixtureRequest) -> NDArray[numpy.float64] | None:
+    return request.param
 
 
 #@pytest.fixture(scope='module',
 #                params=[(25, 5, 5)])
-#def shape(request):
-#    yield (3, *request.param)
+#def shape(request: FixtureRequest):
+#    return (3, *request.param)
 
 
 @pytest.fixture(params=['diag'])        # 'center'
@@ -86,7 +86,7 @@ def j_distribution(
         request: FixtureRequest,
         shape: tuple[int, ...],
         j_mag: float,
-        ) -> Iterable[NDArray[numpy.float64]]:
+        ) -> NDArray[numpy.float64]:
     j = numpy.zeros(shape, dtype=complex)
     center_mask = numpy.zeros(shape, dtype=bool)
     center_mask[:, shape[1] // 2, shape[2] // 2, shape[3] // 2] = True
@@ -96,7 +96,7 @@ def j_distribution(
     elif request.param == 'diag':
         j[numpy.roll(center_mask, [1, 1, 1], axis=(1, 2, 3))] = (1 + 1j) * j_mag
         j[numpy.roll(center_mask, [-1, -1, -1], axis=(1, 2, 3))] = (1 - 1j) * j_mag
-    yield j
+    return j
 
 
 @dataclasses.dataclass()
@@ -145,7 +145,7 @@ def sim(
         omega=omega,
         dxes=dxes,
         epsilon=eps_vec,
-        matrix_solver_opts={'atol': 1e-15, 'tol': 1e-11},
+        matrix_solver_opts={'atol': 1e-15, 'rtol': 1e-11},
         )
     e = unvec(e_vec, shape[1:])
 

meanas/test/test_fdfd_pml.py

@@ -1,4 +1,4 @@
-from typing import Iterable
+# ruff: noqa: ARG001
 import pytest       # type: ignore
 import numpy
 from numpy.typing import NDArray
@@ -44,30 +44,30 @@ def test_pml(sim: FDResult, src_polarity: int) -> None:
 # Also see conftest.py
 
 @pytest.fixture(params=[1 / 1500])
-def omega(request: FixtureRequest) -> Iterable[float]:
-    yield request.param
+def omega(request: FixtureRequest) -> float:
+    return request.param
 
 
 @pytest.fixture(params=[None])
-def pec(request: FixtureRequest) -> Iterable[NDArray[numpy.float64] | None]:
-    yield request.param
+def pec(request: FixtureRequest) -> NDArray[numpy.float64] | None:
+    return request.param
 
 
 @pytest.fixture(params=[None])
-def pmc(request: FixtureRequest) -> Iterable[NDArray[numpy.float64] | None]:
-    yield request.param
+def pmc(request: FixtureRequest) -> NDArray[numpy.float64] | None:
+    return request.param
 
 
 @pytest.fixture(params=[(30, 1, 1),
                         (1, 30, 1),
                         (1, 1, 30)])
-def shape(request: FixtureRequest) -> Iterable[tuple[int, ...]]:
-    yield (3, *request.param)
+def shape(request: FixtureRequest) -> tuple[int, int, int]:
+    return (3, *request.param)
 
 
 @pytest.fixture(params=[+1, -1])
-def src_polarity(request: FixtureRequest) -> Iterable[int]:
-    yield request.param
+def src_polarity(request: FixtureRequest) -> int:
+    return request.param
 
 
 @pytest.fixture()
@@ -78,7 +78,7 @@ def j_distribution(
         dxes: dx_lists_mut,
         omega: float,
         src_polarity: int,
-        ) -> Iterable[NDArray[numpy.complex128]]:
+        ) -> NDArray[numpy.complex128]:
     j = numpy.zeros(shape, dtype=complex)
 
     dim = numpy.where(numpy.array(shape[1:]) > 1)[0][0]    # Propagation axis
@@ -106,7 +106,7 @@ def j_distribution(
 
     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
+    return j
 
 
 @pytest.fixture()
@@ -115,9 +115,9 @@ def epsilon(
         shape: tuple[int, ...],
         epsilon_bg: float,
         epsilon_fg: float,
-        ) -> Iterable[NDArray[numpy.float64]]:
+        ) -> NDArray[numpy.float64]:
     epsilon = numpy.full(shape, epsilon_fg, dtype=float)
-    yield epsilon
+    return epsilon
 
 
 @pytest.fixture(params=['uniform'])
@@ -127,7 +127,7 @@ def dxes(
         dx: float,
         omega: float,
         epsilon_fg: float,
-        ) -> Iterable[list[list[NDArray[numpy.float64]]]]:
+        ) -> list[list[NDArray[numpy.float64]]]:
     if request.param == 'uniform':
         dxes = [[numpy.full(s, dx) for s in shape[1:]] for _ in range(2)]
     dim = numpy.where(numpy.array(shape[1:]) > 1)[0][0]    # Propagation axis
@@ -141,7 +141,7 @@ def dxes(
                 epsilon_effective=epsilon_fg,
                 thickness=10,
                 )
-    yield dxes
+    return dxes
 
 
 @pytest.fixture()
@@ -162,7 +162,7 @@ def sim(
         omega=omega,
         dxes=dxes,
         epsilon=eps_vec,
-        matrix_solver_opts={'atol': 1e-15, 'tol': 1e-11},
+        matrix_solver_opts={'atol': 1e-15, 'rtol': 1e-11},
         )
     e = unvec(e_vec, shape[1:])
 

meanas/test/test_fdtd.py

@@ -1,4 +1,5 @@
-from typing import Iterable, Any
+# ruff: noqa: ARG001
+from typing import Any
 import dataclasses
 import pytest       # type: ignore
 import numpy
@@ -101,7 +102,7 @@ def test_poynting_divergence(sim: 'TDResult') -> None:
 def test_poynting_planes(sim: 'TDResult') -> None:
     mask = (sim.js[0] != 0).any(axis=0)
     if mask.sum() > 1:
-        pytest.skip('test_poynting_planes can only test single point sources, got {}'.format(mask.sum()))
+        pytest.skip(f'test_poynting_planes can only test single point sources, got {mask.sum()}')
 
     args: dict[str, Any] = {
         'dxes': sim.dxes,
@@ -150,8 +151,8 @@ def test_poynting_planes(sim: 'TDResult') -> None:
 
 
 @pytest.fixture(params=[0.3])
-def dt(request: FixtureRequest) -> Iterable[float]:
-    yield request.param
+def dt(request: FixtureRequest) -> float:
+    return request.param
 
 
 @dataclasses.dataclass()
@@ -168,8 +169,8 @@ class TDResult:
 
 
 @pytest.fixture(params=[(0, 4, 8)])  # (0,)
-def j_steps(request: FixtureRequest) -> Iterable[tuple[int, ...]]:
-    yield request.param
+def j_steps(request: FixtureRequest) -> tuple[int, ...]:
+    return request.param
 
 
 @pytest.fixture(params=['center', 'random'])
@@ -177,7 +178,7 @@ def j_distribution(
         request: FixtureRequest,
         shape: tuple[int, ...],
         j_mag: float,
-        ) -> Iterable[NDArray[numpy.float64]]:
+        ) -> NDArray[numpy.float64]:
     j = numpy.zeros(shape)
     if request.param == 'center':
         j[:, shape[1] // 2, shape[2] // 2, shape[3] // 2] = j_mag
@@ -185,7 +186,7 @@ def j_distribution(
         j[:, 0, 0, 0] = j_mag
     elif request.param == 'random':
         j[:] = PRNG.uniform(low=-j_mag, high=j_mag, size=shape)
-    yield j
+    return j
 
 
 @pytest.fixture()
@@ -199,8 +200,7 @@ def sim(
         j_steps: tuple[int, ...],
         ) -> TDResult:
     is3d = (numpy.array(shape) == 1).sum() == 0
-    if is3d:
-        if dt != 0.3:
+    if is3d and dt != 0.3:
         pytest.skip('Skipping dt != 0.3 because test is 3D (for speed)')
 
     sim = TDResult(

meanas/test/utils.py

@@ -1,5 +1,3 @@
-from typing import Any
-
 import numpy
 from numpy.typing import NDArray
 
@@ -10,22 +8,25 @@ PRNG = numpy.random.RandomState(12345)
 def assert_fields_close(
         x: NDArray,
         y: NDArray,
-        *args: Any,
-        **kwargs: Any,
-        ) -> None:
-    numpy.testing.assert_allclose(
-        x, y, verbose=False,            # type: ignore
-        err_msg='Fields did not match:\n{}\n{}'.format(numpy.moveaxis(x, -1, 0),
-                                                       numpy.moveaxis(y, -1, 0)),
         *args,
         **kwargs,
+        ) -> None:
+    x_disp = numpy.moveaxis(x, -1, 0)
+    y_disp = numpy.moveaxis(y, -1, 0)
+    numpy.testing.assert_allclose(
+        x,           # type: ignore
+        y,           # type: ignore
+        *args,
+        verbose=False,
+        err_msg=f'Fields did not match:\n{x_disp}\n{y_disp}',
+        **kwargs,
         )
 
 def assert_close(
         x: NDArray,
         y: NDArray,
-        *args: Any,
-        **kwargs: Any,
+        *args,
+        **kwargs,
         ) -> None:
     numpy.testing.assert_allclose(x, y, *args, **kwargs)
 

pdoc_templates/html_helpers.py

@@ -321,7 +321,6 @@ class _ToMarkdown:
         """Wrap URLs in Python-Markdown-compatible <angle brackets>."""
         return re.sub(r'(?<![<"\'])(\s*)((?:http|ftp)s?://[^>)\s]+)(\s*)', r'\1<\2>\3', text)
 
-import subprocess
 
 class _MathPattern(InlineProcessor):
     NAME = 'pdoc-math'

pyproject.toml

@@ -39,8 +39,8 @@ include = [
     ]
 dynamic = ["version"]
 dependencies = [
-    "numpy~=1.21",
-    "scipy",
+    "numpy>=2.0",
+    "scipy~=1.14",
     ]
 
 
@@ -51,3 +51,48 @@ path = "meanas/__init__.py"
 dev = ["pytest", "pdoc", "gridlock"]
 examples = ["gridlock"]
 test = ["pytest"]
+
+
+[tool.ruff]
+exclude = [
+    ".git",
+    "dist",
+    ]
+line-length = 245
+indent-width = 4
+lint.dummy-variable-rgx = "^(_+|(_+[a-zA-Z0-9_]*[a-zA-Z0-9]+?))$"
+lint.select = [
+    "NPY", "E", "F", "W", "B", "ANN", "UP", "SLOT", "SIM", "LOG",
+    "C4", "ISC", "PIE", "PT", "RET", "TCH", "PTH", "INT",
+    "ARG", "PL", "R", "TRY",
+    "G010", "G101", "G201", "G202",
+    "Q002", "Q003", "Q004",
+    ]
+lint.ignore = [
+    #"ANN001",   # No annotation
+    "ANN002",   # *args
+    "ANN003",   # **kwargs
+    "ANN401",   # Any
+    "ANN101",   # self: Self
+    "SIM108",   # single-line if / else assignment
+    "RET504",   # x=y+z; return x
+    "PIE790",   # unnecessary pass
+    "ISC003",   # non-implicit string concatenation
+    "C408",     # dict(x=y) instead of {'x': y}
+    "PLR09",    # Too many xxx
+    "PLR2004",  # magic number
+    "PLC0414",  # import x as x
+    "TRY003",   # Long exception message
+    "TRY002",   # Exception()
+    ]
+
+
+[[tool.mypy.overrides]]
+module = [
+    "scipy",
+    "scipy.optimize",
+    "scipy.linalg",
+    "scipy.sparse",
+    "scipy.sparse.linalg",
+    ]
+ignore_missing_imports = true