enable numpy 2.0 and recent scipy

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

@@ -6,12 +6,13 @@ See the readme or `import meanas; help(meanas)` for more info.
 
 import pathlib
 
-__version__ = '0.8'
+__version__ = '0.9'
 __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,

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
@@ -200,7 +196,7 @@ def poynting_e_cross_h(dxes: dx_lists_t) -> Callable[[cfdfield_t, cfdfield_t], c
 
     Note:
         If `E` and `H` are peak amplitudes as assumed elsewhere in this code,
-        the time-average of the poynting vector is `<S> = Re(S)/2 = Re(E x H) / 2`.
+        the time-average of the poynting vector is `<S> = Re(S)/2 = Re(E x H*) / 2`.
         The factor of `1/2` can be omitted if root-mean-square quantities are used
         instead.
 

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
@@ -45,14 +45,14 @@ def e_full(
         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,17 +321,18 @@ 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')]
-    Ex, Ey, Ez = [ei * da for ei, da in zip(numpy.split(e, 3), dxag)]
+    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, strict=True))
 
     block_diags = [[ None,     fx @ -Ez, fx @  Ey],
                    [ fy @  Ez, None,     fy @ -Ex],
                    [ fz @ -Ey, fz @  Ex, None]]
     block_matrix = sparse.bmat([[sparse.diags(x) if x is not None else None for x in row]
                                 for row in block_diags])
-    P = block_matrix @ sparse.diags(numpy.concatenate(dxag))
+    P = block_matrix @ sparse.diags(numpy.concatenate(dxbg))
     return P
 
 
@@ -348,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
@@ -43,7 +44,8 @@ def _scipy_qmr(
         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)
+            logger.info(f'Solver residual at iteration {ii} : {cur_norm}')
 
     if 'callback' in kwargs:
         def augmented_callback(xk: ArrayLike) -> None:
@@ -67,12 +69,12 @@ 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,
         ) -> vcfdfield_t:
     """
     Conjugate gradient FDFD solver using CSR sparse matrices.

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,166 +12,166 @@ 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^{-\gamma z} \\
+\vec{H}(x,y,z) &= (\vec{H}_t(x, y) + H_z(x, y)\vec{z}) e^{-\gamma 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 + \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}
 $$
 
 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}
+\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$:
+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}
+\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}
 $$
 
-With a similar approach (but using $\\gamma \\tilde{\\partial}_y$ instead), we can get
+With a similar approach (but using $\gamma \tilde{\partial}_y$ instead), we can get
 
 $$
-\\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}
+\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}
 $$
 
-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
+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
 
 $$
-\\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 \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}
 $$
 
 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} \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}
 $$
 
-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 $\gamma 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} (\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}
 $$
 
 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} (\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}
 $$
 
 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}
+-\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}
 $$
 
 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
+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.
 
 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
+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.
 
 
@@ -182,10 +182,10 @@ from typing import Any
 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
 
 
@@ -419,7 +422,7 @@ def _normalized_fields(
     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)
+    assert Sz_tavg > 0, f'Found a mode propagating in the wrong direction! {Sz_tavg=}'
 
     energy = epsilon * e.conj() * e
 
@@ -717,6 +720,109 @@ 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],
         omega: complex,

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,
         }
         ```

meanas/fdfd/waveguide_cyl.py

@@ -9,9 +9,9 @@ As the z-dependence is known, all the functions in this file assume a 2D grid
 # TODO update module docs
 
 import numpy
-import scipy.sparse as sparse       # type: ignore
+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, cfdfield_t
 from ..fdmath.operators import deriv_forward, deriv_back
 from ..eigensolvers import signed_eigensolve, rayleigh_quotient_iteration
 
@@ -25,6 +25,9 @@ def cylindrical_operator(
     """
     Cylindrical coordinate waveguide operator of the form
 
+    (NOTE: See 10.1364/OL.33.001848)
+    TODO: consider 10.1364/OE.20.021583
+
     TODO
 
     for use with a field vector of the form `[E_r, E_y]`.
@@ -73,8 +76,10 @@ def cylindrical_operator(
 
     omega2 = omega * omega
 
-    op = (omega2 * diag((Tx, Ty)) + pa) @ diag((a0, a1)) + \
+    op = (
+        (omega2 * diag((Tx, Ty)) + pa) @ diag((a0, a1))
         - (sparse.bmat(((None, Ty), (Tx, None))) + pb / omega2) @ diag((b0, b1))
+        )
     return op
 
 

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]] | 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]] | 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]] | 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]] | 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]] | 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]] | 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
+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)]
+    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
@@ -83,7 +83,7 @@ def shift_with_mirror(
         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)]
+    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]],
         ) -> 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]],
         ) -> 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]],
         ) -> 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]],
         ) -> 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]]]
 """
  '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]]]
 """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
 

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/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)
 

nom-eme.py

@@ -1,391 +0,0 @@
-import scipy
-import numpy
-from numpy.typing import ArrayLike, NDArray
-
-
-#from simphony.elements import Model
-#from simphony.netlist import Subcircuit
-#from simphony.simulation import SweepSimulation
-#
-#from matplotlib import pyplot as plt
-#
-#
-#class PeriodicLayer(Model):
-#    def __init__(self, left_modes, right_modes, s_params):
-#        self.left_modes = left_modes
-#        self.right_modes = right_modes
-#        self.left_ports = len(self.left_modes)
-#        self.right_ports = len(self.right_modes)
-#        self.normalize_fields()
-#        self.s_params = s_params
-#
-#    def normalize_fields(self):
-#        for mode in range(len(self.left_modes)):
-#            self.left_modes[mode].normalize()
-#        for mode in range(len(self.right_modes)):
-#            self.right_modes[mode].normalize()
-#
-#
-#class PeriodicEME:
-#    def __init__(self, layers=[], num_periods=1):
-#        self.layers = layers
-#        self.num_periods = num_periods
-#        self.wavelength = wavelength
-#
-#    def propagate(self):
-#        wl = self.wavelength
-#        if not len(self.layers):
-#            raise Exception("Must place layers before propagating")
-#
-#        num_modes = max([l.num_modes for l in self.layers])
-#        iface = InterfaceSingleMode if num_modes == 1 else InterfaceMultiMode
-#
-#        eme = EME(layers=self.layers)
-#        left, right = eme.propagate()
-#        self.single_period = eme.s_matrix
-#
-#        period_layer = PeriodicLayer(left.modes, right.modes, self.single_period)
-#        current_layer = PeriodicLayer(left.modes, right.modes, self.single_period)
-#        interface = iface(right, left)
-#
-#        for _ in range(self.num_periods - 1):
-#            current_layer.s_params = cascade(current_layer, interface, wl)
-#            current_layer.s_params = cascade(current_layer, period_layer, wl)
-#
-#        self.s_params = current_layer.s_params
-#
-#
-#class EME:
-#    def __init__(self, layers=[]):
-#        self.layers = layers
-#        self.wavelength = None
-#
-#    def propagate(self):
-#        layers = self.layers
-#        wl = layers[0].wavelength if self.wavelength is None else self.wavelength
-#        if not len(layers):
-#            raise Exception("Must place layers before propagating")
-#
-#        num_modes = max([l.num_modes for l in layers])
-#        iface = InterfaceSingleMode if num_modes == 1 else InterfaceMultiMode
-#
-#        first_layer = layers[0]
-#        current = Current(wl, first_layer)
-#        interface = iface(first_layer, layers[1])
-#
-#        current.s = cascade(current, interface, wl)
-#        current.right_pins = interface.right_pins
-#
-#        for index in range(1, len(layers) - 1):
-#            layer1 = layers[index]
-#            layer2 = layers[index + 1]
-#            interface = iface(layer1, layer2)
-#
-#            current.s = cascade(current, layer1, wl)
-#            current.right_pins = layer1.right_pins
-#
-#            current.s = cascade(current, interface, wl)
-#            current.right_pins = interface.right_pins
-#
-#        last_layer = layers[-1]
-#        current.s = cascade(current, last_layer, wl)
-#        current.right_pins = last_layer.right_pins
-#
-#        self.s_matrix = current.s
-#        return first_layer, last_layer
-#
-#
-#def stack(sa, sb):
-#    qab = numpy.eye() - sa.r11 @ sb.r11
-#    qba = numpy.eye() - sa.r11 @ sb.r11
-#    #s.t12 = sa.t12 @ numpy.pinv(qab) @ sb.t12
-#    #s.r21 = sa.t12 @ numpy.pinv(qab) @ sb.r22 @ sa.t21 + sa.r22
-#    #s.r12 = sb.t21 @ numpy.pinv(qba) @ sa.r11 @ sb.t12 + sb.r11
-#    #s.t21 = sb.t21 @ numpy.pinv(qba) @ sa.t21
-#    s.t12 = sa.t12 @ numpy.linalg.solve(qab, sb.t12)
-#    s.r21 = sa.t12 @ numpy.linalg.solve(qab, sb.r22 @ sa.t21) + sa.r22
-#    s.r12 = sb.t21 @ numpy.linalg.solve(qba, sa.r11 @ sb.t12) + sb.r11
-#    s.t21 = sb.t21 @ numpy.linalg.solve(qba, sa.t21)
-#    return s
-#
-#
-#def cascade(first, second, wavelength):
-#    circuit = Subcircuit("Device")
-#
-#    circuit.add([(first, "first"), (second, "second")])
-#    for port in range(first.right_ports):
-#        circuit.connect("first", "right" + str(port), "second", "left" + str(port))
-#
-#    simulation = SweepSimulation(circuit, wavelength, wavelength, num=1)
-#    result = simulation.simulate()
-#    return result.s
-#
-#
-#class InterfaceSingleMode(Model):
-#    def __init__(self, layer1, layer2, num_modes=1):
-#        self.num_modes = num_modes
-#        self.num_ports = 2 * num_modes
-#        self.s = self.solve(layer1, layer2, num_modes)
-#
-#    def solve(self, layer1, layer2, num_modes):
-#        nm = num_modes
-#        s = numpy.zeros((2 * nm, 2 * nm), dtype=complex)
-#
-#        for ii, left_mode in enumerate(layer1.modes):
-#            for oo, right_mode in enumerate(layer2.modes):
-#                r, t = get_rt(left_mode, right_mode)
-#                s[     oo, ii] = r
-#                s[nm + oo, ii] = t
-#
-#        for ii, right_mode in enumerate(layer2.modes):
-#            for oo, left_mode in enumerate(layer1.modes):
-#                r, t = get_rt(right_mode, left_mode)
-#                s[     oo, nm + ii] = t
-#                s[nm + oo, nm + ii] = r
-#        return s
-#
-#
-#class InterfaceMultiMode(Model):
-#    def __init__(self, layer1, layer2):
-#        self.s = self.solve(layer1, layer2)
-#
-#    def solve(self, layer1, layer2):
-#        n1p = layer1.num_modes
-#        n2p = layer2.num_modes
-#        num_ports = n1p + n2p
-#        s = numpy.zeros((num_ports, num_ports), dtype=complex)
-#
-#        for l1p in range(n1p):
-#            ts = get_t(l1p, layer1, layer2)
-#            rs = get_r(l1p, layer1, layer2, ts)
-#            s[n1p:, l1p] = ts
-#            s[:n1p, l1p] = rs
-#
-#        for l2p in range(n2p):
-#            ts = get_t(l2p, layer2, layer1)
-#            rs = get_r(l2p, layer2, layer1, ts)
-#            s[:n1p, n1p + l2p] = ts
-#            s[n1p:, n1p + l2p] = rs
-#
-#        return s
-
-
-def get_t(p, left, right):
-    A = numpy.empty(left.num_modes, right.num_modes, dtype=complex)
-    for i in range(left.num_modes):
-        for k in range(right.num_modes):
-            # TODO optimize loop
-            A[i, k] = inner_product(right[k], left[i]) + inner_product(left[i], right[k])
-
-    b = numpy.zeros(left.num_modes)
-    b[p] = 2 * inner_product(left[p], left[p])
-
-    x = numpy.linalg.solve(A, b)
-    # NOTE: `A` does not depend on `p`, so it might make sense to partially precompute
-    # the solution (pinv(A), or LU decomposition?)
-    # Actually solve() can take multiple vectors, so just pass it something with the full diagonal?
-
-    xx = numpy.matmul(numpy.linalg.pinv(A), b)        #TODO verify
-    assert numpy.allclose(xx, x)
-    return x
-
-
-def get_r(p, left, right, t):
-    r = numpy.empty(left.num_modes, dtype=complex)
-    for ii in range(left.num_modes):
-        r[ii] = sum((inner_product(right[kk], left[ii]) - inner_product(left[ii], right[kk])) * t[kk]
-                    for kk in range(right.num_modes)
-                   ) / (2 * inner_product(left[ii], left[ii]))
-    return r
-
-
-def get_rt(left, right):
-    a = 0.5 * (inner_product(left, right) + inner_product(right, left))
-    b = 0.5 * (inner_product(left, right) - inner_product(right, left))
-    t = (a ** 2 - b ** 2) / a
-    r = 1 - t / (a + b)
-    return -r, t
-
-
-def inner_product(left_E, right_H, dxes):
-    #  ExHy' - EyHx'
-    cross_z = left_E[0] * right_H[1].conj() - left_E[1] * right_H[0].conj()
-#    cross_z = numpy.cross(left_E, numpy.conj(right_H), axisa=0, axisb=0, axisc=0)[2]
-    return numpy.trapz(numpy.trapz(cross_z, dxes[0][0]), dxes[0][1]) / 2       # TODO might need cumsum on dxes
-
-
-def propagation_matrix(mode_neffs: ArrayLike, wavelength: float, distance: float):
-    eigenv = numpy.array(mode_neffs, copy=False) * 2 * numpy.pi / wavelength
-    prop_diag = numpy.diag(numpy.exp(distance * 1j * numpy.hstack((eigenv, eigenv))))
-    prop_matrix = numpy.roll(prop_diag, len(eigenv), axis=0)
-    return prop_matrix
-
-
-def connect_s(
-        A: NDArray[numpy.complex128],
-        k: int,
-        B: NDArray[numpy.complex128],
-        l: int,
-        ) -> NDArray[numpy.complex128]:
-    """
-    TODO
-    freq x ...  x n x n
-
-    Based on skrf implementation
-
-    Connect two n-port networks' s-matrices together.
-
-    Specifically, connect port `k` on network `A` to port `l` on network
-    `B`. The resultant network has nports = (A.rank + B.rank-2); first
-    (A.rank - 1) ports are from `A`, remainder are from B.
-
-    Assumes same reference impedance for both `k` and `l`; may need to
-    connect an "impedance mismatch" thru element first!
-
-    Args:
-        A: S-parameter matrix of `A`, shape is fxnxn
-        k: port index on `A` (port indices start from 0)
-        B: S-parameter matrix of `B`, shape is fxnxn
-        l: port index on `B`
-
-    Returns:
-        new S-parameter matrix
-    """
-    if k > A.shape[-1] - 1 or l > B.shape[-1] - 1:
-        raise ValueError("port indices are out of range")
-
-    #C = scipy.sparse.block_diag((A, B), dtype=complex)
-    #return innerconnect_s(C, k, A.shape[0] + l)
-
-    nA = A.shape[-1]
-    nB = B.shape[-1]
-    nC = nA + nB - 2
-    assert numpy.array_equal(A.shape[:-2], B.shape[:-2])
-
-    ll = slice(l, l + 1)
-    kk = slice(k, k + 1)
-
-    denom = 1 - A[..., kk, kk] * B[..., ll, ll]
-    Anew = A + A[..., kk, :] * B[..., ll, ll] * A[..., :, kk] / denom
-    Bnew = A[..., kk, :] * B[..., :, ll] / denom
-    Anew = numpy.delete(Anew, (k,), 1)
-    Anew = numpy.delete(Anew, (k,), 2)
-    Bnew = numpy.delete(Bnew, (l,), 1)
-    Bnew = numpy.delete(Bnew, (l,), 2)
-
-    dtype = (A[0, 0] * B[0, 0]).dtype
-    C = numpy.zeros(tuple(A.shape[:-2]) + (nC, nC), dtype=dtype)
-    C[..., :nA - 1, :nA - 1] = Anew
-    C[..., nA - 1:, nA - 1:] = Bnew
-    return C
-
-
-def innerconnect_s(
-        S: NDArray[numpy.complex128],
-        k: int,
-        l: int,
-        ) -> NDArray[numpy.complex128]:
-    """
-    TODO
-    freq x ...  x n x n
-
-    Based on skrf implementation
-
-
-    Connect two ports of a single n-port network's s-matrix.
-    Specifically, connect port `k`  to port `l` on `S`. This results in
-    a (n-2)-port network.
-
-    Assumes same reference impedance for both `k` and `l`; may need to
-    connect an "impedance mismatch" thru element first!
-
-    Args:
-        S: S-parameter matrix of `S`, shape is fxnxn
-        k: port index on `S` (port indices start from 0)
-        l: port index on `S`
-
-    Returns:
-        new S-parameter matrix
-
-    Notes:
-        - Compton, R.C., "Perspectives in microwave circuit analysis",
-            doi:10.1109/MWSCAS.1989.101955
-        - Filipsson, G., "A New General Computer Algorithm for S-Matrix Calculation
-            of Interconnected Multiports",
-            doi:10.1109/EUMA.1981.332972
-    """
-    if k > S.shape[-1] - 1 or l > S.shape[-1] - 1:
-        raise ValueError("port indices are out of range")
-
-    ll = slice(l, l + 1)
-    kk = slice(k, k + 1)
-
-    mkl = 1 - S[..., kk, ll]
-    mlk = 1 - S[..., ll, kk]
-    C = S + (
-              S[..., kk, :] * S[..., :, l] * mlk
-            + S[..., ll, :] * S[..., :, k] * mkl
-            + S[..., kk, :] * S[..., l, l] * S[..., :, kk]
-            + S[..., ll, :] * S[..., k, k] * S[..., :, ll]
-        ) / (
-            mlk * mkl - S[..., kk, kk] * S[..., ll, ll]
-        )
-
-    # remove connected ports
-    C = numpy.delete(C, (k, l), 1)
-    C = numpy.delete(C, (k, l), 2)
-
-    return C
-
-
-def s2abcd(
-        S: NDArray[numpy.complex128],
-        z0: NDArray[numpy.complex128],
-        ) -> NDArray[numpy.complex128]:
-    assert numpy.array_equal(S.shape[:2] == (2, 2))
-    Z1, Z2 = z0
-    cross = S[0, 1] * S[1, 0]
-
-    T = numpy.empty_like(S, dtype=complex)
-    T[0, 0, :] = (Z1.conj() + S[0, 0] * Z1) * (1 - S[1, 1]) + cross * Z1    # A numerator
-    T[0, 1, :] = (Z1.conj() + S[0, 0] * Z1) * (Z1.conj() + S[1, 1] * Z2) - cross * Z1 * Z2  # B numerator
-    T[1, 0, :] = (1 - S[0, 0]) * (1 - S[1, 1]) - cross                      # C numerator
-    T[1, 1, :] = (1 - S[0, 0]) * (Z2.conj() + S[1, 1] * Z2) + cross * Z2    # D numerator
-    det = 2 * S[1, 0] * numpy.sqrt(Z1.real * Z2.real)
-    T /= det
-    return T
-
-
-def generalize_S(
-        S: NDArray[numpy.complex128],
-        r0: float,
-        z0: NDArray[numpy.complex128],
-        ) -> NDArray[numpy.complex128]:
-    g = (z0 - r0) / (z0 + r0)
-    D = numpy.diag((1 - g) / numpy.abs(1 - g.conj()) * numpy.sqrt(1 - numpy.abs(g * g.conj())))
-    G = numpy.diag(g)
-    U = numpy.eye(S.shape[-1]).reshape((S.ndim - 2) * (1,) + (S.shape[-2], S.shape[-1]))
-    S_gen = pinv(D.conj()) @ (S - G.conj()) @ pinv(U - G @ S) @ D
-    return S_gen
-
-
-def change_R0(
-        S: NDArray[numpy.complex128],
-        r1: float,
-        r2: float,
-        ) -> NDArray[numpy.complex128]:
-    g = (r2 - r1) / (r2 + r1)
-    U = numpy.eye(S.shape[-1]).reshape((S.ndim - 2) * (1,) + (S.shape[-2], S.shape[-1]))
-    G = U * g
-    S_r2 = (S - G) @ pinv(U - G @ S)
-    return S_r2
-
-# Zc = numpy.sqrt(B / C)
-# gamma = numpy.arccosh(A) / L_TL
-# n_eff = -1j * gamma * c_light / (2 * pi * f)
-# n_eff_grp = n_eff + f * diff(n_eff) / diff(f)
-# attenuation = (1 - S[0, 0] * S[0, 0].conj()) / (S[1, 0] * S[1, 0].conj())
-# R = numpy.real(gamma * Zc)
-# C = numpy.real(gamma / Zc)
-# L = numpy.imag(gamma * Zc) / (-1j * 2 * pi * f)
-# G = numpy.imag(gamma / Zc) / (-1j * 2 * pi * f)

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>=1.26",
+    "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

spar.py

@@ -1,387 +0,0 @@
-# Based on scripts from Andy H. va rfcafe
-
-# IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES. VOL 42, NO 2. FEBRUARY 1994
-# Conversions Between S, Z, Y, h, ABCD, and T Parameters which are Valid for Complex Source and Load Impedances
-# Dean A. Frickey, Member, EEE
-# Tables I and II
-
-import numpy
-
-
-def s_to_z(s, z0):
-    """
-    Scattering (S) to Impedance (Z)
-
-    Args:
-        s: The scattering matrix.
-        z0: The port impedances (Ohms).
-
-    Returns:
-        The impedance matrix.
-    """
-    z0c = numpy.conj(z0)
-
-    z = numpy.empty([2, 2], dtype=complex)
-    z[0, 0] = (z0c[0] + s[0, 0] * z0[0]) * (1 - s[1, 1]) + s[0, 1] * s[1, 0] * z0[0]
-    z[0, 1] = 2 * s[0, 1] * numpy.sqrt(z0[0].real * z0[1].real)
-    z[1, 0] = 2 * s[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
-    z[1, 1] = (1 - s[0, 0]) * (z0c[1] + s[1, 1] * z0[1]) + s[0, 1] * s[1, 0] * z0[1]
-
-    z /= (1 - s[0, 0]) * (1 - s[1, 1]) - s[0, 1] * s[1, 0]
-    return z
-
-
-def z_to_s(z, z0):
-    """
-    Impedance (Z) to Scattering (S)
-
-    Args:
-        z: The impedance matrix.
-        z0: The port impedances (Ohms).
-
-    Returns:
-        The scattering matrix.
-    """
-    z0c = numpy.conj(z0)
-
-    s = numpy.empty([2, 2], dtype=complex)
-    s[0, 0] = (z[0, 0] - z0c[0]) * (z[1, 1] + z0[1]) - z[0, 1] * z[1, 0]
-    s[0, 1] = 2 * z[0, 1] * numpy.sqrt(z0[0].real * z0[1].real)
-    s[1, 0] = 2 * z[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
-    s[1, 1] = (z[0, 0] + z0[0]) * (z[1, 1] - z0c[1]) - z[0, 1] * z[1, 0]
-
-    s /= (z[0, 0] + z0[0]) * (z[1, 1] + z0[1]) - z[0, 1] * z[1, 0]
-    return s
-
-
-def s_to_y(s, z0):
-    """
-    Scattering (S) to Admittance (Y)
-
-    Args:
-        s: The scattering matrix.
-        z0: The port impedances (Ohms).
-
-    Returns:
-        The admittance matrix.
-    """
-    z0c = numpy.conj(z0)
-
-    y = numpy.empty([2, 2], dtype=complex)
-    y[0, 0] = (1 - s[0, 0]) * (z0c[1] + s[1, 1] * z0[1]) + s[0,1] * s[1, 0] * z0[1]
-    y[0, 1] = -2 * s[0,1] * numpy.sqrt(z0[0].real * z0[1].real)
-    y[1, 0] = -2 * s[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
-    y[1, 1] = (z0c[0] + s[0, 0] * z0[0]) * (1 - s[1,1]) + s[0,1] * s[1, 0] * z0[0]
-
-    y /= (z0c[0] + s[0, 0] * z0[0]) * (z0c[1] + s[1, 1] * z0[1]) - s[0,1] * s[1, 0] * z0[0] * z0[1]
-    return y
-
-
-def y_to_s(y, z0):
-    """
-    Admittance (Y) to Scattering (S)
-
-    Args:
-        y: The admittance matrix.
-        z0: The port impedances (Ohms).
-
-    Returns:
-        The scattering matrix.
-    """
-    z0c = numpy.conj(z0)
-
-    s = numpy.empty([2, 2], dtype=complex)
-    s[0, 0] = (1 - y[0, 0] * z0c[0]) * (1 + y[1, 1] * z0[1]) + y[0,1] * y[1, 0] * z0c[0] * z0[1]
-    s[0, 1] = -2 * y[0,1] * numpy.sqrt(z0[0].real * z0[1].real)
-    s[1, 0] = -2 * y[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
-    s[1, 1] = (1 + y[0, 0] * z0[0]) * (1 - y[1,1] * z0c[1]) + y[0,1] * y[1, 0] * z0[0] * z0c[1]
-
-    s /= (1 + y[0, 0] * z0[0]) * (1 + y[1, 1] * z0[1]) - y[0,1] * y[1, 0] * z0[0] * z0[1]
-    return s
-
-
-def s_to_h(s, z0):
-    """
-    Scattering (S) to Hybrid (H)
-
-    Args:
-        s: The scattering matrix.
-        z0: The port impedances (Ohms).
-
-    Returns:
-        The hybrid matrix.
-    """
-    z0c = numpy.conj(z0)
-
-    h = numpy.empty([2, 2], dtype=complex)
-    h[0, 0] = (z0c[0] + s[0, 0] * z0[0]) * (z0c[1] + s[1, 1] * z0[1]) - s[0,1] * s[1, 0] * z0[0] * z0[1]
-    h[0, 1] = 2 * s[0,1] * numpy.sqrt(z0[0].real * z0[1].real)
-    h[1, 0] = -2 * s[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
-    h[1, 1] = (1 - s[0, 0]) * (1 - s[1,1]) - s[0,1] * s[1, 0]
-
-    h /= (1 - s[0, 0]) * (z0c[1] + s[1, 1] * z0[1]) + s[0,1] * s[1, 0] * z0[1]
-    return h
-
-
-def h_to_s(h, z0):
-    """
-    Hybrid (H) to Scattering (S)
-
-    Args:
-        h: The hybrid matrix.
-        z0: The port impedances (Ohms).
-
-    Returns:
-        The scattering matrix.
-    """
-    z0c = numpy.conj(z0)
-
-    s = numpy.empty([2, 2], dtype=complex)
-    s[0, 0] = (h[0, 0] - z0c[0]) * (1 + h[1, 1] * z0[1]) - h[0,1] * h[1, 0] * z0[1]
-    s[0, 1] = 2 * h[0,1] * numpy.sqrt(z0[0].real * z0[1].real)
-    s[1, 0] = -2 * h[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
-    s[1, 1] = (z0[0] + h[0, 0]) * (1 - h[1,1] * z0c[1]) + h[0,1] * h[1, 0] * z0c[1]
-
-    s /= (z0[0] + h[0, 0]) * (1 + h[1, 1] * z0[1]) - h[0,1] * h[1, 0] * z0[1]
-    return s
-
-
-def s_to_abcd(s, z0):
-    """
-    Scattering to Chain (ABCD)
-
-    Args:
-        s: The scattering matrix.
-        z0: The port impedances (Ohms).
-
-    Returns:
-        The chain matrix.
-    """
-    z0c = numpy.conj(z0)
-
-    ans = numpy.empty([2, 2], dtype=complex)
-    ans[0, 0] = (z0c[0] + s[0, 0] * z0[0]) * (1 - s[1, 1]) + s[0,1] * s[1, 0] * z0[0]
-    ans[0, 1] = (z0c[0] + s[0, 0] * z0[0]) * (z0c[1] + s[1,1] * z0[1]) - s[0,1] * s[1, 0] * z0[0] * z0[1]
-    ans[1, 0] = (1 - s[0, 0]) * (1 - s[1, 1]) - s[0,1] * s[1, 0]
-    ans[1, 1] = (1 - s[0, 0]) * (z0c[1] + s[1,1] * z0[1]) + s[0,1] * s[1, 0] * z0[1]
-
-    ans /= 2 * s[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
-    return ans
-
-
-def abcd_to_s(abcd, z0):
-    """
-    Chain (ABCD) to Scattering (S)
-
-    Args:
-        abcd: The chain matrix.
-        z0: The port impedances (Ohms).
-
-    Return:
-        The scattering matrix.
-    """
-    A = abcd[0, 0]
-    B = abcd[0, 1]
-    C = abcd[1, 0]
-    D = abcd[1, 1]
-
-    z0c = numpy.conj(z0)
-
-    s = numpy.empty([2, 2], dtype=complex)
-    s[0, 0] = A * z0[1] + B - C * z0c[0] * z0[1] - D * z0c[0]
-    s[0, 1] = 2 * (A * D - B * C) * numpy.sqrt(z0[0].real * z0[1].real)
-    s[1, 0] = 2 * numpy.sqrt(z0[0].real * z0[1].real)
-    s[1, 1] = -A * z0c[1] + B - C * z0[0] * z0c[1] + D * z0[0]
-
-    s /= A * z0[1] + B + C * z0[0] * z0[1] + D * z0[0]
-    return s
-
-
-def t_to_z(t, z0):
-    """
-    Chain Transfer (T) to Impedance (Z)
-
-    Args:
-        t: The chain transfer matrix.
-        z0: The port impedances (Ohms).
-
-    Returns:
-        The impedance matrix.
-    """
-    z0c = numpy.conj(z0)
-
-    z = numpy.empty([2, 2], dtype=complex)
-    z[0, 0] = z0c[0] * (t[0, 0] + t[0, 1]) + z0[0] * (t[1, 0] + t[1,1])
-    z[0, 1] = 2 * numpy.sqrt(z0[0].real * z0[1].real) * (t[0, 0] * t[1,1] - t[0,1] * t[1, 0])
-    z[1, 0] = 2 * numpy.sqrt(z0[0].real * z0[1].real)
-    z[1, 1] = z0c[1] * (t[0, 0] - t[1, 0]) - z0[1] * (t[0,1] - t[1,1])
-
-    z /= t[0, 0] + t[0, 1] - t[1, 0] - t[1,1]
-    return z
-
-
-def z_to_t(z, z0):
-    """
-    Impedance (Z) to Chain Transfer (T)
-
-    Args:
-        z: The impedance matrix.
-        z0: The port impedances (Ohms).
-
-    Returns:
-        The chain transfer matrix.
-    """
-    z0c = numpy.conj(z0)
-
-    t = numpy.empty([2, 2], dtype=complex)
-    t[0, 0] = (z[0, 0] + z0[0]) * (z[1, 1] + z0[1]) - z[0,1] * z[1, 0]
-    t[0, 1] = (z[0, 0] + z0[0]) * (z0c[1] - z[1,1]) + z[0,1] * z[1, 0]
-    t[1, 0] = (z[0, 0] - z0c[0]) * (z[1, 1] + z0[1]) - z[0,1] * z[1, 0]
-    t[1, 1] = (z0c[0] - z[0, 0]) * (z[1,1] - z0c[1]) + z[0,1] * z[1, 0]
-
-    t /= 2 * z[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
-    return t
-
-
-def t_to_y(t, z0):
-    """
-    Chain Transfer (T) to Admittance (Y)
-
-    Args:
-        t: The chain transfer matrix.
-        z0: The port impedances (Ohms).
-
-    Returns:
-        The admittance matrix.
-    """
-    z0c = numpy.conj(z0)
-
-    y = numpy.empty([2, 2], dtype=complex)
-    y[0, 0] = z0c[1] * (t[0, 0] - t[1, 0]) - z0[1] * (t[0, 1] - t[1,1])
-    y[0, 1] = -2 * numpy.sqrt(z0[0].real * z0[1].real) * (t[0, 0] * t[1,1] - t[0,1] * t[1, 0])
-    y[1, 0] = -2 * numpy.sqrt(z0[0].real * z0[1].real)
-    y[1, 1] = z0c[0] * (t[0, 0] + t[0,1]) + z0[0] * (t[1, 0] + t[1,1])
-
-    y /= t[0, 0] * z0c[0] * z0c[1] - t[0, 1] * z0c[0] * z0[1] + t[1, 0] * z0[0] * z0c[1] - t[1,1] * z0[0] * z0[1]
-    return y
-
-
-def y_to_t(y, z0):
-    """
-    Admittance (Y) to Chain Transfer (T)
-
-    Args:
-        y: The admittance matrix.
-        z0: The port impedances (Ohms).
-
-    Returns:
-        The chain transfer matrix.
-    """
-    z0c = numpy.conj(z0)
-
-    t = numpy.empty([2, 2], dtype=complex)
-    t[0, 0] = (-1 - y[0, 0] * z0[0]) * (1 + y[1, 1] * z0[1]) + y[0,1] * y[1, 0] * z0[0] * z0[1]
-    t[0, 1] = (1 + y[0, 0] * z0[0]) * (1 - y[1,1] * z0c[1]) + y[0,1] * y[1, 0] * z0[0] * z0c[1]
-    t[1, 0] = (y[0, 0] * z0c[0] - 1) * (1 + y[1, 1] * z0[1]) - y[0,1] * y[1, 0] * z0c[0] * z0[1]
-    t[1, 1] = (1 - y[0, 0] * z0c[0]) * (1 - y[1,1] * z0c[1]) - y[0,1] * y[1, 0] * z0c[0] * z0c[1]
-
-    t /= 2 * y[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
-    return t
-
-
-def t_to_h(t, z0):
-    """
-    Chain Transfer (T) to Hybrid (H)
-
-    Args:
-        t: The chain transfer matrix.
-        z0: The port impedances (Ohms).
-
-    Returns:
-        The hybrid matrix.
-    """
-    z0c = numpy.conj(z0)
-
-
-    h = numpy.empty([2, 2], dtype=complex)
-    h[0, 0] = z0c[1]*(t[0, 0] * z0c[0] + t[1, 0] * z0[0]) - z0[1] * (t[0, 1] * z0c[0] + t[1,1] * z0[0])
-    h[0, 1] = 2 * numpy.sqrt(z0[0].real * z0[1].real) * (t[0, 0] * t[1,1] - t[0,1] * t[1, 0])
-    h[1, 0] = -2 * numpy.sqrt(z0[0].real * z0[1].real)
-    h[1, 1] = t[0, 0] + t[0,1] - t[1, 0] - t[1,1]
-
-    h /= z0c[1] * (t[0, 0] - t[1, 0]) - z0[1] * (t[0, 1] - t[1,1])
-    return h
-
-
-def h_to_t(h, z0):
-    """
-    Hybrid (H) to Chain Transfer (T)
-
-    Args:
-        t: The hybrid matrix.
-        z0: The port impedances (Ohms).
-
-    Returns:
-        The chain transfer matrix.
-    """
-    z0c = numpy.conj(z0)
-
-    t = numpy.empty([2, 2], dtype=complex)
-    t[0, 0] = (-h[0, 0] - z0[0]) * (1 + h[1, 1] * z0[1]) + h[0,1] * h[1, 0] * z0[1]
-    t[0, 1] = (h[0, 0] + z0[0]) * (1 - h[1,1] * z0c[1]) + h[0,1] * h[1, 0] * z0c[1]
-    t[1, 0] = (z0c[0] - h[0, 0]) * (1 + h[1, 1] * z0[1]) + h[0,1] * h[1, 0] * z0[1]
-    t[1, 1] = (h[0, 0] - z0c[0]) * (1 - h[1,1] * z0c[1]) + h[0,1] * h[1, 0] * z0c[1]
-
-    t /= 2 * h[1, 0] * numpy.sqrt(z0[0].real * z0[1].real)
-    return t
-
-
-def t_to_abcd(t, z0):
-    """
-    Chain Transfer (T) to Chain (ABCD)
-
-    Args:
-        t: The chain transfer matrix.
-        z0: The port impedances (Ohms).
-
-    Returns:
-        The chain matrix.
-    """
-    z0c = numpy.conj(z0)
-    ans = numpy.empty([2, 2], dtype=complex)
-    ans[0, 0] = z0c[0] * (t[0, 0] + t[0, 1]) + z0[0] * (t[1, 0] + t[1, 1])
-    ans[0, 1] = z0c[1] * (t[0, 0] * z0c[0] + t[1, 0] * z0[0]) - z0[1] * (t[0, 1] * z0c[0] + t[1, 1] * z0[0])
-    ans[1, 0] = t[0, 0] + t[0, 1] - t[1, 0] - t[1, 1]
-    ans[1, 1] = z0c[1] * (t[0, 0] - t[1, 0]) - z0[1] * (t[0, 1] - t[1, 1])
-
-    ans /= 2 * numpy.sqrt(z0[0].real * z0[1].real)
-    return ans
-
-
-def abcd_to_t(abcd, z0):
-    """
-    Chain (ABCD) to Chain Transfer (T)
-
-    Args:
-        abcd: The chain matrix.
-        z0: The port impedances (Ohms).
-
-    Returns:
-        The chain transfer matrix.
-    """
-    # Break out the components
-    A = abcd[0, 0]
-    B = abcd[0, 1]
-    C = abcd[1, 0]
-    D = abcd[1, 1]
-
-    z0c = numpy.conj(z0)
-
-    t = numpy.empty([2, 2], dtype=complex)
-    t[0, 0] = A * z0[1] + B + C * z0[0] * z0[1] + D * z0[0]
-    t[0, 1] = A * z0c[1] - B + C * z0[0] * z0c[1] - D * z0[0]
-    t[1, 0] = A * z0[1] + B - C * z0c[0] * z0[1] - D * z0c[0]
-    t[1, 1] = A * z0c[1] - B - C * z0c[0] * z0c[1] + D * z0c[0]
-
-    t /= 2 * numpy.sqrt(z0[0].real * z0[1].real)
-    return t

spar_tests.py

@@ -1,58 +0,0 @@
-# IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES. VOL 42, NO 2. FEBRUARY 1994
-# Conversions Between S, Z, Y, h, ABCD, and T Parameters which are Valid for Complex Source and Load Impedances
-# Dean A. Frickey, Member, EEE
-
-# Tables I and II
-
-import numpy as np
-from two_port_conversions import *
-
-""" Testing """
-# Analog Devices - HMC455LP3 - S parameters
-# High output IP3 GaAs InGaP Heterojunction Bipolar Transistor
-
-# MHz S (Magntidue and Angle (deg))
-# 1487.273 0.409 160.117 4.367 163.864 0.063 115.967 0.254 -132.654
-
-s11 = 0.409 * np.exp(1j * np.radians(160.117))
-s12 = 4.367 * np.exp(1j * np.radians(163.864))
-s21 = 0.063 * np.exp(1j * np.radians(115.967))
-s22 = 0.254 * np.exp(1j * np.radians(-132.654))
-
-s_orig = np.array([[s11, s12], [s21, s22]])
-
-# Data specified at 50 Ohms (adding small complex component to test conversions)
-z1, z2 = 50 + 0.01j, 50 - 0.02j
-z0 = np.array([z1, z2])
-
-""" Conversions """
-print(f'Original S: \n{s_orig}\n')
-
-# S --> Z --> T --> Z --> S
-z = s_to_z(s_orig, z0)
-t = z_to_t(z, z0)
-z = t_to_z(t, z0)
-s = z_to_s(z, z0)
-print(f'Test (S --> Z --> T --> Z --> S): \n{s}\n')
-
-# S --> Y --> T --> Y --> S
-y = s_to_y(s_orig, z0)
-t = y_to_t(y, z0)
-y = t_to_y(t, z0)
-s = y_to_s(y, z0)
-print(f'Test (S --> Y --> T --> Y --> S): \n{s}\n')
-
-# S --> H --> T --> H --> S
-h = s_to_h(s_orig, z0)
-t = h_to_t(h, z0)
-h = t_to_h(t, z0)
-s = h_to_s(h, z0)
-print(f'Test (S --> H --> T --> H --> S): \n{s}\n')
-
-# S --> ABCD --> T --> ABCD --> S
-abcd = s_to_abcd(s_orig, z0)
-t = abcd_to_t(abcd, z0)
-abcd = t_to_abcd(t, z0)
-s = abcd_to_s(abcd, z0)
-print(f'Test (S --> ABCD --> T --> ABCD --> S): \n{s}\n')
-

spconv.py

@@ -1,268 +0,0 @@
-import scipy
-import numpy
-from numpy.typing import ArrayLike, NDArray
-from numpy.linalg import pinv
-from numpy import sqrt, real, abs, pi
-
-
-def diag(twod):
-    # numpy.diag() but makes a stack of diagonal matrices
-    return numpy.einsum('ij,jk->ijk', twod, numpy.eye(twod.shape[-1], dtype=twod.dtype))
-
-
-def s2z(s, zref):
-    # G0_inv @ inv(I - S) @ (S Z0 + Z0*) @ G0
-    #  Where Z0 is diag(zref) and G0 = diag(1/sqrt(abs(real(zref))))
-    nf = s.shape[-1]
-    I = numpy.eye(nf)[None, :, :]
-    zref = numpy.array(zref, copy=False)
-    gref = 1 / sqrt(abs(zref.real))
-    z = diag(1 / gref) @ pinv(I - s) @ ( S @ diag(zref) + diag(zref).conj()) @ diag(gref)
-    return z
-
-
-def change_of_zref(
-        s,  # (nf, np, np)
-        zref_old,  # (nf, np)
-        zref_new,  # (nf, np)
-        ):
-    # Change-of-Z0 to Z0'
-    # S' = inv(A) @ (S - rho*) @ inv(I - rho @ S) @ A*
-    #  A = inv(G0') @ G0 @ inv(I - rho*)   (diagonal)
-    #  rho = (Z0' - Z0) @ inv(Z0' + Z0)    (diagonal)
-
-    I = numpy.zeros_like(SL)
-    numpy.einsum('...jj->...j', I)[...] = 1
-
-    zref_old = numpy.array(zref_old, copy=False)
-    zref_new = numpy.array(zref_new, copy=False)
-
-    g_old = 1 / sqrt(abs(zref_old.real))
-    r_new = sqrt(abs(zref_new.real))
-
-    rhov = (zref_new - zref_old) / (zref_new + zref_old)
-    av = r_new * g_old / (1 - rhov.conj())
-
-    s_new = diag(1 / av) @ (s - diag(rhov.conj())) @ pinv(I[None, :] - diag(rhov) @ S) @ diag(av.conj())
-    return s_new
-
-
-def embedding(
-        See,  # (nf, ne, ne)
-        Sei,  # (nf, ne, ni)
-        Sie,  # (nf, ni, ne)
-        Sii,  # (nf, ni, ni)
-        SL,   # (nf, ni, ni)
-        ):
-    # Reveyrand, doi:10.1109/INMMIC.2018.8430023
-
-    I = numpy.zeros_like(SL)
-    numpy.einsum('...jj->...j', I)[...] = 1
-
-    S_tot = See + Sei @ pinv(I - SL @ Sii) @ SL @ Sie
-    return S_tot
-
-
-def deembedding(
-        Sei,  # (nf, ne, ni)
-        Sie,  # (nf, ni, ne)
-        Sext, # (nf, ne, ne)
-        See,  # (nf, ne, ne)
-        Si,   # (nf, ni, ni)
-        ):
-    # Reveyrand, doi:10.1109/INMMIC.2018.8430023
-    # Requires balanced number of ports, similar to VNA calibration
-    Sei_inv = pinv(Sei)
-    Sdif = Sext - See
-    return Sei_inv @ Sdif @ pinv(Sie + Sii @ Sei_inv @ Sdif)
-
-
-def thru_with_Zref_change(
-        zref0,  # (nf,)
-        zref1,  # (nf,)
-        ):
-    s = numpy.empty(tuple(zref0.shape) + (2, 2), dtype=complex)
-    s[..., 0, 0] = zref1 - zref0
-    s[..., 0, 1] = 2 * sqrt(zref0 * zref1)
-    s[..., 1, 0] = s[..., 0, 1]
-    s[..., 1, 1] = -s[..., 0, 0]
-
-    s /= zref0 + zref1
-    return s
-
-
-def propagation_matrix(mode_neffs: ArrayLike, wavelength: float, distance: float):
-    eigenv = numpy.array(mode_neffs, copy=False) * 2 * pi / wavelength
-    prop_diag = numpy.diag(numpy.exp(distance * 1j * numpy.hstack((eigenv, eigenv))))
-    prop_matrix = numpy.roll(prop_diag, len(eigenv), axis=0)
-    return prop_matrix
-
-
-def connect_s(
-        A: NDArray[numpy.complex128],
-        k: int,
-        B: NDArray[numpy.complex128],
-        l: int,
-        ) -> NDArray[numpy.complex128]:
-    """
-    TODO
-    freq x ...  x n x n
-
-    Based on skrf implementation
-
-    Connect two n-port networks' s-matrices together.
-
-    Specifically, connect port `k` on network `A` to port `l` on network
-    `B`. The resultant network has nports = (A.rank + B.rank-2); first
-    (A.rank - 1) ports are from `A`, remainder are from B.
-
-    Assumes same reference impedance for both `k` and `l`; may need to
-    connect an "impedance mismatch" thru element first!
-
-    Args:
-        A: S-parameter matrix of `A`, shape is fxnxn
-        k: port index on `A` (port indices start from 0)
-        B: S-parameter matrix of `B`, shape is fxnxn
-        l: port index on `B`
-
-    Returns:
-        new S-parameter matrix
-    """
-    if k > A.shape[-1] - 1 or l > B.shape[-1] - 1:
-        raise ValueError("port indices are out of range")
-
-    #C = scipy.sparse.block_diag((A, B), dtype=complex)
-    #return innerconnect_s(C, k, A.shape[0] + l)
-
-    nA = A.shape[-1]
-    nB = B.shape[-1]
-    nC = nA + nB - 2
-    assert numpy.array_equal(A.shape[:-2], B.shape[:-2])
-
-    ll = slice(l, l + 1)
-    kk = slice(k, k + 1)
-
-    denom = 1 - A[..., kk, kk] * B[..., ll, ll]
-    Anew = A + A[..., kk, :] * B[..., ll, ll] * A[..., :, kk] / denom
-    Bnew = A[..., kk, :] * B[..., :, ll] / denom
-    Anew = numpy.delete(Anew, (k,), 1)
-    Anew = numpy.delete(Anew, (k,), 2)
-    Bnew = numpy.delete(Bnew, (l,), 1)
-    Bnew = numpy.delete(Bnew, (l,), 2)
-
-    dtype = (A[0, 0] * B[0, 0]).dtype
-    C = numpy.zeros(tuple(A.shape[:-2]) + (nC, nC), dtype=dtype)
-    C[..., :nA - 1, :nA - 1] = Anew
-    C[..., nA - 1:, nA - 1:] = Bnew
-    return C
-
-
-def innerconnect_s(
-        S: NDArray[numpy.complex128],
-        k: int,
-        l: int,
-        ) -> NDArray[numpy.complex128]:
-    """
-    TODO
-    freq x ...  x n x n
-
-    Based on skrf implementation
-
-
-    Connect two ports of a single n-port network's s-matrix.
-    Specifically, connect port `k`  to port `l` on `S`. This results in
-    a (n-2)-port network.
-
-    Assumes same reference impedance for both `k` and `l`; may need to
-    connect an "impedance mismatch" thru element first!
-
-    Args:
-        S: S-parameter matrix of `S`, shape is fxnxn
-        k: port index on `S` (port indices start from 0)
-        l: port index on `S`
-
-    Returns:
-        new S-parameter matrix
-
-    Notes:
-        - Compton, R.C., "Perspectives in microwave circuit analysis",
-            doi:10.1109/MWSCAS.1989.101955
-        - Filipsson, G., "A New General Computer Algorithm for S-Matrix Calculation
-            of Interconnected Multiports",
-            doi:10.1109/EUMA.1981.332972
-    """
-    if k > S.shape[-1] - 1 or l > S.shape[-1] - 1:
-        raise ValueError("port indices are out of range")
-
-    ll = slice(l, l + 1)
-    kk = slice(k, k + 1)
-
-    mkl = 1 - S[..., kk, ll]
-    mlk = 1 - S[..., ll, kk]
-    C = S + (
-              S[..., kk, :] * S[..., :, l] * mlk
-            + S[..., ll, :] * S[..., :, k] * mkl
-            + S[..., kk, :] * S[..., l, l] * S[..., :, kk]
-            + S[..., ll, :] * S[..., k, k] * S[..., :, ll]
-        ) / (
-            mlk * mkl - S[..., kk, kk] * S[..., ll, ll]
-        )
-
-    # remove connected ports
-    C = numpy.delete(C, (k, l), 1)
-    C = numpy.delete(C, (k, l), 2)
-
-    return C
-
-
-def s2abcd(
-        S: NDArray[numpy.complex128],
-        z0: NDArray[numpy.complex128],
-        ) -> NDArray[numpy.complex128]:
-    assert numpy.array_equal(S.shape[:2] == (2, 2))
-    Z1, Z2 = z0
-    cross = S[0, 1] * S[1, 0]
-
-    T = numpy.empty_like(S, dtype=complex)
-    T[0, 0, :] = (Z1.conj() + S[0, 0] * Z1) * (1 - S[1, 1]) + cross * Z1    # A numerator
-    T[0, 1, :] = (Z1.conj() + S[0, 0] * Z1) * (Z1.conj() + S[1, 1] * Z2) - cross * Z1 * Z2  # B numerator
-    T[1, 0, :] = (1 - S[0, 0]) * (1 - S[1, 1]) - cross                      # C numerator
-    T[1, 1, :] = (1 - S[0, 0]) * (Z2.conj() + S[1, 1] * Z2) + cross * Z2    # D numerator
-    det = 2 * S[1, 0] * numpy.sqrt(Z1.real * Z2.real)
-    T /= det
-    return T
-
-
-def generalize_S(
-        S: NDArray[numpy.complex128],
-        r0: float,
-        z0: NDArray[numpy.complex128],
-        ) -> NDArray[numpy.complex128]:
-    g = (z0 - r0) / (z0 + r0)
-    D = numpy.diag((1 - g) / numpy.abs(1 - g.conj()) * numpy.sqrt(1 - numpy.abs(g * g.conj())))
-    G = numpy.diag(g)
-    U = numpy.eye(S.shape[-1]).reshape((S.ndim - 2) * (1,) + (S.shape[-2], S.shape[-1]))
-    S_gen = pinv(D.conj()) @ (S - G.conj()) @ pinv(U - G @ S) @ D
-    return S_gen
-
-
-def change_R0(
-        S: NDArray[numpy.complex128],
-        r1: float,
-        r2: float,
-        ) -> NDArray[numpy.complex128]:
-    g = (r2 - r1) / (r2 + r1)
-    U = numpy.eye(S.shape[-1]).reshape((S.ndim - 2) * (1,) + (S.shape[-2], S.shape[-1]))
-    G = U * g
-    S_r2 = (S - G) @ pinv(U - G @ S)
-    return S_r2
-
-# Zc = numpy.sqrt(B / C)
-# gamma = numpy.arccosh(A) / L_TL
-# n_eff = -1j * gamma * c_light / (2 * pi * f)
-# n_eff_grp = n_eff + f * diff(n_eff) / diff(f)
-# attenuation = (1 - S[0, 0] * S[0, 0].conj()) / (S[1, 0] * S[1, 0].conj())
-# R = numpy.real(gamma * Zc)
-# C = numpy.real(gamma / Zc)
-# L = numpy.imag(gamma * Zc) / (-1j * 2 * pi * f)
-# G = numpy.imag(gamma / Zc) / (-1j * 2 * pi * f)