[docs] publish MkDocs site from Forgejo Actions

.forgejo/workflows/docs.yml

@@ -0,0 +1,57 @@
+name: Publish Docs
+
+on:
+  push:
+    branches:
+      - master
+    paths:
+      - ".forgejo/workflows/docs.yml"
+      - "README.md"
+      - "make_docs.sh"
+      - "mkdocs.yml"
+      - "pyproject.toml"
+      - "docs/**"
+      - "meanas/**"
+  workflow_dispatch:
+
+jobs:
+  publish-docs:
+    runs-on: docker
+    container:
+      image: python:3.13-bookworm
+    env:
+      DOCS_SITE_URL: ${{ vars.DOCS_SITE_URL }}
+    steps:
+      - name: Check out the repository
+        uses: https://data.forgejo.org/actions/checkout@v4
+
+      - name: Install build dependencies
+        run: |
+          apt-get update
+          apt-get install -y --no-install-recommends git
+
+      - name: Install docs dependencies
+        run: |
+          pip install -e '.[docs]'
+
+      - name: Build documentation
+        run: |
+          ./make_docs.sh
+
+      - name: Publish docs branch
+        run: |
+          ./scripts/publish_docs_branch.sh site docs-site
+
+      - name: Write job summary
+        run: |
+          {
+            echo "## Published docs"
+            echo
+            echo "- Branch: \`docs-site\`"
+            if [[ -n "${DOCS_SITE_URL:-}" ]]; then
+              echo "- URL: ${DOCS_SITE_URL}"
+            else
+              echo "- URL: set the \`DOCS_SITE_URL\` repository variable to advertise the published site"
+            fi
+            echo "- Recommended repository setting: configure the Wiki tab to point at the published docs URL"
+          } >> "$GITHUB_STEP_SUMMARY"

.gitignore

@@ -54,6 +54,10 @@ coverage.xml
 
 # documentation
 doc/
+site/
+_doc_mathimg/
+doc.md
+doc.htex
 
 # PyBuilder
 target/

README.md

@@ -94,6 +94,99 @@ python3 -m pytest -rsxX | tee test_results.txt
 
 ## Use
 
-See `examples/` for some simple examples; you may need additional
+`meanas` is organized around a few core workflows:
-packages such as [gridlock](https://mpxd.net/code/jan/gridlock)
+
-to run the examples.
+- `meanas.fdfd`: frequency-domain wave equations, sparse operators, SCPML, and
+  iterative solves for driven problems.
+- `meanas.fdfd.waveguide_2d` / `meanas.fdfd.waveguide_3d`: waveguide mode
+  solvers, mode-source construction, and overlap windows for port-based
+  excitation and analysis.
+- `meanas.fdtd`: Yee-step updates, CPML boundaries, flux/energy accounting, and
+  on-the-fly phasor extraction for comparing time-domain runs against FDFD.
+- `meanas.fdmath`: low-level finite-difference operators, vectorization helpers,
+  and derivations shared by the FDTD and FDFD layers.
+
+The most mature user-facing workflows are:
+
+1. Build an FDFD operator or waveguide port source, then solve a driven
+   frequency-domain problem.
+2. Run an FDTD simulation, extract one or more frequency-domain phasors with
+   `meanas.fdtd.accumulate_phasor(...)`, and compare those phasors against an
+   FDFD reference on the same Yee grid.
+
+## Documentation
+
+API and workflow docs are generated from the package docstrings with
+[MkDocs](https://www.mkdocs.org/), [Material for MkDocs](https://squidfunk.github.io/mkdocs-material/),
+and [mkdocstrings](https://mkdocstrings.github.io/).
+
+When hosted on a Forgejo instance, the intended setup is:
+
+- publish the generated site from a dedicated `docs-site` branch
+- serve that branch from the instance's static-pages host
+- point the repository's **Wiki** tab at the published docs URL
+
+The repository contains a Forgejo Actions workflow for publishing the docs
+branch automatically. Set the repository variable `DOCS_SITE_URL` to the final
+published URL so MkDocs can generate canonical links correctly.
+
+Install the docs toolchain with:
+
+```bash
+pip3 install -e './meanas[docs]'
+```
+
+Then build the docs site with:
+
+```bash
+./make_docs.sh
+```
+
+This produces:
+
+- a normal multi-page site under `site/`
+- a combined printable single-page HTML site under `site/print_page/`
+- an optional fully inlined `site/standalone.html` when `htmlark` is available
+
+The same build output is what the Forgejo Actions workflow publishes to the
+`docs-site` branch.
+
+The docs build uses a local MathJax bundle vendored under `docs/assets/`, so
+the rendered HTML does not rely on external services for equation rendering.
+
+Tracked examples under `examples/` are the intended starting points:
+
+- `examples/fdtd.py`: broadband FDTD pulse excitation, phasor extraction, and a
+  residual check against the matching FDFD operator.
+- `examples/waveguide.py`: waveguide mode solving, unidirectional mode-source
+  construction, overlap readout, and FDTD/FDFD comparison on a guided structure.
+- `examples/fdfd.py`: direct frequency-domain waveguide excitation and overlap /
+  Poynting analysis without a time-domain run.
+
+Several examples rely on optional packages such as
+[gridlock](https://mpxd.net/code/jan/gridlock).
+
+### Frequency-domain waveguide workflow
+
+For a structure with a constant cross-section in one direction:
+
+1. Build `dxes` and the diagonal `epsilon` / `mu` distributions on the Yee grid.
+2. Solve the port mode with `meanas.fdfd.waveguide_3d.solve_mode(...)`.
+3. Build a unidirectional source with `compute_source(...)`.
+4. Build a matching overlap window with `compute_overlap_e(...)`.
+5. Solve the full FDFD problem and project the result onto the overlap window or
+   evaluate plane flux with `meanas.fdfd.functional.poynting_e_cross_h(...)`.
+
+### Time-domain phasor workflow
+
+For a broadband or continuous-wave FDTD run:
+
+1. Advance the fields with `meanas.fdtd.maxwell_e/maxwell_h` or
+   `updates_with_cpml(...)`.
+2. Inject electric current using the same sign convention used throughout the
+   examples and library: `E -= dt * J / epsilon`.
+3. Accumulate the desired phasor with `accumulate_phasor(...)` or the Yee-aware
+   wrappers `accumulate_phasor_e/h/j(...)`.
+4. Build the matching FDFD operator on the stretched `dxes` if CPML/SCPML is
+   part of the simulation, and compare the extracted phasor to the FDFD field or
+   residual.

docs/api/eigensolvers.md

@@ -0,0 +1,3 @@
+# eigensolvers
+
+::: meanas.eigensolvers

docs/api/fdfd.md

@@ -0,0 +1,15 @@
+# fdfd
+
+::: meanas.fdfd
+
+## Core operator layers
+
+::: meanas.fdfd.functional
+
+::: meanas.fdfd.operators
+
+::: meanas.fdfd.solvers
+
+::: meanas.fdfd.scpml
+
+::: meanas.fdfd.farfield

docs/api/fdmath.md

@@ -0,0 +1,13 @@
+# fdmath
+
+::: meanas.fdmath
+
+## Functional and sparse operators
+
+::: meanas.fdmath.functional
+
+::: meanas.fdmath.operators
+
+::: meanas.fdmath.vectorization
+
+::: meanas.fdmath.types

docs/api/fdtd.md

@@ -0,0 +1,15 @@
+# fdtd
+
+::: meanas.fdtd
+
+## Core update and analysis helpers
+
+::: meanas.fdtd.base
+
+::: meanas.fdtd.pml
+
+::: meanas.fdtd.boundaries
+
+::: meanas.fdtd.energy
+
+::: meanas.fdtd.phasor

docs/api/index.md

@@ -0,0 +1,14 @@
+# API Overview
+
+The package is documented directly from its docstrings. The most useful entry
+points are:
+
+- [meanas](meanas.md): top-level package overview
+- [eigensolvers](eigensolvers.md): generic eigenvalue utilities used by the mode solvers
+- [fdfd](fdfd.md): frequency-domain operators, sources, PML, solvers, and far-field transforms
+- [waveguides](waveguides.md): straight, cylindrical, and 3D waveguide mode helpers
+- [fdtd](fdtd.md): timestepping, CPML, energy/flux helpers, and phasor extraction
+- [fdmath](fdmath.md): shared discrete operators, vectorization helpers, and derivation background
+
+The waveguide and FDTD pages are the best places to start if you want the
+mathematical derivations rather than just the callable reference.

docs/api/meanas.md

@@ -0,0 +1,3 @@
+# meanas
+
+::: meanas

docs/api/waveguides.md

@@ -0,0 +1,7 @@
+# waveguides
+
+::: meanas.fdfd.waveguide_2d
+
+::: meanas.fdfd.waveguide_3d
+
+::: meanas.fdfd.waveguide_cyl

docs/assets/vendor/mathjax/manifest.json

@@ -0,0 +1,38 @@
+{
+  "etag": "7ab8013dfc2c395304109fe189c1772ef87b366a",
+  "files": {
+    "core.js": 212262,
+    "loader.js": 16370,
+    "startup.js": 22870,
+    "input/tex-full.js": 270711,
+    "input/asciimath.js": 107690,
+    "input/mml.js": 13601,
+    "input/mml/entities.js": 32232,
+    "output/chtml.js": 211041,
+    "output/chtml/fonts/tex.js": 104359,
+    "output/chtml/fonts/woff-v2/MathJax_Zero.woff": 1368,
+    "output/chtml/fonts/woff-v2/MathJax_Vector-Regular.woff": 1136,
+    "output/chtml/fonts/woff-v2/MathJax_Vector-Bold.woff": 1116,
+    "output/chtml/fonts/woff-v2/MathJax_Typewriter-Regular.woff": 17604,
+    "output/chtml/fonts/woff-v2/MathJax_Size4-Regular.woff": 5148,
+    "output/chtml/fonts/woff-v2/MathJax_Size3-Regular.woff": 3244,
+    "output/chtml/fonts/woff-v2/MathJax_Size2-Regular.woff": 5464,
+    "output/chtml/fonts/woff-v2/MathJax_Size1-Regular.woff": 5792,
+    "output/chtml/fonts/woff-v2/MathJax_Script-Regular.woff": 11852,
+    "output/chtml/fonts/woff-v2/MathJax_SansSerif-Regular.woff": 12660,
+    "output/chtml/fonts/woff-v2/MathJax_SansSerif-Italic.woff": 14628,
+    "output/chtml/fonts/woff-v2/MathJax_SansSerif-Bold.woff": 15944,
+    "output/chtml/fonts/woff-v2/MathJax_Math-Regular.woff": 19288,
+    "output/chtml/fonts/woff-v2/MathJax_Math-Italic.woff": 19360,
+    "output/chtml/fonts/woff-v2/MathJax_Math-BoldItalic.woff": 19776,
+    "output/chtml/fonts/woff-v2/MathJax_Main-Regular.woff": 34160,
+    "output/chtml/fonts/woff-v2/MathJax_Main-Italic.woff": 20832,
+    "output/chtml/fonts/woff-v2/MathJax_Main-Bold.woff": 34464,
+    "output/chtml/fonts/woff-v2/MathJax_Fraktur-Regular.woff": 21480,
+    "output/chtml/fonts/woff-v2/MathJax_Fraktur-Bold.woff": 22340,
+    "output/chtml/fonts/woff-v2/MathJax_Calligraphic-Regular.woff": 9600,
+    "output/chtml/fonts/woff-v2/MathJax_Calligraphic-Bold.woff": 9908,
+    "output/chtml/fonts/woff-v2/MathJax_AMS-Regular.woff": 40808
+  },
+  "version": "3.1.4"
+}

docs/index.md

@@ -0,0 +1,33 @@
+# meanas
+
+`meanas` is a Python package for finite-difference electromagnetic simulation.
+It combines:
+
+- `meanas.fdfd` for frequency-domain operators, sources, waveguide modes, and SCPML
+- `meanas.fdtd` for Yee-grid timestepping, CPML, energy/flux accounting, and phasor extraction
+- `meanas.fdmath` for the shared discrete operators and derivations underneath both solvers
+
+This documentation is built directly from the package docstrings. The API pages
+are the source of truth for the mathematical derivations and calling
+conventions.
+
+## Recommended starting points
+
+- Use the [FDTD API](api/fdtd.md) when you need time-domain stepping, CPML, or
+  phasor extraction.
+- Use the [FDFD API](api/fdfd.md) when you need driven frequency-domain solves
+  or operator algebra.
+- Use the [Waveguide API](api/waveguides.md) for mode solving, port sources, and
+  overlap windows.
+- Use the [fdmath API](api/fdmath.md) when you need the lower-level finite-difference
+  operators or the derivation background shared across the package.
+
+## Build outputs
+
+The docs build generates two HTML views from the same source:
+
+- a normal multi-page site
+- a print-oriented combined page under `site/print_page/`
+
+If `htmlark` is installed, `./make_docs.sh` also writes a fully inlined
+`site/standalone.html`.

docs/javascripts/mathjax.js

@@ -0,0 +1,19 @@
+window.MathJax = {
+  loader: {
+    load: ["input/tex-full", "output/chtml"]
+  },
+  tex: {
+    processEscapes: true,
+    processEnvironments: true
+  },
+  options: {
+    ignoreHtmlClass: ".*|",
+    processHtmlClass: "arithmatex"
+  }
+};
+
+document$.subscribe(() => {
+  MathJax.typesetPromise?.(
+    document.querySelectorAll(".arithmatex")
+  ).catch((error) => console.error(error));
+});

docs/stylesheets/extra.css

@@ -0,0 +1,13 @@
+.md-typeset .arithmatex {
+  overflow-x: auto;
+}
+
+.md-typeset .doc-contents {
+  overflow-wrap: anywhere;
+}
+
+.md-typeset h1 code,
+.md-typeset h2 code,
+.md-typeset h3 code {
+  word-break: break-word;
+}

examples/fdfd1.py

@@ -1,4 +1,5 @@
 import importlib
+import logging
 import numpy
 from numpy.linalg import norm
 from matplotlib import pyplot, colors
@@ -6,12 +7,14 @@ import logging
 
 import meanas
 from meanas import fdtd
-from meanas.fdmath import vec, unvec
+from meanas.fdmath import vec, unvec, fdfield_t
 from meanas.fdfd import waveguide_3d, functional, scpml, operators
 from meanas.fdfd.solvers import generic as generic_solver
 
 import gridlock
 
+from matplotlib import pyplot
+
 
 logging.basicConfig(level=logging.DEBUG)
 logging.getLogger('matplotlib').setLevel(logging.WARNING)

examples/fdtd.py

@@ -1,18 +1,30 @@
 """
-Example code for running an OpenCL FDTD simulation
+Example code for a broadband FDTD run with phasor extraction.
 
-See main() for simulation setup.
+This script shows the intended low-level workflow for:
+
+1. building a Yee-grid simulation with CPML on all faces,
+2. driving it with an electric-current pulse,
+3. extracting a single-frequency phasor on the fly, and
+4. checking that phasor against the matching stretched-grid FDFD operator.
 """
 
 import sys
 import time
+import copy
 
 import numpy
 import h5py
+from numpy.linalg import norm
 
 from meanas import fdtd
 from meanas.fdtd import cpml_params, updates_with_cpml
-from masque import Pattern, shapes
+from meanas.fdtd.misc import gaussian_packet
+
+from meanas.fdfd.operators import e_full
+from meanas.fdfd.scpml import stretch_with_scpml
+from meanas.fdmath import vec
+from masque import Pattern, Circle, Polygon
 import gridlock
 import pcgen
 
@@ -41,50 +53,51 @@ def perturbed_l3(a: float, radius: float, **kwargs) -> Pattern:
         `masque.Pattern` object containing the L3 design
     """
 
-    default_args = {'hole_dose':    1,
+    default_args = {
-                    'trench_dose':  1,
+        'hole_layer':   0,
-                    'hole_layer':   0,
+        'trench_layer': 1,
-                    'trench_layer': 1,
+        'shifts_a':     (0.15, 0, 0.075),
-                    'shifts_a':     (0.15, 0, 0.075),
+        'shifts_r':     (1.0, 1.0, 1.0),
-                    'shifts_r':     (1.0, 1.0, 1.0),
+        'xy_size':      (10, 10),
-                    'xy_size':      (10, 10),
+        'perturbed_radius': 1.1,
-                    'perturbed_radius': 1.1,
+        'trench_width': 1.2e3,
-                    'trench_width': 1.2e3,
+        }
-                    }
     kwargs = {**default_args, **kwargs}
 
-    xyr = pcgen.l3_shift_perturbed_defect(mirror_dims=kwargs['xy_size'],
+    xyr = pcgen.l3_shift_perturbed_defect(
-                                          perturbed_radius=kwargs['perturbed_radius'],
+        mirror_dims=kwargs['xy_size'],
-                                          shifts_a=kwargs['shifts_a'],
+        perturbed_radius=kwargs['perturbed_radius'],
-                                          shifts_r=kwargs['shifts_r'])
+        shifts_a=kwargs['shifts_a'],
+        shifts_r=kwargs['shifts_r'],
+        )
     xyr *= a
     xyr[:, 2] *= radius
 
     pat = Pattern()
-    pat.name = f'L3p-a{a:g}r{radius:g}rp{kwargs["perturbed_radius"]:g}'
+    #pat.name = f'L3p-a{a:g}r{radius:g}rp{kwargs["perturbed_radius"]:g}'
-    pat.shapes += [shapes.Circle(radius=r, offset=(x, y),
+    pat.shapes[(kwargs['hole_layer'], 0)] += [
-                                 dose=kwargs['hole_dose'],
+        Circle(radius=r, offset=(x, y))
-                                 layer=kwargs['hole_layer'])
+        for x, y, r in xyr]
-                   for x, y, r in xyr]
 
     maxes = numpy.max(numpy.fabs(xyr), axis=0)
-    pat.shapes += [shapes.Polygon.rectangle(
+    pat.shapes[(kwargs['trench_layer'], 0)] += [
-        lx=(2 * maxes[0]), ly=kwargs['trench_width'],
+        Polygon.rectangle(
-        offset=(0, s * (maxes[1] + a + kwargs['trench_width'] / 2)),
+            lx=(2 * maxes[0]), ly=kwargs['trench_width'],
-        dose=kwargs['trench_dose'], layer=kwargs['trench_layer'])
+            offset=(0, s * (maxes[1] + a + kwargs['trench_width'] / 2))
-                   for s in (-1, 1)]
+            )
+        for s in (-1, 1)]
     return pat
 
 
 def main():
     dtype = numpy.float32
-    max_t = 8000            # number of timesteps
+    max_t = 3600            # number of timesteps
 
     dx = 40                 # discretization (nm/cell)
     pml_thickness = 8       # (number of cells)
 
     wl = 1550               # Excitation wavelength and fwhm
-    dwl = 200
+    dwl = 100
 
     # Device design parameters
     xy_size = numpy.array([10, 10])
@@ -107,69 +120,106 @@ def main():
 
     # #### Create the grid, mask, and draw the device ####
     grid = gridlock.Grid(edge_coords)
-    epsilon = grid.allocate(n_air**2, dtype=dtype)
+    epsilon = grid.allocate(n_air ** 2, dtype=dtype)
-    grid.draw_slab(epsilon,
+    grid.draw_slab(
-                   surface_normal=2,
+        epsilon,
-                   center=[0, 0, 0],
+        slab = dict(axis='z', center=0, span=th),
-                   thickness=th,
+        foreground = n_slab ** 2,
-                   eps=n_slab**2)
+        )
+
     mask = perturbed_l3(a, r)
+    grid.draw_polygons(
+        epsilon,
+        slab = dict(axis='z', center=0, span=2 * th),
+        foreground = n_air ** 2,
+        offset2d = (0, 0),
+        polygons = mask.as_polygons(library=None),
+        )
 
-    grid.draw_polygons(epsilon,
+    print(f'{grid.shape=}')
-                       surface_normal=2,
-                       center=[0, 0, 0],
-                       thickness=2 * th,
-                       eps=n_air**2,
-                       polygons=mask.as_polygons())
 
-    print(grid.shape)
+    dt = dx * 0.99 / numpy.sqrt(3)
-
+    ee = numpy.zeros_like(epsilon, dtype=dtype)
-    dt = .99/numpy.sqrt(3)
+    hh = numpy.zeros_like(epsilon, dtype=dtype)
-    e = [numpy.zeros_like(epsilon[0], dtype=dtype) for _ in range(3)]
-    h = [numpy.zeros_like(epsilon[0], dtype=dtype) for _ in range(3)]
 
     dxes = [grid.dxyz, grid.autoshifted_dxyz()]
 
     # PMLs in every direction
-    pml_params = [[cpml_params(axis=dd, polarity=pp, dt=dt,
+    pml_params = [
-                               thickness=pml_thickness, epsilon_eff=1.0**2)
+        [cpml_params(axis=dd, polarity=pp, dt=dt, thickness=pml_thickness, epsilon_eff=n_air ** 2)
-                   for pp in (-1, +1)]
+         for pp in (-1, +1)]
-                  for dd in range(3)]
+        for dd in range(3)]
-    update_E, update_H = updates_with_cpml(cpml_params=pml_params, dt=dt,
+    update_E, update_H = updates_with_cpml(cpml_params=pml_params, dt=dt, dxes=dxes, epsilon=epsilon)
-                                           dxes=dxes, epsilon=epsilon)
 
-    # Source parameters and function
+    # sample_interval = numpy.floor(1 / (2 * 1 / wl * dt)).astype(int)
-    w = 2 * numpy.pi * dx / wl
+    # print(f'Save time interval would be {sample_interval} * dt = {sample_interval * dt:3g}')
-    fwhm = dwl * w * w / (2 * numpy.pi * dx)
-    alpha = (fwhm ** 2) / 8 * numpy.log(2)
-    delay = 7/numpy.sqrt(2 * alpha)
 
-    def field_source(i):
+
-        t0 = i * dt - delay
+    # Source parameters and function. The pulse normalization is kept outside
-        return numpy.sin(w * t0) * numpy.exp(-alpha * t0**2)
+    # accumulate_phasor(); the helper only performs the Fourier sum.
+    source_phasor, delay = gaussian_packet(wl=wl, dwl=100, dt=dt, turn_on=1e-5)
+    aa, cc, ss = source_phasor(numpy.arange(max_t))
+    srca_real = aa * cc
+    src_maxt = numpy.argwhere(numpy.diff(aa < 1e-5))[-1]
+    assert aa[src_maxt - 1] >= 1e-5
+    phasor_norm = dt / (aa * cc * cc).sum()
+
+    Jph = numpy.zeros_like(epsilon, dtype=complex)
+    Jph[1, *(grid.shape // 2)] = epsilon[1, *(grid.shape // 2)]
+    omega = 2 * numpy.pi / wl
+    eph = numpy.zeros((1, *epsilon.shape), dtype=complex)
 
     # #### Run a bunch of iterations ####
     output_file = h5py.File('simulation_output.h5', 'w')
     start = time.perf_counter()
-    for t in range(max_t):
+    for tt in range(max_t):
-        update_E(e, h, epsilon)
+        update_E(ee, hh, epsilon)
 
-        e[1][tuple(grid.shape//2)] += field_source(t)
+        if tt < src_maxt:
-        update_H(e, h)
+            # Electric-current injection uses E -= dt * J / epsilon, which is
+            # the same sign convention used by the matching FDFD right-hand side.
+            ee[1, *(grid.shape // 2)] -= srca_real[tt]
+        update_H(ee, hh)
 
-        avg_rate = (t + 1)/(time.perf_counter() - start))
+        avg_rate = (tt + 1) / (time.perf_counter() - start)
-        print(f'iteration {t}: average {avg_rate} iterations per sec')
         sys.stdout.flush()
 
-        if t % 20 == 0:
+        if tt % 200 == 0:
-            r = sum([(f * f * e).sum() for f, e in zip(e, epsilon)])
+            print(f'iteration {tt}: average {avg_rate} iterations per sec')
-            print('E sum', r)
+            E_energy_sum = (ee * ee * epsilon).sum()
+            print(f'{E_energy_sum=}')
 
         # Save field slices
-        if (t % 20 == 0 and (max_t - t <= 1000 or t <= 2000)) or t == max_t-1:
+        if (tt % 20 == 0 and (max_t - tt <= 1000 or tt <= 2000)) or tt == max_t - 1:
-            print('saving E-field')
+            print(f'saving E-field at iteration {tt}')
-            for j, f in enumerate(e):
+            output_file[f'/E_t{tt}'] = ee[:, :, :, ee.shape[3] // 2]
-                output_file['/E{}_t{}'.format('xyz'[j], t)] = f[:, :, round(f.shape[2]/2)]
+
+        fdtd.accumulate_phasor(
+            eph,
+            omega,
+            dt,
+            ee,
+            tt,
+            # The pulse is delayed relative to t=0, so the extracted phasor
+            # needs the same phase offset in its sample times.
+            offset_steps=0.5 - delay / dt,
+            # accumulate_phasor() already multiplies by dt, so pass the
+            # discrete-sum normalization without its extra dt factor.
+            weight=phasor_norm / dt,
+        )
+
+    Eph = eph[0]
+    b = -1j * omega * Jph
+    dxes_fdfd = copy.deepcopy(dxes)
+    for pp in (-1, +1):
+        for dd in range(3):
+            stretch_with_scpml(dxes_fdfd, axis=dd, polarity=pp, omega=omega, epsilon_effective=n_air ** 2, thickness=pml_thickness)
+    # Compare the extracted phasor to the FDFD operator on the stretched grid,
+    # not the unstretched Yee spacings used by the raw time-domain update.
+    A = e_full(omega=omega, dxes=dxes_fdfd, epsilon=epsilon)
+    residual = norm(A @ vec(Eph) - vec(b)) / norm(vec(b))
+    print(f'FDFD residual is {residual}')
+
 
 if __name__ == '__main__':
     main()

examples/waveguide.py

@@ -0,0 +1,351 @@
+"""
+Example code for guided-wave FDFD and FDTD comparison.
+
+This example is the reference workflow for:
+
+1. solving a waveguide port mode,
+2. turning that mode into a one-sided source and overlap window,
+3. comparing a direct FDFD solve against a time-domain phasor extracted from FDTD.
+"""
+from typing import Callable
+import logging
+import time
+import copy
+
+import numpy
+import h5py
+from numpy.linalg import norm
+
+import gridlock
+import meanas
+from meanas import fdtd, fdfd
+from meanas.fdtd import cpml_params, updates_with_cpml
+from meanas.fdtd.misc import gaussian_packet
+
+from meanas.fdmath import vec, unvec, vcfdfield_t, cfdfield_t, fdfield_t, dx_lists_t
+from meanas.fdfd import waveguide_3d, functional, scpml, operators
+from meanas.fdfd.solvers import generic as generic_solver
+from meanas.fdfd.operators import e_full
+from meanas.fdfd.scpml import stretch_with_scpml
+
+
+logging.basicConfig(level=logging.DEBUG)
+for pp in ('matplotlib', 'PIL'):
+    logging.getLogger(pp).setLevel(logging.WARNING)
+
+logger = logging.getLogger(__name__)
+
+
+def pcolor(vv, fig=None, ax=None) -> None:
+    if fig is None:
+        assert ax is None
+        fig, ax = pyplot.subplots()
+    mb = ax.pcolor(vv, cmap='seismic', norm=colors.CenteredNorm())
+    fig.colorbar(mb)
+    ax.set_aspect('equal')
+
+
+def draw_grid(
+        *,
+        dx: float,
+        pml_thickness: int,
+        n_wg: float = 3.476,        # Si index @ 1550
+        n_cladding: float = 1.00,    # Air index
+        wg_w: float = 400,
+        wg_th: float = 200,
+        ) -> tuple[gridlock.Grid, fdfield_t]:
+    """ Create the grid and draw the device  """
+    # Half-dimensions of the simulation grid
+    xyz_max = numpy.array([800, 900, 600]) + (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]
+    edge_coords = [numpy.hstack((-h[::-1], h)) for h in half_edge_coords]
+
+    grid = gridlock.Grid(edge_coords)
+    epsilon = grid.allocate(n_cladding**2, dtype=numpy.float32)
+    grid.draw_cuboid(
+        epsilon,
+        x = dict(center=0, span=8e3),
+        y = dict(center=0, span=wg_w),
+        z = dict(center=0, span=wg_th),
+        foreground = n_wg ** 2,
+        )
+
+    return grid, epsilon
+
+
+def get_waveguide_mode(
+        *,
+        grid: gridlock.Grid,
+        dxes: dx_lists_t,
+        omega: float,
+        epsilon: fdfield_t,
+        ) -> tuple[vcfdfield_t, vcfdfield_t]:
+    """Create a mode source and overlap window for one forward-going port."""
+    dims = numpy.array([[-10, -20, -15],
+                        [-10,  20,  15]]) * [[numpy.median(numpy.real(dx)) for dx in dxes[0]]]
+    ind_dims = (grid.pos2ind(dims[0], which_shifts=None).astype(int),
+                grid.pos2ind(dims[1], which_shifts=None).astype(int))
+    wg_args = dict(
+        slices = [slice(i, f+1) for i, f in zip(*ind_dims)],
+        dxes = dxes,
+        axis = 0,
+        polarity = +1,
+        )
+
+    wg_results = waveguide_3d.solve_mode(mode_number=0, omega=omega, epsilon=epsilon, **wg_args)
+    J = waveguide_3d.compute_source(E=wg_results['E'], wavenumber=wg_results['wavenumber'],
+                                    omega=omega, epsilon=epsilon, **wg_args)
+
+    # compute_overlap_e() returns the normalized upstream overlap window used to
+    # project another field onto this same guided mode.
+    e_overlap = waveguide_3d.compute_overlap_e(E=wg_results['E'], wavenumber=wg_results['wavenumber'], **wg_args, omega=omega)
+    return J, e_overlap
+
+
+def main(
+        *,
+        solver: Callable = generic_solver,
+        dx: float = 40,                  # discretization (nm / cell)
+        pml_thickness: int = 10,         # (number of cells)
+        wl: float = 1550,                # Excitation wavelength
+        wg_w: float = 600,               # Waveguide width
+        wg_th: float = 220,              # Waveguide thickness
+        ):
+    omega = 2 * numpy.pi / wl
+
+    grid, epsilon = draw_grid(dx=dx, pml_thickness=pml_thickness)
+
+    # Add SCPML stretching to the FDFD grid; this matches the CPML-backed FDTD
+    # run below so the two solvers see the same absorbing boundary profile.
+    dxes = [grid.dxyz, grid.autoshifted_dxyz()]
+    for a in (0, 1, 2):
+        for p in (-1, 1):
+            dxes = scpml.stretch_with_scpml(dxes, omega=omega, axis=a, polarity=p, thickness=pml_thickness)
+
+
+    J, e_overlap = get_waveguide_mode(grid=grid, dxes=dxes, omega=omega, epsilon=epsilon)
+
+
+    pecg = numpy.zeros_like(epsilon)
+    # pecg.draw_cuboid(pecg, center=[700, 0, 0], dimensions=[80, 1e8, 1e8], eps=1)
+    # pecg.visualize_isosurface(pecg)
+
+    pmcg = numpy.zeros_like(epsilon)
+    # grid.draw_cuboid(pmcg, center=[700, 0, 0], dimensions=[80, 1e8, 1e8], eps=1)
+    # grid.visualize_isosurface(pmcg)
+
+
+#    ss = (1, slice(None), J.shape[2]//2+6, slice(None))
+#    pcolor(J3[ss].T.imag)
+#    pcolor((numpy.abs(J3).sum(axis=(0, 2)) > 0).astype(float).T)
+#    pyplot.show(block=True)
+
+    E_fdfd = fdfd_solve(
+        omega = omega,
+        dxes = dxes,
+        epsilon = epsilon,
+        J = J,
+        pec = pecg,
+        pmc = pmcg,
+        )
+
+
+    #
+    # Plot results
+    #
+    center = grid.pos2ind([0, 0, 0], None).astype(int)
+    fig, axes = pyplot.subplots(2, 2)
+    pcolor(numpy.real(E[1][center[0], :, :]).T, fig=fig, ax=axes[0, 0])
+    axes[0, 1].plot(numpy.log10(numpy.abs(E[1][:, center[1], center[2]]) + 1e-10))
+    axes[0, 1].grid(alpha=0.6)
+    axes[0, 1].set_ylabel('log10 of field')
+    pcolor(numpy.real(E[1][:, :, center[2]]).T, fig=fig, ax=axes[1, 0])
+
+    def poyntings(E):
+        H = functional.e2h(omega, dxes)(E)
+        poynting = fdtd.poynting(e=E, h=H.conj(), dxes=dxes)
+        cross1 = operators.poynting_e_cross(vec(E), dxes) @ vec(H).conj()
+        cross2 = operators.poynting_h_cross(vec(H), dxes) @ vec(E).conj() * -1
+        s1 = 0.5 * unvec(numpy.real(cross1), grid.shape)
+        s2 = 0.5 * unvec(numpy.real(cross2), grid.shape)
+        s0 = 0.5 * poynting.real
+#        s2 = poynting.imag
+        return s0, s1, s2
+
+    s0x, s1x, s2x = poyntings(E)
+    axes[1, 1].plot(s0x[0].sum(axis=2).sum(axis=1), label='s0', marker='.')
+    axes[1, 1].plot(s1x[0].sum(axis=2).sum(axis=1), label='s1', marker='.')
+    axes[1, 1].plot(s2x[0].sum(axis=2).sum(axis=1), label='s2', marker='.')
+    axes[1, 1].plot(E[1][:, center[1], center[2]].real.T, label='Ey', marker='x')
+    axes[1, 1].grid(alpha=0.6)
+    axes[1, 1].legend()
+
+    q = []
+    for i in range(-5, 30):
+        e_ovl_rolled = numpy.roll(e_overlap, i, axis=1)
+        q += [numpy.abs(vec(E) @ vec(e_ovl_rolled).conj())]
+    fig, ax = pyplot.subplots()
+    ax.plot(q, marker='.')
+    ax.grid(alpha=0.6)
+    ax.set_title('Overlap with mode')
+
+    logger.info('Average overlap with mode:', sum(q[8:32]) / len(q[8:32]))
+
+    pyplot.show(block=True)
+
+
+def fdfd_solve(
+        *,
+        omega: float,
+        dxes = dx_lists_t,
+        epsilon: fdfield_t,
+        J: cfdfield_t,
+        pec: fdfield_t,
+        pmc: fdfield_t,
+        ) -> cfdfield_t:
+    """ Construct and run the solve """
+    sim_args = dict(
+        omega = omega,
+        dxes = dxes,
+        epsilon = vec(epsilon),
+        pec = vec(pecg),
+        pmc = vec(pmcg),
+        )
+
+    x = solver(J=vec(J), **sim_args)
+
+    b = -1j * omega * vec(J)
+    A = operators.e_full(**sim_args).tocsr()
+    logger.info('Norm of the residual is ', norm(A @ x - b) / norm(b))
+
+    E = unvec(x, epsilon.shape[1:])
+    return E
+
+
+def main2():
+    dtype = numpy.float32
+    max_t = 3600            # number of timesteps
+
+    dx = 40                 # discretization (nm/cell)
+    pml_thickness = 8       # (number of cells)
+
+    wl = 1550               # Excitation wavelength and fwhm
+    dwl = 100
+
+    # Device design parameters
+    xy_size = numpy.array([10, 10])
+    a = 430
+    r = 0.285
+    th = 170
+
+    # refractive indices
+    n_slab = 3.408  # InGaAsP(80, 50) @ 1550nm
+    n_cladding = 1.0        # air
+
+    # Half-dimensions of the simulation grid
+    xy_max = (xy_size + 1) * a * [1, numpy.sqrt(3)/2]
+    z_max = 1.6 * a
+    xyz_max = numpy.hstack((xy_max, z_max)) + pml_thickness * dx
+
+    # Coordinates of the edges of the cells. The fdtd package can only do square grids at the moment.
+    half_edge_coords = [numpy.arange(dx/2, m + dx, step=dx) for m in xyz_max]
+    edge_coords = [numpy.hstack((-h[::-1], h)) for h in half_edge_coords]
+
+    # #### Create the grid, mask, and draw the device ####
+    grid = gridlock.Grid(edge_coords)
+    epsilon = grid.allocate(n_cladding ** 2, dtype=dtype)
+    grid.draw_slab(
+        epsilon,
+        slab = dict(axis='z', center=0, span=th),
+        foreground = n_slab ** 2,
+        )
+
+
+    print(f'{grid.shape=}')
+
+    dt = dx * 0.99 / numpy.sqrt(3)
+    ee = numpy.zeros_like(epsilon, dtype=dtype)
+    hh = numpy.zeros_like(epsilon, dtype=dtype)
+
+    dxes = [grid.dxyz, grid.autoshifted_dxyz()]
+
+    # PMLs in every direction
+    pml_params = [
+        [cpml_params(axis=dd, polarity=pp, dt=dt, thickness=pml_thickness, epsilon_eff=n_cladding ** 2)
+         for pp in (-1, +1)]
+        for dd in range(3)]
+    update_E, update_H = updates_with_cpml(cpml_params=pml_params, dt=dt, dxes=dxes, epsilon=epsilon)
+
+    # sample_interval = numpy.floor(1 / (2 * 1 / wl * dt)).astype(int)
+    # print(f'Save time interval would be {sample_interval} * dt = {sample_interval * dt:3g}')
+
+
+    # Source parameters and function. The phasor helper only performs the
+    # Fourier accumulation; the pulse normalization stays explicit here.
+    source_phasor, delay = gaussian_packet(wl=wl, dwl=100, dt=dt, turn_on=1e-5)
+    aa, cc, ss = source_phasor(numpy.arange(max_t))
+    srca_real = aa * cc
+    src_maxt = numpy.argwhere(numpy.diff(aa < 1e-5))[-1]
+    assert aa[src_maxt - 1] >= 1e-5
+    phasor_norm = dt / (aa * cc * cc).sum()
+
+    Jph = numpy.zeros_like(epsilon, dtype=complex)
+    Jph[1, *(grid.shape // 2)] = epsilon[1, *(grid.shape // 2)]
+    omega = 2 * numpy.pi / wl
+    eph = numpy.zeros((1, *epsilon.shape), dtype=complex)
+
+    # #### Run a bunch of iterations ####
+    output_file = h5py.File('simulation_output.h5', 'w')
+    start = time.perf_counter()
+    for tt in range(max_t):
+        update_E(ee, hh, epsilon)
+
+        if tt < src_maxt:
+            # Electric-current injection uses E -= dt * J / epsilon, which is
+            # the sign convention matched by the FDFD source term -1j * omega * J.
+            ee[1, *(grid.shape // 2)] -= srca_real[tt]
+        update_H(ee, hh)
+
+        avg_rate = (tt + 1) / (time.perf_counter() - start)
+
+        if tt % 200 == 0:
+            print(f'iteration {tt}: average {avg_rate} iterations per sec')
+            E_energy_sum = (ee * ee * epsilon).sum()
+            print(f'{E_energy_sum=}')
+
+        # Save field slices
+        if (tt % 20 == 0 and (max_t - tt <= 1000 or tt <= 2000)) or tt == max_t - 1:
+            print(f'saving E-field at iteration {tt}')
+            output_file[f'/E_t{tt}'] = ee[:, :, :, ee.shape[3] // 2]
+
+        fdtd.accumulate_phasor(
+            eph,
+            omega,
+            dt,
+            ee,
+            tt,
+            # The pulse is delayed relative to t=0, so the extracted phasor must
+            # apply the same delay to its sample times.
+            offset_steps=0.5 - delay / dt,
+            # accumulate_phasor() already contributes dt, so remove the extra dt
+            # from the externally computed normalization.
+            weight=phasor_norm / dt,
+        )
+
+    Eph = eph[0]
+    b = -1j * omega * Jph
+    dxes_fdfd = copy.deepcopy(dxes)
+    for pp in (-1, +1):
+        for dd in range(3):
+            stretch_with_scpml(dxes_fdfd, axis=dd, polarity=pp, omega=omega, epsilon_effective=n_cladding ** 2, thickness=pml_thickness)
+    # Residuals must be checked on the stretched FDFD grid, because the FDTD run
+    # already includes those same absorbing layers through CPML.
+    A = e_full(omega=omega, dxes=dxes_fdfd, epsilon=epsilon)
+    residual = norm(A @ vec(Eph) - vec(b)) / norm(vec(b))
+    print(f'FDFD residual is {residual}')
+
+
+if __name__ == '__main__':
+    main()

make_docs.sh

@@ -2,18 +2,12 @@
 
 set -Eeuo pipefail
 
-cd ~/projects/meanas
+ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
+cd "$ROOT"
 
-# Approach 1: pdf to html?
+mkdocs build --clean
-#pdoc3 --pdf --force --template-dir pdoc_templates -o doc . | \
-#    pandoc --metadata=title:"meanas" --toc --toc-depth=4 --from=markdown+abbreviations --to=html --output=doc.html --gladtex -s -
 
-# Approach 2: pdf to html with gladtex
+PRINT_PAGE='site/print_page/index.html'
-rm -rf _doc_mathimg
+if [[ -f "$PRINT_PAGE" ]] && command -v htmlark >/dev/null 2>&1; then
-pdoc --pdf --force --template-dir pdoc_templates -o doc . > doc.md
+    htmlark "$PRINT_PAGE" -o site/standalone.html
-pandoc --metadata=title:"meanas" --from=markdown+abbreviations --to=html --output=doc.htex --gladtex -s --css pdoc_templates/pdoc.css doc.md
+fi
-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 .
-#find doc -iname '*.html' -exec gladtex -a -n -d _mathimg -c white {} \;

meanas/fdfd/__init__.py

@@ -9,9 +9,12 @@ Submodules:
 
 - `operators`, `functional`: General FDFD problem setup.
 - `solvers`: Solver interface and reference implementation.
-- `scpml`: Stretched-coordinate perfectly matched layer (scpml) boundary conditions
+- `scpml`: Stretched-coordinate perfectly matched layer (SCPML) boundary conditions.
 - `waveguide_2d`: Operators and mode-solver for waveguides with constant cross-section.
-- `waveguide_3d`: Functions for transforming `waveguide_2d` results into 3D.
+- `waveguide_3d`: Functions for transforming `waveguide_2d` results into 3D,
+  including mode-source and overlap-window construction.
+- `farfield`, `bloch`, `eme`: specialized helper modules for near/far transforms,
+  Bloch-periodic problems, and eigenmode expansion.
 
 
 ================================================================
@@ -86,10 +89,6 @@ $$
     -\omega^2 \epsilon_{\vec{r}} \cdot \tilde{E}_{\vec{r}} = -\imath \omega \tilde{J}_{\vec{r}} \\
 $$
 
-# TODO FDFD?
-# TODO PML
-
-
 """
 from . import (
     solvers as solvers,

meanas/fdfd/bloch.py

@@ -136,6 +136,14 @@ except ImportError:
     logger.info('Using numpy fft')
 
 
+def _assemble_hmn_vector(
+        h_m: NDArray[numpy.complex128],
+        h_n: NDArray[numpy.complex128],
+        ) -> NDArray[numpy.complex128]:
+    stacked = numpy.concatenate((numpy.ravel(h_m), numpy.ravel(h_n)))
+    return stacked[:, None]
+
+
 def generate_kmn(
         k0: ArrayLike,
         G_matrix: ArrayLike,
@@ -253,8 +261,8 @@ def maxwell_operator(
             h_m, h_n = b_m, b_n
         else:
             # transform from mn to xyz
-            b_xyz = (m * b_m[:, :, :, None]
+            b_xyz = (m * b_m
-                   + n * b_n[:, :, :, None])
+                   + n * b_n)    # noqa: E128
 
             # divide by mu
             temp = ifftn(b_xyz, axes=range(3))
@@ -265,10 +273,7 @@ def maxwell_operator(
             h_m = numpy.sum(h_xyz * m, axis=3)
             h_n = numpy.sum(h_xyz * n, axis=3)
 
-        h.shape = (h.size,)
+        return _assemble_hmn_vector(h_m, h_n)
-        h = numpy.concatenate((h_m.ravel(), h_n.ravel()), axis=None, out=h)     # ravel and merge
-        h.shape = (h.size, 1)
-        return h
 
     return operator
 
@@ -403,8 +408,8 @@ def inverse_maxwell_operator_approx(
             b_m, b_n = hin_m, hin_n
         else:
             # transform from mn to xyz
-            h_xyz = (m * hin_m[:, :, :, None]
+            h_xyz = (m * hin_m
-                   + n * hin_n[:, :, :, None])
+                   + n * hin_n)  # noqa: E128
 
             # multiply by mu
             temp = ifftn(h_xyz, axes=range(3))
@@ -412,8 +417,8 @@ def inverse_maxwell_operator_approx(
             b_xyz = fftn(temp, axes=range(3))
 
             # transform back to mn
-            b_m = numpy.sum(b_xyz * m, axis=3)
+            b_m = numpy.sum(b_xyz * m, axis=3, keepdims=True)
-            b_n = numpy.sum(b_xyz * n, axis=3)
+            b_n = numpy.sum(b_xyz * n, axis=3, keepdims=True)
 
         # cross product and transform into xyz basis
         e_xyz = (n * b_m
@@ -428,10 +433,7 @@ def inverse_maxwell_operator_approx(
         h_m = numpy.sum(d_xyz * n, axis=3, keepdims=True) / +k_mag
         h_n = numpy.sum(d_xyz * m, axis=3, keepdims=True) / -k_mag
 
-        h.shape = (h.size,)
+        return _assemble_hmn_vector(h_m, h_n)
-        h = numpy.concatenate((h_m, h_n), axis=None, out=h)
-        h.shape = (h.size, 1)
-        return h
 
     return operator
 
@@ -449,7 +451,7 @@ def find_k(
         solve_callback: Callable[..., None] | None = None,
         iter_callback: Callable[..., None] | None = None,
         v0: NDArray[numpy.complex128] | None = None,
-        ) -> tuple[float, float, NDArray[numpy.complex128], NDArray[numpy.complex128]]:
+        ) -> tuple[NDArray[numpy.float64], float, NDArray[numpy.complex128], NDArray[numpy.complex128]]:
     """
     Search for a bloch vector that has a given frequency.
 
@@ -494,15 +496,15 @@ def find_k(
 
     res = scipy.optimize.minimize_scalar(
         lambda x: abs(get_f(x, band) - frequency),
-        k_guess,
+        method='bounded',
-        method='Bounded',
         bounds=k_bounds,
         options={'xatol': abs(tolerance)},
         )
 
     assert n is not None
     assert v is not None
-    return float(res.x * direction), float(res.fun + frequency), n, v
+    actual_frequency = get_f(float(res.x), band)
+    return direction * float(res.x), float(actual_frequency), n, v
 
 
 def eigsolve(
@@ -723,7 +725,12 @@ def eigsolve(
                 amax=pi,
                 )
 
-        result = scipy.optimize.minimize_scalar(trace_func, bounds=(0, pi), tol=tolerance)
+        result = scipy.optimize.minimize_scalar(
+            trace_func,
+            method='bounded',
+            bounds=(0, pi),
+            options={'xatol': tolerance},
+            )
         new_E = result.fun
         theta = result.x
 
@@ -752,7 +759,7 @@ def eigsolve(
         v = eigvecs[:, i]
         n = eigvals[i]
         v /= norm(v)
-        Av = (scipy_op @ v.copy())[:, 0]
+        Av = numpy.asarray(scipy_op @ v.copy()).reshape(-1)
         eigness = norm(Av - (v.conj() @ Av) * v)
         f = numpy.sqrt(-numpy.real(n))
         df = numpy.sqrt(-numpy.real(n) + eigness)
@@ -821,18 +828,18 @@ def inner_product(
     # eRxhR_x = numpy.cross(eR.reshape(3, -1), hR.reshape(3, -1), axis=0).reshape(eR.shape)[0] / normR
     # logger.info(f'power {eRxhR_x.sum() / 2})
 
-    eR /= numpy.sqrt(norm2R)
+    eR_norm = eR / numpy.sqrt(abs(norm2R))
-    hR /= numpy.sqrt(norm2R)
+    hR_norm = hR / numpy.sqrt(abs(norm2R))
-    eL /= numpy.sqrt(norm2L)
+    eL_norm = eL / numpy.sqrt(abs(norm2L))
-    hL /= numpy.sqrt(norm2L)
+    hL_norm = hL / numpy.sqrt(abs(norm2L))
 
     # (eR x hL)[0] and (eL x hR)[0]
-    eRxhL_x = eR[1] * hL[2] - eR[2] - hL[1]
+    eRxhL_x = eR_norm[1] * hL_norm[2] - eR_norm[2] * hL_norm[1]
-    eLxhR_x = eL[1] * hR[2] - eL[2] - hR[1]
+    eLxhR_x = eL_norm[1] * hR_norm[2] - eL_norm[2] * hR_norm[1]
 
     #return 1j *  (eRxhL_x - eLxhR_x).sum() / numpy.sqrt(norm2R * norm2L)
     #return (eRxhL_x.sum() - eLxhR_x.sum()) / numpy.sqrt(norm2R * norm2L)
-    return eRxhL_x.sum() - eLxhR_x.sum()
+    return eLxhR_x.sum() - eRxhL_x.sum()
 
 
 def trq(
@@ -846,8 +853,8 @@ def trq(
     np = inner_product(eO, -hO, eI,  hI)
     nn = inner_product(eO, -hO, eI, -hI)
 
-    assert pp == -nn
+    assert numpy.allclose(pp, -nn, atol=1e-12, rtol=1e-12)
-    assert pn == -np
+    assert numpy.allclose(pn, -np, atol=1e-12, rtol=1e-12)
 
     logger.info(f'''
         {pp=:4g} {pn=:4g}

meanas/fdfd/eme.py

@@ -55,7 +55,13 @@ def get_abcd(
     B = r21 @ t21i
     C = -t21i @ r12
     D = t21i
-    return sparse.block_array(((A, B), (C, D)))
+    return sparse.block_array(
+        [
+            [sparse.csr_array(A), sparse.csr_array(B)],
+            [sparse.csr_array(C), sparse.csr_array(D)],
+            ],
+        format='csr',
+        )
 
 
 def get_s(
@@ -75,8 +81,8 @@ def get_s(
 
     if force_nogain:
         # force S @ S.H diagonal
-        U, sing, V = numpy.linalg.svd(ss)
+        U, sing, Vh = numpy.linalg.svd(ss)
-        ss = numpy.diag(sing) @ U @ V
+        ss = U @ numpy.diag(numpy.minimum(sing, 1.0)) @ Vh
 
     if force_reciprocal:
         ss = 0.5 * (ss + ss.T)

meanas/fdfd/farfield.py

@@ -78,15 +78,12 @@ def near_to_farfield(
     kx, ky = numpy.meshgrid(kxx, kyy, indexing='ij')
     kxy2 = kx * kx + ky * ky
     kxy = numpy.sqrt(kxy2)
-    kz = numpy.sqrt(k * k - kxy2)
+    kz = numpy.sqrt(numpy.maximum(0, k * k - kxy2))
 
-    sin_th = ky / kxy
+    sin_th = numpy.divide(ky, kxy, out=numpy.zeros_like(ky), where=kxy != 0)
-    cos_th = kx / kxy
+    cos_th = numpy.divide(kx, kxy, out=numpy.ones_like(kx), where=kxy != 0)
     cos_phi = kz / k
 
-    sin_th[numpy.logical_and(kx == 0, ky == 0)] = 0
-    cos_th[numpy.logical_and(kx == 0, ky == 0)] = 1
-
     # Normalized vector potentials N, L
     N = [-Hn_fft[1] * cos_phi * cos_th + Hn_fft[0] * cos_phi * sin_th,
           Hn_fft[1] * sin_th + Hn_fft[0] * cos_th]                      # noqa: E127
@@ -114,8 +111,8 @@ def near_to_farfield(
     outputs = {
         'E': E_far,
         'H': H_far,
-        'dkx': kx[1] - kx[0],
+        'dkx': float(kxx[1] - kxx[0]),
-        'dky': ky[1] - ky[0],
+        'dky': float(kyy[1] - kyy[0]),
         'kx': kx,
         'ky': ky,
         'theta': theta,
@@ -177,22 +174,19 @@ def far_to_nearfield(
     padded_shape = cast('Sequence[int]', padded_size)
 
     k = 2 * pi
-    kxs = fftshift(fftfreq(s[0], 1 / (s[0] * dkx)))
+    kxs = dkx * fftshift(fftfreq(s[0], d=1 / s[0]))
-    kys = fftshift(fftfreq(s[0], 1 / (s[1] * dky)))
+    kys = dky * fftshift(fftfreq(s[1], d=1 / s[1]))
 
     kx, ky = numpy.meshgrid(kxs, kys, indexing='ij')
     kxy2 = kx * kx + ky * ky
     kxy = numpy.sqrt(kxy2)
 
-    kz = numpy.sqrt(k * k - kxy2)
+    kz = numpy.sqrt(numpy.maximum(0, k * k - kxy2))
 
-    sin_th = ky / kxy
+    sin_th = numpy.divide(ky, kxy, out=numpy.zeros_like(ky), where=kxy != 0)
-    cos_th = kx / kxy
+    cos_th = numpy.divide(kx, kxy, out=numpy.ones_like(kx), where=kxy != 0)
     cos_phi = kz / k
 
-    sin_th[numpy.logical_and(kx == 0, ky == 0)] = 0
-    cos_th[numpy.logical_and(kx == 0, ky == 0)] = 1
-
     theta = numpy.arctan2(ky, kx)
     phi = numpy.arccos(cos_phi)
     theta[numpy.logical_and(kx == 0, ky == 0)] = 0
@@ -212,21 +206,41 @@ def far_to_nearfield(
     N = [L[1],
         -L[0]]  # noqa: E128
 
-    En_fft = [-( L[0] * sin_th + L[1] * cos_phi * cos_th) / cos_phi,
+    En_fft = [
-              -(-L[0] * cos_th + L[1] * cos_phi * sin_th) / cos_phi]
+        numpy.divide(
+            -(L[0] * sin_th + L[1] * cos_phi * cos_th),
+            cos_phi,
+            out=numpy.zeros_like(L[0]),
+            where=cos_phi != 0,
+            ),
+        numpy.divide(
+            -(-L[0] * cos_th + L[1] * cos_phi * sin_th),
+            cos_phi,
+            out=numpy.zeros_like(L[0]),
+            where=cos_phi != 0,
+            ),
+        ]
 
-    Hn_fft = [( N[0] * sin_th + N[1] * cos_phi * cos_th) / cos_phi,
+    Hn_fft = [
-              (-N[0] * cos_th + N[1] * cos_phi * sin_th) / cos_phi]
+        numpy.divide(
-
+            N[0] * sin_th + N[1] * cos_phi * cos_th,
-    for i in range(2):
+            cos_phi,
-        En_fft[i][cos_phi == 0] = 0
+            out=numpy.zeros_like(N[0]),
-        Hn_fft[i][cos_phi == 0] = 0
+            where=cos_phi != 0,
+            ),
+        numpy.divide(
+            -N[0] * cos_th + N[1] * cos_phi * sin_th,
+            cos_phi,
+            out=numpy.zeros_like(N[0]),
+            where=cos_phi != 0,
+            ),
+        ]
 
     E_near = [ifftshift(ifft2(ifftshift(Ei), s=padded_shape)) for Ei in En_fft]
     H_near = [ifftshift(ifft2(ifftshift(Hi), s=padded_shape)) for Hi in Hn_fft]
 
     dx = 2 * pi / (s[0] * dkx)
-    dy = 2 * pi / (s[0] * dky)
+    dy = 2 * pi / (s[1] * dky)
 
     outputs = {
         'E': E_near,
@@ -236,4 +250,3 @@ def far_to_nearfield(
     }
 
     return outputs
-

meanas/fdfd/functional.py

@@ -41,8 +41,8 @@ def e_full(
         curls = ch(ce(e))
         return cfdfield_t(curls - omega ** 2 * epsilon * e)
 
-    def op_mu(e: cfdfield) -> cfdfield_t:
+    def op_mu(e: cfdfield_t) -> cfdfield_t:
-        curls = ch(mu * ce(e))          # type: ignore   # mu = None ok because we don't return the function
+        curls = ch(ce(e) / mu)          # type: ignore   # mu = None ok because we don't return the function
         return cfdfield_t(curls - omega ** 2 * epsilon * e)
 
     if mu is None:
@@ -138,12 +138,12 @@ def m2j(
     """
     ch = curl_back(dxes[1])
 
-    def m2j_mu(m: cfdfield) -> cfdfield_t:
+    def m2j_mu(m: cfdfield_t) -> cfdfield_t:
-        J = ch(m / mu) / (-1j * omega)          # type: ignore  # mu=None ok
+        J = ch(m / mu) / (1j * omega)           # type: ignore  # mu=None ok
         return cfdfield_t(J)
 
-    def m2j_1(m: cfdfield) -> cfdfield_t:
+    def m2j_1(m: cfdfield_t) -> cfdfield_t:
-        J = ch(m) / (-1j * omega)
+        J = ch(m) / (1j * omega)
         return cfdfield_t(J)
 
     if mu is None:
@@ -158,10 +158,23 @@ def e_tfsf_source(
         epsilon: fdfield,
         mu: fdfield | None = None,
         ) -> cfdfield_updater_t:
-    """
+    r"""
     Operator that turns an E-field distribution into a total-field/scattered-field
     (TFSF) source.
 
+    If `A` is the full wave operator from `e_full(...)` and `Q` is the diagonal
+    mask selecting the total-field region, then the TFSF source is the commutator
+
+    $$
+    \frac{A Q - Q A}{-i \omega} E.
+    $$
+
+    This vanishes in the interior of the total-field and scattered-field regions
+    and is supported only at their shared boundary, where the mask discontinuity
+    makes `A` and `Q` fail to commute. The returned current is therefore the
+    distributed source needed to inject the desired total field without also
+    forcing the scattered-field region.
+
     Args:
         TF_region: mask which is set to 1 in the total-field region, and 0 elsewhere
                    (i.e. in the scattered-field region).
@@ -175,7 +188,6 @@ def e_tfsf_source(
         Function `f` which takes an E field and returns a current distribution,
         `f(E)` -> `J`
     """
-    # TODO documentation
     A = e_full(omega, dxes, epsilon, mu)
 
     def op(e: cfdfield) -> cfdfield_t:
@@ -188,7 +200,13 @@ def poynting_e_cross_h(dxes: dx_lists_t) -> Callable[[cfdfield, cfdfield], cfdfi
     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$
+    for the Poynting vector, $S = E \times H$.
+
+    On the Yee grid, the electric and magnetic components are not stored at the
+    same locations. This helper therefore applies the same one-cell electric-field
+    shifts used by the sparse `operators.poynting_e_cross(...)` construction so
+    that the discrete cross product matches the face-centered energy flux used in
+    `meanas.fdtd.energy.poynting(...)`.
 
     Note:
         This function also shifts the input `E` field by one cell as required
@@ -204,7 +222,8 @@ def poynting_e_cross_h(dxes: dx_lists_t) -> Callable[[cfdfield, cfdfield], cfdfi
         dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types`
 
     Returns:
-        Function `f` that returns E x H as required for the poynting vector.
+        Function `f` that returns the staggered-grid cross product `E \times H`.
+        For time-average power, call it as `f(E, H.conj())` and take `Re(...) / 2`.
     """
     def exh(e: cfdfield, h: cfdfield) -> cfdfield_t:
         s = numpy.empty_like(e)

meanas/fdfd/operators.py

@@ -236,10 +236,12 @@ def eh_full(
     else:
         pm = sparse.diags_array(numpy.where(pmc, 0, 1))    # set pm to (not PMC)
 
-    iwe = pe @ (1j * omega * sparse.diags_array(epsilon)) @ pe
+    iwe = pe @ (1j * omega * sparse.diags(epsilon)) @ pe
-    iwm = 1j * omega
+    if mu is None:
-    if mu is not None:
+        iwm = 1j * omega * sparse.eye(epsilon.size)
-        iwm *= sparse.diags_array(mu)
+    else:
+        iwm = 1j * omega * sparse.diags(mu)
+
     iwm = pm @ iwm @ pm
 
     A1 = pe @ curl_back(dxes[1]) @ pm
@@ -308,16 +310,22 @@ def m2j(
 
 
 def poynting_e_cross(e: vcfdfield, dxes: dx_lists_t) -> sparse.sparray:
-    """
+    r"""
-    Operator for computing the Poynting vector, containing the
+    Operator for computing the staggered-grid `(E \times)` part of the Poynting vector.
-    (E x) portion of the Poynting vector.
+
+    On the Yee grid the E and H components live on different edges, so the
+    electric field must be shifted by one cell in the appropriate direction
+    before forming the discrete cross product. This sparse operator encodes that
+    shifted cross product directly and is the matrix equivalent of
+    `functional.poynting_e_cross_h(...)`.
 
     Args:
         e: Vectorized E-field for the ExH cross product
         dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types`
 
     Returns:
-        Sparse matrix containing (E x) portion of Poynting cross product.
+        Sparse matrix containing the `(E \times)` part of the staggered Poynting
+        cross product.
     """
     shape = [len(dx) for dx in dxes[0]]
 
@@ -337,15 +345,26 @@ def poynting_e_cross(e: vcfdfield, dxes: dx_lists_t) -> sparse.sparray:
 
 
 def poynting_h_cross(h: vcfdfield, dxes: dx_lists_t) -> sparse.sparray:
-    """
+    r"""
-    Operator for computing the Poynting vector, containing the (H x) portion of the Poynting vector.
+    Operator for computing the staggered-grid `(H \times)` part of the Poynting vector.
+
+    Together with `poynting_e_cross(...)`, this provides the matrix form of the
+    Yee-grid cross product used in the flux helpers. The two are related by the
+    usual antisymmetry of the cross product,
+
+    $$
+    H \times E = -(E \times H),
+    $$
+
+    once the same staggered field placement is used on both sides.
 
     Args:
         h: Vectorized H-field for the HxE cross product
         dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types`
 
     Returns:
-        Sparse matrix containing (H x) portion of Poynting cross product.
+        Sparse matrix containing the `(H \times)` part of the staggered Poynting
+        cross product.
     """
     shape = [len(dx) for dx in dxes[0]]
 
@@ -370,11 +389,23 @@ def e_tfsf_source(
         epsilon: vfdfield,
         mu: vfdfield | None = None,
         ) -> sparse.sparray:
-    """
+    r"""
     Operator that turns a desired E-field distribution into a
      total-field/scattered-field (TFSF) source.
 
-    TODO: Reference Rumpf paper
+    Let `A` be the full wave operator from `e_full(...)`, and let
+    `Q = \mathrm{diag}(TF_region)` be the projector onto the total-field region.
+    Then the TFSF current operator is the commutator
+
+    $$
+    \frac{A Q - Q A}{-i \omega}.
+    $$
+
+    Inside regions where `Q` is locally constant, `A` and `Q` commute and the
+    source vanishes. Only cells at the TF/SF boundary contribute nonzero current,
+    which is exactly the desired distributed source for injecting the chosen
+    field into the total-field region without directly forcing the
+    scattered-field region.
 
     Args:
         TF_region: Mask, which is set to 1 inside the total-field region and 0 in the
@@ -386,9 +417,7 @@ def e_tfsf_source(
 
     Returns:
         Sparse matrix that turns an E-field into a current (J) distribution.
-
     """
-    # TODO documentation
     A = e_full(omega, dxes, epsilon, mu)
     Q = sparse.diags_array(TF_region)
     return (A @ Q - Q @ A) / (-1j * omega)
@@ -402,11 +431,17 @@ def e_boundary_source(
         mu: vfdfield | None = None,
         periodic_mask_edges: bool = False,
         ) -> sparse.sparray:
-    """
+    r"""
     Operator that turns an E-field distrubtion into a current (J) distribution
       along the edges (external and internal) of the provided mask. This is just an
       `e_tfsf_source()` with an additional masking step.
 
+    Equivalently, this helper first constructs the TFSF commutator source for the
+    full mask and then zeroes out all cells except the mask boundary. The
+    boundary is defined as the set of cells whose mask value changes under a
+    one-cell shift in any Cartesian direction. With `periodic_mask_edges=False`
+    the shifts mirror at the domain boundary; with `True` they wrap periodically.
+
     Args:
         mask: The current distribution is generated at the edges of the mask,
               i.e. any points where shifting the mask by one cell in any direction

meanas/fdfd/scpml.py

@@ -128,6 +128,11 @@ def stretch_with_scpml(
     dx_ai = dxes[0][axis].astype(complex)
     dx_bi = dxes[1][axis].astype(complex)
 
+    if thickness == 0:
+        dxes[0][axis] = dx_ai
+        dxes[1][axis] = dx_bi
+        return dxes
+
     pos = numpy.hstack((0, dx_ai.cumsum()))
     pos_a = (pos[:-1] + pos[1:]) / 2
     pos_b = pos[:-1]
@@ -153,10 +158,7 @@ def stretch_with_scpml(
         def l_d(x: NDArray[numpy.float64]) -> NDArray[numpy.float64]:
             return (x - bound) / (pos[-1] - bound)
 
-        if thickness == 0:
+        slc = slice(-thickness, None)
-            slc = slice(None)
-        else:
-            slc = slice(-thickness, None)
 
     dx_ai[slc] *= 1 + 1j * s_function(l_d(pos_a[slc])) / d / s_correction
     dx_bi[slc] *= 1 + 1j * s_function(l_d(pos_b[slc])) / d / s_correction

meanas/fdfd/solvers.py

@@ -48,9 +48,11 @@ def _scipy_qmr(
             logger.info(f'Solver residual at iteration {ii} : {cur_norm}')
 
     if 'callback' in kwargs:
+        callback = kwargs['callback']
+
         def augmented_callback(xk: ArrayLike) -> None:
             log_residual(xk)
-            kwargs['callback'](xk)
+            callback(xk)
 
         kwargs['callback'] = augmented_callback
     else:
@@ -118,15 +120,15 @@ def generic(
     Pl, Pr = operators.e_full_preconditioners(dxes)
 
     if adjoint:
-        A = (Pl @ A0 @ Pr).H
+        A = (Pl @ A0 @ Pr).T.conjugate()
-        b = Pr.H @ b0
+        b = Pr.T.conjugate() @ b0
     else:
         A = Pl @ A0 @ Pr
         b = Pl @ b0
 
     if E_guess is not None:
         if adjoint:
-            x0 = Pr.H @ E_guess
+            x0 = Pr.T.conjugate() @ E_guess
         else:
             x0 = Pl @ E_guess
         matrix_solver_opts['x0'] = x0
@@ -134,7 +136,7 @@ def generic(
     x = matrix_solver(A.tocsr(), b, **matrix_solver_opts)
 
     if adjoint:
-        x0 = Pl.H @ x
+        x0 = Pl.T.conjugate() @ x
     else:
         x0 = Pr @ x
 

meanas/fdfd/waveguide_2d.py

@@ -175,8 +175,6 @@ if the result is introduced into a space with a discretized z-axis.
 
 
 """
-# TODO update module docs
-
 from typing import Any
 from collections.abc import Sequence
 import numpy
@@ -339,7 +337,7 @@ def normalized_fields_e(
         mu: vfdslice | None = None,
         prop_phase: float = 0,
         ) -> tuple[vcfdslice_t, vcfdslice_t]:
-    """
+    r"""
     Given a vector `e_xy` containing the vectorized E_x and E_y fields,
      returns normalized, vectorized E and H fields for the system.
 
@@ -357,6 +355,21 @@ def normalized_fields_e(
     Returns:
         `(e, h)`, where each field is vectorized, normalized,
         and contains all three vector components.
+
+    Notes:
+        `e_xy` is only the transverse electric eigenvector. This helper first
+        reconstructs the full three-component `E` and `H` fields with `exy2e(...)`
+        and `exy2h(...)`, then normalizes them to unit forward power using
+        `_normalized_fields(...)`.
+
+        The normalization target is
+
+        $$
+        \Re\left[\mathrm{inner\_product}(e, h, \mathrm{conj\_h}=True)\right] = 1,
+        $$
+
+        so the returned fields represent a unit-power forward mode under the
+        discrete Yee-grid Poynting inner product.
     """
     e = exy2e(wavenumber=wavenumber, dxes=dxes, epsilon=epsilon) @ e_xy
     h = exy2h(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) @ e_xy
@@ -374,7 +387,7 @@ def normalized_fields_h(
         mu: vfdslice | None = None,
         prop_phase: float = 0,
         ) -> tuple[vcfdslice_t, vcfdslice_t]:
-    """
+    r"""
     Given a vector `h_xy` containing the vectorized H_x and H_y fields,
      returns normalized, vectorized E and H fields for the system.
 
@@ -392,6 +405,13 @@ def normalized_fields_h(
     Returns:
         `(e, h)`, where each field is vectorized, normalized,
         and contains all three vector components.
+
+    Notes:
+        This is the `H_x/H_y` analogue of `normalized_fields_e(...)`. The final
+        normalized mode should describe the same physical solution, but because
+        the overall complex phase and sign are chosen heuristically,
+        `normalized_fields_e(...)` and `normalized_fields_h(...)` need not return
+        identical representatives for nearly symmetric modes.
     """
     e = hxy2e(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) @ h_xy
     h = hxy2h(wavenumber=wavenumber, dxes=dxes, mu=mu) @ h_xy
@@ -409,7 +429,25 @@ def _normalized_fields(
         mu: vfdslice | None = None,
         prop_phase: float = 0,
         ) -> tuple[vcfdslice_t, vcfdslice_t]:
-    # TODO documentation
+    r"""
+    Normalize a reconstructed waveguide mode to unit forward power.
+
+    The eigenproblem solved by `solve_mode(s)` determines only the mode shape and
+    propagation constant. The overall complex amplitude and sign are still free.
+    This helper fixes those remaining degrees of freedom in two steps:
+
+    1. Compute the discrete longitudinal Poynting flux with
+       `inner_product(e, h, conj_h=True)`, including the half-cell longitudinal
+       phase adjustment controlled by `prop_phase`.
+    2. Multiply both fields by a scalar chosen so that the real forward power is
+       `1`, then choose a reproducible phase/sign representative by making a
+       dominant-energy sample real and using a weighted quadrant sum to break
+       mirror-symmetry ties.
+
+    The sign heuristic is intentionally pragmatic rather than fundamental: it is
+    only there to make downstream tests and source/overlap construction choose a
+    consistent representative when the physical mode is symmetric.
+    """
     shape = [s.size for s in dxes[0]]
 
     # Find time-averaged Sz and normalize to it
@@ -921,7 +959,7 @@ def solve_mode(
     return vcfdfield2_t(e_xys[0]), wavenumbers[0]
 
 
-def inner_product(    # TODO documentation
+def inner_product(
         e1: vcfdfield2,
         h2: vcfdfield2,
         dxes: dx_lists2_t,
@@ -929,6 +967,36 @@ def inner_product(    # TODO documentation
         conj_h: bool = False,
         trapezoid: bool = False,
         ) -> complex:
+    r"""
+    Compute the discrete waveguide overlap / Poynting inner product.
+
+    This is the 2D transverse integral corresponding to the time-averaged
+    longitudinal Poynting flux,
+
+    $$
+    \frac{1}{2}\int (E_x H_y - E_y H_x) \, dx \, dy
+    $$
+
+    with the Yee-grid staggering and optional propagation-phase adjustment used
+    by the waveguide helpers in this module.
+
+    Args:
+        e1: Vectorized electric field, typically from `exy2e(...)` or
+            `normalized_fields_e(...)`.
+        h2: Vectorized magnetic field, typically from `hxy2h(...)`,
+            `exy2h(...)`, or one of the normalization helpers.
+        dxes: Two-dimensional Yee-grid spacing lists `[dx_e, dx_h]`.
+        prop_phase: Phase advance over one propagation cell. This is used to
+            shift the H field into the same longitudinal reference plane as the
+            E field.
+        conj_h: Whether to conjugate `h2` before forming the overlap. Use
+            `True` for the usual time-averaged power normalization.
+        trapezoid: Whether to use trapezoidal quadrature instead of the default
+            rectangular Yee-cell sum.
+
+    Returns:
+        Complex overlap / longitudinal power integral.
+    """
 
     shape = [s.size for s in dxes[0]]
 
@@ -951,5 +1019,3 @@ def inner_product(    # TODO documentation
         Sz_b = E1[1] * H2[0] * dxes_real[0][0][:, None] * dxes_real[1][1][None, :]
     Sz = 0.5 * (Sz_a.sum() - Sz_b.sum())
     return Sz
-
-

meanas/fdfd/waveguide_3d.py

@@ -3,8 +3,25 @@ Tools for working with waveguide modes in 3D domains.
 
 This module relies heavily on `waveguide_2d` and mostly just transforms
 its parameters into 2D equivalents and expands the results back into 3D.
+
+The intended workflow is:
+
+1. Select a single-cell slice normal to the propagation axis.
+2. Solve the corresponding 2D mode problem with `solve_mode(...)`.
+3. Turn that mode into a one-sided source with `compute_source(...)`.
+4. Build an overlap window with `compute_overlap_e(...)` for port readout.
+
+`polarity` is part of the public convention throughout this module:
+
+- `+1` means forward propagation toward increasing index along `axis`
+- `-1` means backward propagation toward decreasing index along `axis`
+
+That same convention controls which side of the selected slice is used for the
+overlap window and how the expanded field is phased.
 """
 from typing import Any, cast
+import warnings
+from typing import Any
 from collections.abc import Sequence
 import numpy
 from numpy.typing import NDArray
@@ -24,9 +41,9 @@ def solve_mode(
         epsilon: fdfield,
         mu: fdfield | None = None,
         ) -> dict[str, complex | NDArray[complexfloating]]:
-    """
+    r"""
     Given a 3D grid, selects a slice from the grid and attempts to
-     solve for an eigenmode propagating through that slice.
+    solve for an eigenmode propagating through that slice.
 
     Args:
         mode_number: Number of the mode, 0-indexed
@@ -35,19 +52,22 @@ def solve_mode(
         axis: Propagation axis (0=x, 1=y, 2=z)
         polarity: Propagation direction (+1 for +ve, -1 for -ve)
         slices: `epsilon[tuple(slices)]` is used to select the portion of the grid to use
-                as the waveguide cross-section. `slices[axis]` should select only one item.
+                as the waveguide cross-section. `slices[axis]` must select exactly one item.
         epsilon: Dielectric constant
         mu: Magnetic permeability (default 1 everywhere)
 
     Returns:
-        ```
+        Dictionary containing:
-        {
+
-            'E': NDArray[complexfloating],
+        - `E`: full-grid electric field for the solved mode
-            'H': NDArray[complexfloating],
+        - `H`: full-grid magnetic field for the solved mode
-            'wavenumber': complex,
+        - `wavenumber`: propagation constant corrected for the discretized
-            'wavenumber_2d': complex,
+          propagation axis
-        }
+        - `wavenumber_2d`: propagation constant of the reduced 2D eigenproblem
-        ```
+
+    Notes:
+        The returned fields are normalized through the `waveguide_2d`
+        normalization convention before being expanded back to 3D.
     """
     if mu is None:
         mu = numpy.ones_like(epsilon)
@@ -137,7 +157,14 @@ def compute_source(
         mu: Magnetic permeability (default 1 everywhere)
 
     Returns:
-        J distribution for the unidirectional source
+        `J` distribution for a one-sided electric-current source.
+
+    Notes:
+        The source is built from the expanded mode field and a boundary-source
+        operator. The resulting current is intended to be injected with the
+        same sign convention used elsewhere in the package:
+
+        `E -= dt * J / epsilon`
     """
     E_expanded = expand_e(E=E, dxes=dxes, wavenumber=wavenumber, axis=axis,
                           polarity=polarity, slices=slices)
@@ -157,19 +184,52 @@ def compute_source(
 
 
 def compute_overlap_e(
-        E: cfdfield,
+        E: cfdfield_t,
         wavenumber: complex,
         dxes: dx_lists_t,
         axis: int,
         polarity: int,
         slices: Sequence[slice],
+        omega: float,
         ) -> cfdfield_t:
-    """
+    r"""
-    Given an eigenmode obtained by `solve_mode`, calculates an overlap_e for the
+    Build an overlap field for projecting another 3D electric field onto a mode.
-    mode orthogonality relation Integrate(((E x H_mode) + (E_mode x H)) dot dn)
+
-    [assumes reflection symmetry].
+    The returned field is intended for the discrete overlap expression
+
+    $$
+    \sum \mathrm{overlap\_e} \; E_\mathrm{other}^*
+    $$
+
+    where the sum is over the full Yee-grid field storage.
+
+    The construction uses a two-cell window immediately upstream of the selected
+    slice:
+
+    - for `polarity=+1`, the two cells just before `slices[axis].start`
+    - for `polarity=-1`, the two cells just after `slices[axis].stop`
+
+    The window is clipped to the simulation domain if necessary. A clipped but
+    non-empty window raises `RuntimeWarning`; an empty window raises
+    `ValueError`.
+
+    The derivation below assumes reflection symmetry and the standard waveguide
+    overlap relation involving
+
+    $$
+    \int ((E \times H_\mathrm{mode}) + (E_\mathrm{mode} \times H)) \cdot dn.
+    $$
+
+    E x H_mode + E_mode x H
+    -> Ex Hmy - EyHmx + Emx Hy - Emy Hx  (Z-prop)
+    Ex we/B Emx + Ex i/B dy Hmz - Ey (-we/B Emy) - Ey i/B dx Hmz
+    we/B (Ex Emx + Ey Emy) + i/B (Ex dy Hmz - Ey dx Hmz)
+    we/B (Ex Emx + Ey Emy) + i/B (Ex dy (dx Emy - dy Emx) - Ey dx (dx Emy - dy Emx))
+    we/B (Ex Emx + Ey Emy) + i/B (Ex dy dx Emy - Ex dy dy Emx - Ey dx dx Emy -  Ey dx dy Emx)
+
+    Ex j/wu (-jB Emx - dx Emz) - Ey j/wu (dy Emz + jB Emy)
+    B/wu (Ex Emx + Ey Emy) - j/wu (Ex dx Emz + Ey dy Emz)
 
-    TODO: add reference or derivation for compute_overlap_e
 
     Args:
         E: E-field of the mode
@@ -181,25 +241,41 @@ def compute_overlap_e(
                 as the waveguide cross-section. slices[axis] should select only one item.
 
     Returns:
-        overlap_e such that `numpy.sum(overlap_e * other_e.conj())` computes the overlap integral
+        `overlap_e` normalized so that `numpy.sum(overlap_e * E.conj()) == 1`
+        over the retained overlap window.
     """
     slices = tuple(slices)
 
     Ee = expand_e(E=E, wavenumber=wavenumber, dxes=dxes,
                   axis=axis, polarity=polarity, slices=slices)
 
-    start, stop = sorted((slices[axis].start, slices[axis].start - 2 * polarity))
+    axis_size = E.shape[axis + 1]
+    if polarity > 0:
+        start = slices[axis].start - 2
+        stop = slices[axis].start
+    else:
+        start = slices[axis].stop
+        stop = slices[axis].stop + 2
+
+    clipped_start = max(0, start)
+    clipped_stop = min(axis_size, stop)
+    if clipped_start >= clipped_stop:
+        raise ValueError('Requested overlap window lies outside the domain')
+    if clipped_start != start or clipped_stop != stop:
+        warnings.warn('Requested overlap window was clipped to fit within the domain', RuntimeWarning)
 
     slices2_l = list(slices)
-    slices2_l[axis] = slice(start, stop)
+    slices2_l[axis] = slice(clipped_start, clipped_stop)
     slices2 = (slice(None), *slices2_l)
 
     Etgt = numpy.zeros_like(Ee)
     Etgt[slices2] = Ee[slices2]
 
-    # note no sqrt() when normalizing below since we want to get 1.0 after overlapping with the
+    # Note: We normalize so that (Etgt @ E.conj()) == 1, so (Etgt @ Etgt.conj) != 1
-    # original field, not the normalized one
+    norm = (Etgt.conj() * Etgt).sum()
-    Etgt /= (Etgt.conj() * Etgt).sum()    # type: ignore
+    if norm == 0:
+        raise ValueError('Requested overlap window contains no overlap field support')
+    Etgt /= norm
     return cfdfield_t(Etgt)
 
 
@@ -211,7 +287,7 @@ def expand_e(
         polarity: int,
         slices: Sequence[slice],
         ) -> cfdfield_t:
-    """
+    r"""
     Given an eigenmode obtained by `solve_mode`, expands the E-field from the 2D
     slice where the mode was calculated to the entire domain (along the propagation
     axis). This assumes the epsilon cross-section remains constant throughout the
@@ -229,6 +305,16 @@ def expand_e(
 
     Returns:
         `E`, with the original field expanded along the specified `axis`.
+
+    Notes:
+        This helper assumes that the waveguide cross-section remains constant
+        along the propagation axis and applies the phase factor
+
+        $$
+        e^{-i \, \mathrm{polarity} \, wavenumber \, \Delta z}
+        $$
+
+        to each copied slice.
     """
     slices = tuple(slices)
 

meanas/fdfd/waveguide_cyl.py

@@ -2,63 +2,125 @@ r"""
 Operators and helper functions for cylindrical waveguides with unchanging cross-section.
 
 Waveguide operator is derived according to 10.1364/OL.33.001848.
-The curl equations in the complex coordinate system become
+
+As in `waveguide_2d`, the propagation dependence is separated from the
+transverse solve. Here the propagation coordinate is the bend angle `\theta`,
+and the fields are assumed to have the form
+
+$$
+\vec{E}(r, y, \theta), \vec{H}(r, y, \theta) \propto e^{-\imath m \theta},
+$$
+
+where `m` is the angular wavenumber returned by `solve_mode(s)`. It is often
+convenient to introduce the corresponding linear wavenumber
+
+$$
+\beta = \frac{m}{r_{\min}},
+$$
+
+so that the cylindrical problem resembles the straight-waveguide problem with
+additional metric factors.
+
+Those metric factors live on the staggered radial Yee grids. If the left edge of
+the computational window is at `r = r_{\min}`, define the electric-grid and
+magnetic-grid radial sample locations by
 
 $$
 \begin{aligned}
--\imath \omega \mu_{xx} H_x &= \tilde{\partial}_y E_z + \imath \beta frac{E_y}{\tilde{t}_x} \\
+r_a(n) &= r_{\min} + \sum_{j \le n} \Delta r_{e, j}, \\
--\imath \omega \mu_{yy} H_y &= -\imath \beta E_x - \frac{1}{\hat{t}_x} \tilde{\partial}_x \tilde{t}_x E_z \\
+r_b\!\left(n + \tfrac{1}{2}\right) &= r_{\min} + \tfrac{1}{2}\Delta r_{e, n}
--\imath \omega \mu_{zz} H_z &= \tilde{\partial}_x E_y - \tilde{\partial}_y E_x \\
+  + \sum_{j < n} \Delta r_{h, j},
-\imath \omega \epsilon_{xx} E_x &= \hat{\partial}_y H_z + \imath \beta \frac{H_y}{\hat{T}} \\
-\imath \omega \epsilon_{yy} E_y &= -\imath \beta H_x - \{1}{\tilde{t}_x} \hat{\partial}_x \hat{t}_x} H_z \\
-\imath \omega \epsilon_{zz} E_z &= \hat{\partial}_x H_y - \hat{\partial}_y H_x \\
 \end{aligned}
 $$
 
-where $t_x = 1 + \frac{\Delta_{x, m}}{R_0}$ is the grid spacing adjusted by the nominal radius $R0$.
+and from them the diagonal metric matrices
-
-Rewrite the last three equations as
 
 $$
 \begin{aligned}
-\imath \beta H_y &=  \imath \omega \hat{t}_x \epsilon_{xx} E_x - \hat{t}_x \hat{\partial}_y H_z \\
+T_a &= \operatorname{diag}(r_a / r_{\min}), \\
-\imath \beta H_x &= -\imath \omega \hat{t}_x \epsilon_{yy} E_y - \hat{t}_x \hat{\partial}_x H_z \\
+T_b &= \operatorname{diag}(r_b / r_{\min}).
-\imath \omega E_z &= \frac{1}{\epsilon_{zz}} \hat{\partial}_x H_y - \frac{1}{\epsilon_{zz}} \hat{\partial}_y H_x \\
 \end{aligned}
 $$
 
-The derivation then follows the same steps as the straight waveguide, leading to the eigenvalue problem
+With the same forward/backward derivative notation used in `waveguide_2d`, the
-
+coordinate-transformed discrete curl equations used here are
-$$
-\beta^2 \begin{bmatrix} E_x \\
-                        E_y \end{bmatrix} =
-    (\omega^2 \begin{bmatrix} T_b T_b \mu_{yy} \epsilon_{xx} & 0 \\
-                                                            0 & T_a T_a \mu_{xx} \epsilon_{yy} \end{bmatrix} +
-              \begin{bmatrix} -T_b \mu_{yy} \hat{\partial}_y \\
-                               T_a \mu_{xx} \hat{\partial}_x \end{bmatrix} T_b \mu_{zz}^{-1}
-              \begin{bmatrix} -\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} +
-      \begin{bmatrix} \tilde{\partial}_x \\
-                      \tilde{\partial}_y \end{bmatrix} T_a \epsilon_{zz}^{-1}
-                 \begin{bmatrix} \hat{\partial}_x T_b \epsilon_{xx} & \hat{\partial}_y T_a \epsilon_{yy} \end{bmatrix})
-    \begin{bmatrix} E_x \\
-                    E_y \end{bmatrix}
-$$
-
-which resembles the straight waveguide eigenproblem with additonal $T_a$ and $T_b$ terms. These
-are diagonal matrices containing the $t_x$ values:
 
 $$
 \begin{aligned}
-T_a &=  1 + \frac{\Delta_{x, m               }}{R_0}
+-\imath \omega \mu_{rr} H_r &= \tilde{\partial}_y E_z + \imath \beta T_a^{-1} E_y, \\
-T_b &=  1 + \frac{\Delta_{x, m + \frac{1}{2} }}{R_0}
+-\imath \omega \mu_{yy} H_y &= -\imath \beta T_b^{-1} E_r
+  - T_b^{-1} \tilde{\partial}_r (T_a E_z), \\
+-\imath \omega \mu_{zz} H_z &= \tilde{\partial}_r E_y - \tilde{\partial}_y E_r, \\
+\imath \beta H_y &= -\imath \omega T_b \epsilon_{rr} E_r - T_b \hat{\partial}_y H_z, \\
+\imath \beta H_r &= \imath \omega T_a \epsilon_{yy} E_y
+  - T_b T_a^{-1} \hat{\partial}_r (T_b H_z), \\
+\imath \omega E_z &= T_a \epsilon_{zz}^{-1}
+  \left(\hat{\partial}_r H_y - \hat{\partial}_y H_r\right).
 \end{aligned}
-
-
-TODO: consider 10.1364/OE.20.021583 for an alternate approach
 $$
 
-As in the straight waveguide case, all the functions in this file assume a 2D grid
+The first three equations are the cylindrical analogue of the straight-guide
- (i.e. `dxes = [[[dr_e_0, dx_e_1, ...], [dy_e_0, ...]], [[dr_h_0, ...], [dy_h_0, ...]]]`).
+relations for `H_r`, `H_y`, and `H_z`. The next two are the metric-weighted
+versions of the straight-guide identities for `\imath \beta H_y` and
+`\imath \beta H_r`, and the last equation plays the same role as the
+longitudinal `E_z` reconstruction in `waveguide_2d`.
+
+Following the same elimination steps as in `waveguide_2d`, apply
+`\imath \beta \tilde{\partial}_r` and `\imath \beta \tilde{\partial}_y` to the
+equation for `E_z`, substitute for `\imath \beta H_r` and `\imath \beta H_y`,
+and then eliminate `H_z` with
+
+$$
+H_z = \frac{1}{-\imath \omega \mu_{zz}}
+\left(\tilde{\partial}_r E_y - \tilde{\partial}_y E_r\right).
+$$
+
+This yields the transverse electric eigenproblem implemented by
+`cylindrical_operator(...)`:
+
+$$
+\beta^2
+\begin{bmatrix} E_r \\ E_y \end{bmatrix}
+=
+\left(
+\omega^2
+\begin{bmatrix}
+T_b^2 \mu_{yy} \epsilon_{xx} & 0 \\
+0 & T_a^2 \mu_{xx} \epsilon_{yy}
+\end{bmatrix}
++
+\begin{bmatrix}
+-T_b \mu_{yy} \hat{\partial}_y \\
+ T_a \mu_{xx} \hat{\partial}_x
+\end{bmatrix}
+T_b \mu_{zz}^{-1}
+\begin{bmatrix}
+-\tilde{\partial}_y & \tilde{\partial}_x
+\end{bmatrix}
++
+\begin{bmatrix}
+\tilde{\partial}_x \\
+\tilde{\partial}_y
+\end{bmatrix}
+T_a \epsilon_{zz}^{-1}
+\begin{bmatrix}
+\hat{\partial}_x T_b \epsilon_{xx} &
+\hat{\partial}_y T_a \epsilon_{yy}
+\end{bmatrix}
+\right)
+\begin{bmatrix} E_r \\ E_y \end{bmatrix}.
+$$
+
+Since `\beta = m / r_{\min}`, the solver implemented in this file returns the
+angular wavenumber `m`, while the operator itself is most naturally written in
+terms of the linear quantity `\beta`. The helpers below reconstruct the full
+field components from the solved transverse eigenvector and then normalize the
+mode to unit forward power with the same discrete longitudinal Poynting inner
+product used by `waveguide_2d`.
+
+As in the straight-waveguide case, all functions here assume a 2D grid:
+
+`dxes = [[[dr_e_0, dr_e_1, ...], [dy_e_0, ...]], [[dr_h_0, ...], [dy_h_0, ...]]]`.
 """
 from typing import Any, cast
 from collections.abc import Sequence
@@ -94,17 +156,18 @@ def cylindrical_operator(
           \begin{bmatrix} \tilde{\partial}_x \\
                           \tilde{\partial}_y \end{bmatrix} T_a \epsilon_{zz}^{-1}
                      \begin{bmatrix} \hat{\partial}_x T_b \epsilon_{xx} & \hat{\partial}_y T_a \epsilon_{yy} \end{bmatrix})
-        \begin{bmatrix} E_x \\
+        \begin{bmatrix} E_r \\
                         E_y \end{bmatrix}
     $$
 
     for use with a field vector of the form `[E_r, E_y]`.
 
     This operator can be used to form an eigenvalue problem of the form
-        A @ [E_r, E_y] = wavenumber**2 * [E_r, E_y]
+        A @ [E_r, E_y] = beta**2 * [E_r, E_y]
 
     which can then be solved for the eigenmodes of the system
-    (an `exp(-i * wavenumber * theta)` theta-dependence is assumed for the fields).
+    (an `exp(-i * angular_wavenumber * theta)` theta-dependence is assumed for
+    the fields, with `beta = angular_wavenumber / rmin`).
 
     (NOTE: See module docs and 10.1364/OL.33.001848)
 
@@ -270,7 +333,7 @@ def exy2h(
         mu: vfdslice | None = None
         ) -> sparse.sparray:
     """
-    Operator which transforms the vector `e_xy` containing the vectorized E_x and E_y fields,
+    Operator which transforms the vector `e_xy` containing the vectorized E_r and E_y fields,
      into a vectorized H containing all three H components
 
     Args:
@@ -298,11 +361,11 @@ def exy2e(
         epsilon: vfdslice,
         ) -> sparse.sparray:
     """
-    Operator which transforms the vector `e_xy` containing the vectorized E_x and E_y fields,
+    Operator which transforms the vector `e_xy` containing the vectorized E_r and E_y fields,
      into a vectorized E containing all three E components
 
     Unlike the straight waveguide case, the H_z components do not cancel and must be calculated
-    from E_x and E_y in order to then calculate E_z.
+    from E_r and E_y in order to then calculate E_z.
 
     Args:
         angular_wavenumber: Wavenumber assuming fields have theta-dependence of
@@ -360,9 +423,10 @@ def e2h(
     This operator is created directly from the initial coordinate-transformed equations:
     $$
     \begin{aligned}
-    \imath \omega \epsilon_{xx} E_x &= \hat{\partial}_y H_z + \imath \beta \frac{H_y}{\hat{T}} \\
+    -\imath \omega \mu_{rr} H_r &= \tilde{\partial}_y E_z + \imath \beta T_a^{-1} E_y, \\
-    \imath \omega \epsilon_{yy} E_y &= -\imath \beta H_x - \{1}{\tilde{t}_x} \hat{\partial}_x \hat{t}_x} H_z \\
+    -\imath \omega \mu_{yy} H_y &= -\imath \beta T_b^{-1} E_r
-    \imath \omega \epsilon_{zz} E_z &= \hat{\partial}_x H_y - \hat{\partial}_y H_x \\
+      - T_b^{-1} \tilde{\partial}_r (T_a E_z), \\
+    -\imath \omega \mu_{zz} H_z &= \tilde{\partial}_r E_y - \tilde{\partial}_y E_r,
     \end{aligned}
     $$
 
@@ -397,15 +461,18 @@ def dxes2T(
         rmin: float,
         ) -> tuple[NDArray[numpy.float64], NDArray[numpy.float64]]:
     r"""
-    Returns the $T_a$ and $T_b$ diagonal matrices which are used to apply the cylindrical
+    Construct the cylindrical metric matrices $T_a$ and $T_b$.
-      coordinate transformation in various operators.
+
+    `T_a` is sampled on the E-grid radial locations, while `T_b` is sampled on
+    the staggered H-grid radial locations. These are the diagonal matrices that
+    convert the straight-waveguide algebra into its cylindrical counterpart.
 
     Args:
         dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
         rmin: Radius at the left edge of the simulation domain (at minimum 'x')
 
     Returns:
-        Sparse matrix representations of the operators Ta and Tb
+        Sparse diagonal matrices `(T_a, T_b)`.
     """
     ra = rmin + numpy.cumsum(dxes[0][0])                      # Radius at Ey points
     rb = rmin + dxes[0][0] / 2.0 + numpy.cumsum(dxes[1][0])   # Radius at Ex points
@@ -427,12 +494,12 @@ def normalized_fields_e(
         mu: vfdslice | None = None,
         prop_phase: float = 0,
         ) -> tuple[vcfdslice_t, vcfdslice_t]:
-    """
+    r"""
-    Given a vector `e_xy` containing the vectorized E_x and E_y fields,
+    Given a vector `e_xy` containing the vectorized E_r and E_y fields,
-     returns normalized, vectorized E and H fields for the system.
+    returns normalized, vectorized E and H fields for the system.
 
     Args:
-        e_xy: Vector containing E_x and E_y fields
+        e_xy: Vector containing E_r and E_y fields
         angular_wavenumber: Wavenumber assuming fields have theta-dependence of
             `exp(-i * angular_wavenumber * theta)`. It should satisfy
             `operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy`
@@ -447,6 +514,18 @@ def normalized_fields_e(
     Returns:
         `(e, h)`, where each field is vectorized, normalized,
         and contains all three vector components.
+
+    Notes:
+        The normalization step is delegated to `_normalized_fields(...)`, which
+        enforces unit forward power under the discrete inner product
+
+        $$
+        \frac{1}{2}\int (E_r H_y^* - E_y H_r^*) \, dr \, dy.
+        $$
+
+        The angular wavenumber `m` is first converted into the full three-component
+        fields, then the overall complex phase and sign are fixed so the result is
+        reproducible for symmetric modes.
     """
     e = exy2e(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon) @ e_xy
     h = exy2h(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon, mu=mu) @ e_xy
@@ -465,8 +544,30 @@ def _normalized_fields(
         mu: vfdslice | None = None,
         prop_phase: float = 0,
         ) -> tuple[vcfdslice_t, vcfdslice_t]:
+    r"""
+    Normalize a cylindrical waveguide mode to unit forward power.
+
+    The cylindrical helpers reuse the straight-waveguide inner product after the
+    field reconstruction step. The extra metric factors have already been folded
+    into the reconstructed `e`/`h` fields through `dxes2T(...)` and the
+    cylindrical `exy2e(...)` / `exy2h(...)` operators, so the same discrete
+    longitudinal Poynting integral can be used here.
+
+    The normalization procedure is:
+
+    1. Flip the magnetic field sign so the reconstructed `(e, h)` pair follows
+       the same forward-power convention as `waveguide_2d`.
+    2. Compute the time-averaged forward power with
+       `waveguide_2d.inner_product(..., conj_h=True)`.
+    3. Scale by `1 / sqrt(S_z)` so the resulting mode has unit forward power.
+    4. Remove the arbitrary complex phase and apply a quadrant-sum sign heuristic
+       so symmetric modes choose a stable representative.
+
+    `prop_phase` has the same meaning as in `waveguide_2d`: it compensates for
+    the half-cell longitudinal staggering between the E and H fields when the
+    propagation direction is itself discretized.
+    """
     h *= -1
-    # TODO documentation for normalized_fields
     shape = [s.size for s in dxes[0]]
 
     # Find time-averaged Sz and normalize to it

meanas/fdmath/vectorization.py

@@ -18,35 +18,35 @@ from .types import (
 
 @overload
 def vec(f: None) -> None:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def vec(f: fdfield_t) -> vfdfield_t:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def vec(f: cfdfield_t) -> vcfdfield_t:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def vec(f: fdfield2_t) -> vfdfield2_t:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def vec(f: cfdfield2_t) -> vcfdfield2_t:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def vec(f: fdslice_t) -> vfdslice_t:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def vec(f: cfdslice_t) -> vcfdslice_t:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def vec(f: ArrayLike) -> NDArray:
-    pass
+    pass  # pragma: no cover
 
 def vec(
         f: fdfield_t | cfdfield_t | fdfield2_t | cfdfield2_t | fdslice_t | cfdslice_t | ArrayLike | None,
@@ -70,15 +70,15 @@ def vec(
 
 @overload
 def unvec(v: None, shape: Sequence[int], nvdim: int = 3) -> None:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def unvec(v: vfdfield_t, shape: Sequence[int], nvdim: int = 3) -> fdfield_t:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def unvec(v: vcfdfield_t, shape: Sequence[int], nvdim: int = 3) -> cfdfield_t:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def unvec(v: vfdfield2_t, shape: Sequence[int], nvdim: int = 3) -> fdfield2_t:

meanas/fdtd/__init__.py

@@ -144,6 +144,18 @@ It is often useful to excite the simulation with an arbitrary broadband pulse an
 extract the frequency-domain response by performing an on-the-fly Fourier transform
 of the time-domain fields.
 
+`accumulate_phasor` in `meanas.fdtd.phasor` performs the phase accumulation for one
+or more target frequencies, while leaving source normalization and simulation-loop
+policy to the caller. Convenience wrappers `accumulate_phasor_e`,
+`accumulate_phasor_h`, and `accumulate_phasor_j` apply the usual Yee time offsets.
+The helpers accumulate
+
+$$ \Delta_t \sum_l w_l e^{-i \omega t_l} f_l $$
+
+with caller-provided sample times and weights. In this codebase the matching
+electric-current convention is typically `E -= dt * J / epsilon` in FDTD and
+`-i \omega J` on the right-hand side of the FDFD wave equation.
+
 The Ricker wavelet (normalized second derivative of a Gaussian) is commonly used for the pulse
 shape. It can be written
 
@@ -156,7 +168,44 @@ t=0 (assuming the source is off for t<0 this gives $\sim 10^{-3}$ error at t=0).
 
 Boundary conditions
 ===================
-# TODO notes about boundaries / PMLs
+
+`meanas.fdtd` exposes two boundary-related building blocks:
+
+- `conducting_boundary(...)` for simple perfect-electric-conductor style field
+  clamping at one face of the domain.
+- `cpml_params(...)` and `updates_with_cpml(...)` for convolutional perfectly
+  matched layers (CPMLs) attached to one or more faces of the Yee grid.
+
+`updates_with_cpml(...)` accepts a three-by-two table of CPML parameter blocks:
+
+```
+cpml_params[axis][polarity_index]
+```
+
+where `axis` is `0`, `1`, or `2` and `polarity_index` corresponds to `(-1, +1)`.
+Passing `None` for one entry disables CPML on that face while leaving the other
+directions unchanged. This is how mixed boundary setups such as "absorbing in x,
+periodic in y/z" are expressed.
+
+When comparing an FDTD run against an FDFD solve, use the same stretched
+coordinate system in both places:
+
+1. Build the FDTD update with the desired CPML parameters.
+2. Stretch the FDFD `dxes` with the matching SCPML transform.
+3. Compare the extracted phasor against the FDFD field or residual on those
+   stretched `dxes`.
+
+The electric-current sign convention used throughout the examples and tests is
+
+$$
+E \leftarrow E - \Delta_t J / \epsilon
+$$
+
+which matches the FDFD right-hand side
+
+$$
+-i \omega J.
+$$
 """
 
 from .base import (
@@ -178,3 +227,9 @@ from .energy import (
 from .boundaries import (
     conducting_boundary as conducting_boundary,
     )
+from .phasor import (
+    accumulate_phasor as accumulate_phasor,
+    accumulate_phasor_e as accumulate_phasor_e,
+    accumulate_phasor_h as accumulate_phasor_h,
+    accumulate_phasor_j as accumulate_phasor_j,
+    )

meanas/fdtd/boundaries.py

@@ -28,17 +28,19 @@ def conducting_boundary(
         shifted1_slice = [slice(None)] * 3
         boundary_slice[direction] = 0
         shifted1_slice[direction] = 1
+        boundary = tuple(boundary_slice)
+        shifted1 = tuple(shifted1_slice)
 
         def en(e: fdfield_t) -> fdfield_t:
-            e[direction][boundary_slice] = 0
+            e[direction][boundary] = 0
-            e[u][boundary_slice] = e[u][shifted1_slice]
+            e[u][boundary] = e[u][shifted1]
-            e[v][boundary_slice] = e[v][shifted1_slice]
+            e[v][boundary] = e[v][shifted1]
             return e
 
         def hn(h: fdfield_t) -> fdfield_t:
-            h[direction][boundary_slice] = h[direction][shifted1_slice]
+            h[direction][boundary] = h[direction][shifted1]
-            h[u][boundary_slice] = 0
+            h[u][boundary] = 0
-            h[v][boundary_slice] = 0
+            h[v][boundary] = 0
             return h
 
         return en, hn
@@ -50,20 +52,23 @@ def conducting_boundary(
         boundary_slice[direction] = -1
         shifted1_slice[direction] = -2
         shifted2_slice[direction] = -3
+        boundary = tuple(boundary_slice)
+        shifted1 = tuple(shifted1_slice)
+        shifted2 = tuple(shifted2_slice)
 
         def ep(e: fdfield_t) -> fdfield_t:
-            e[direction][boundary_slice] = -e[direction][shifted2_slice]
+            e[direction][boundary] = -e[direction][shifted2]
-            e[direction][shifted1_slice] = 0
+            e[direction][shifted1] = 0
-            e[u][boundary_slice] = e[u][shifted1_slice]
+            e[u][boundary] = e[u][shifted1]
-            e[v][boundary_slice] = e[v][shifted1_slice]
+            e[v][boundary] = e[v][shifted1]
             return e
 
         def hp(h: fdfield_t) -> fdfield_t:
-            h[direction][boundary_slice] = h[direction][shifted1_slice]
+            h[direction][boundary] = h[direction][shifted1]
-            h[u][boundary_slice] = -h[u][shifted2_slice]
+            h[u][boundary] = -h[u][shifted2]
-            h[u][shifted1_slice] = 0
+            h[u][shifted1] = 0
-            h[v][boundary_slice] = -h[v][shifted2_slice]
+            h[v][boundary] = -h[v][shifted2]
-            h[v][shifted1_slice] = 0
+            h[v][shifted1] = 0
             return h
 
         return ep, hp

meanas/fdtd/energy.py

@@ -4,7 +4,21 @@ from ..fdmath import dx_lists_t, fdfield_t, fdfield
 from ..fdmath.functional import deriv_back
 
 
-# TODO documentation
+"""
+Energy- and flux-accounting helpers for Yee-grid FDTD fields.
+
+These functions complement the derivation in `meanas.fdtd`:
+
+- `poynting(...)` and `poynting_divergence(...)` evaluate the discrete flux terms
+  from the exact time-domain Poynting identity.
+- `energy_hstep(...)` / `energy_estep(...)` evaluate the two staggered energy
+  expressions.
+- `delta_energy_*` helpers evaluate the corresponding energy changes between
+  adjacent half-steps.
+
+The helpers are intended for diagnostics, regression tests, and consistency
+checks between source work, field energy, and flux through cell faces.
+"""
 
 
 def poynting(
@@ -252,13 +266,23 @@ def delta_energy_j(
         e1: fdfield,
         dxes: dx_lists_t | None = None,
         ) -> fdfield_t:
-    """
+    r"""
-    Calculate
+    Calculate the electric-current work term $J \cdot E$ on the Yee grid.
 
-    Note that each value of $J$ contributes to the energy twice (i.e. once per field update)
+    This is the source contribution that appears beside the flux divergence in
-    despite only causing the value of $E$ to change once (same for $M$ and $H$).
+    the discrete Poynting identities documented in `meanas.fdtd`.
 
+    Note that each value of `J` contributes twice in a full Yee cycle (once per
+    half-step energy balance) even though it directly changes `E` only once.
 
+    Args:
+        j0: Electric-current density sampled at the same half-step as the
+            current work term.
+        e1: Electric field sampled at the matching integer timestep.
+        dxes: Grid description; defaults to unit spacing.
+
+    Returns:
+        Per-cell source-work contribution with the scalar field shape.
     """
     if dxes is None:
         dxes = tuple(tuple(numpy.ones(1) for _ in range(3)) for _ in range(2))
@@ -277,6 +301,20 @@ def dxmul(
         mu: fdfield | float | None = None,
         dxes: dx_lists_t | None = None,
         ) -> fdfield_t:
+    """
+    Multiply E- and H-like field products by material weights and cell volumes.
+
+    Args:
+        ee: Three-component electric-field product, such as `e0 * e2`.
+        hh: Three-component magnetic-field product, such as `h1 * h1`.
+        epsilon: Electric material weight; defaults to `1`.
+        mu: Magnetic material weight; defaults to `1`.
+        dxes: Grid description; defaults to unit spacing.
+
+    Returns:
+        Scalar field containing the weighted electric plus magnetic contribution
+        for each Yee cell.
+    """
     if epsilon is None:
         epsilon = 1
     if mu is None:

meanas/fdtd/misc.py

@@ -53,8 +53,8 @@ def gaussian_packet(
         t0 = ii * dt - delay
         envelope = numpy.sqrt(numpy.sqrt(2 * alpha / pi)) * numpy.exp(-alpha * t0 * t0)
 
-        if one_sided and t0 > 0:
+        if one_sided:
-            envelope = 1
+            envelope = numpy.where(t0 > 0, 1.0, envelope)
 
         cc = numpy.cos(omega * t0)
         ss = numpy.sin(omega * t0)
@@ -94,10 +94,14 @@ def ricker_pulse(
     logger.warning('meanas.fdtd.misc functions are still very WIP!')    # TODO
     omega = 2 * pi / wl
     freq = 1 / wl
-    # r0 = omega / 2
 
     from scipy.optimize import root_scalar
-    delay_results = root_scalar(lambda tt: (1 - omega * omega * tt * tt / 2) * numpy.exp(-omega * omega / 4 * tt * tt) - turn_on, x0=0, x1=-2 / omega)
+    delay_results = root_scalar(
+        lambda tt: (1 - omega * omega * tt * tt / 2) * numpy.exp(-omega * omega * tt * tt / 4) - turn_on,
+        x0=0,
+        x1=-2 / omega,
+        )
+
     delay = delay_results.root
     delay = numpy.ceil(delay * freq) / freq     # force delay to integer number of periods to maintain phase
 
@@ -156,11 +160,11 @@ def gaussian_beam(
     wz = numpy.sqrt(wz2)        # == fwhm(z) / sqrt(2 * ln(2))
 
     kk = 2 * pi / wl
-    Rz = zz * (1 + zr2 / z2)
+    inv_Rz = numpy.divide(zz, z2 + zr2, out=numpy.zeros_like(zz), where=(z2 + zr2) != 0)
     gouy = numpy.arctan(zz / zr)
 
-    gaussian = w0 / wz * numpy.exp(-r2 / wz2) * numpy.exp(1j * (kk * zz + kk * r2 / 2 / Rz - gouy))
+    gaussian = w0 / wz * numpy.exp(-r2 / wz2) * numpy.exp(1j * (kk * zz + kk * r2 * inv_Rz / 2 - gouy))
 
     row = gaussian[:, :, gaussian.shape[2] // 2]
-    norm = numpy.sqrt((row * row.conj()).sum())
+    norm = numpy.linalg.norm(row)
     return gaussian / norm

meanas/fdtd/phasor.py

@@ -0,0 +1,126 @@
+"""
+Helpers for extracting single- or multi-frequency phasors from FDTD samples.
+
+These helpers are intentionally low-level: callers own the accumulator arrays and
+the sampling policy. The accumulated quantity is
+
+    dt * sum(weight * exp(-1j * omega * t_step) * sample_step)
+
+where `t_step = (step + offset_steps) * dt`.
+
+The usual Yee offsets are:
+
+- `accumulate_phasor_e(..., step=l)` for `E_l`
+- `accumulate_phasor_h(..., step=l)` for `H_{l + 1/2}`
+- `accumulate_phasor_j(..., step=l)` for `J_{l + 1/2}`
+
+These helpers do not choose warmup/accumulation windows or pulse-envelope
+normalization. They also do not impose a current sign convention. In this
+codebase, electric-current injection normally appears as `E -= dt * J / epsilon`,
+which matches the FDFD right-hand side `-1j * omega * J`.
+"""
+from collections.abc import Sequence
+
+import numpy
+from numpy.typing import ArrayLike, NDArray
+
+
+def _normalize_omegas(
+        omegas: float | complex | Sequence[float | complex] | NDArray,
+        ) -> NDArray[numpy.complexfloating]:
+    omega_array = numpy.atleast_1d(numpy.asarray(omegas, dtype=complex))
+    if omega_array.ndim != 1 or omega_array.size == 0:
+        raise ValueError('omegas must be a scalar or non-empty 1D sequence')
+    return omega_array
+
+
+def _normalize_weight(
+        weight: ArrayLike,
+        omega_shape: tuple[int, ...],
+        ) -> NDArray[numpy.complexfloating]:
+    weight_array = numpy.asarray(weight, dtype=complex)
+    if weight_array.ndim == 0:
+        return numpy.full(omega_shape, weight_array, dtype=complex)
+    if weight_array.shape == omega_shape:
+        return weight_array
+    raise ValueError(f'weight must be scalar or have shape {omega_shape}, got {weight_array.shape}')
+
+
+def accumulate_phasor(
+        accumulator: NDArray[numpy.complexfloating],
+        omegas: float | complex | Sequence[float | complex] | NDArray,
+        dt: float,
+        sample: ArrayLike,
+        step: int,
+        *,
+        offset_steps: float = 0.0,
+        weight: ArrayLike = 1.0,
+        ) -> NDArray[numpy.complexfloating]:
+    """
+    Add one time-domain sample into a phasor accumulator.
+
+    The added quantity is
+
+        dt * weight * exp(-1j * omega * t_step) * sample
+
+    where `t_step = (step + offset_steps) * dt`.
+
+    Note:
+        This helper already multiplies by `dt`. If the caller's normalization
+        factor was derived from a discrete sum that already includes `dt`, pass
+        `weight / dt` here.
+    """
+    if dt <= 0:
+        raise ValueError('dt must be positive')
+
+    omega_array = _normalize_omegas(omegas)
+    sample_array = numpy.asarray(sample)
+    expected_shape = (omega_array.size, *sample_array.shape)
+    if accumulator.shape != expected_shape:
+        raise ValueError(f'accumulator must have shape {expected_shape}, got {accumulator.shape}')
+
+    weight_array = _normalize_weight(weight, omega_array.shape)
+    time = (step + offset_steps) * dt
+    phase = numpy.exp(-1j * omega_array * time)
+    scaled = dt * (weight_array * phase).reshape((-1,) + (1,) * sample_array.ndim)
+    accumulator += scaled * sample_array
+    return accumulator
+
+
+def accumulate_phasor_e(
+        accumulator: NDArray[numpy.complexfloating],
+        omegas: float | complex | Sequence[float | complex] | NDArray,
+        dt: float,
+        sample: ArrayLike,
+        step: int,
+        *,
+        weight: ArrayLike = 1.0,
+        ) -> NDArray[numpy.complexfloating]:
+    """Accumulate an E-field sample taken at integer timestep `step`."""
+    return accumulate_phasor(accumulator, omegas, dt, sample, step, offset_steps=0.0, weight=weight)
+
+
+def accumulate_phasor_h(
+        accumulator: NDArray[numpy.complexfloating],
+        omegas: float | complex | Sequence[float | complex] | NDArray,
+        dt: float,
+        sample: ArrayLike,
+        step: int,
+        *,
+        weight: ArrayLike = 1.0,
+        ) -> NDArray[numpy.complexfloating]:
+    """Accumulate an H-field sample corresponding to `H_{step + 1/2}`."""
+    return accumulate_phasor(accumulator, omegas, dt, sample, step, offset_steps=0.5, weight=weight)
+
+
+def accumulate_phasor_j(
+        accumulator: NDArray[numpy.complexfloating],
+        omegas: float | complex | Sequence[float | complex] | NDArray,
+        dt: float,
+        sample: ArrayLike,
+        step: int,
+        *,
+        weight: ArrayLike = 1.0,
+        ) -> NDArray[numpy.complexfloating]:
+    """Accumulate a current sample corresponding to `J_{step + 1/2}`."""
+    return accumulate_phasor(accumulator, omegas, dt, sample, step, offset_steps=0.5, weight=weight)

meanas/fdtd/pml.py

@@ -1,9 +1,19 @@
 """
-PML implementations
+Convolutional perfectly matched layer (CPML) support for FDTD updates.
 
-#TODO discussion of PMLs
+The helpers in this module construct per-face CPML parameters and then wrap the
-#TODO cpml documentation
+standard Yee updates with the additional auxiliary `psi` fields needed by the
+CPML recurrence.
 
+The intended call pattern is:
+
+1. Build a `cpml_params[axis][polarity_index]` table with `cpml_params(...)`.
+2. Pass that table into `updates_with_cpml(...)` together with `dt`, `dxes`, and
+   `epsilon`.
+3. Advance the returned `update_E` / `update_H` closures in the simulation loop.
+
+Each face can be enabled or disabled independently by replacing one table entry
+with `None`.
 """
 # TODO retest pmls!
 
@@ -32,6 +42,29 @@ def cpml_params(
         ma: float = 1,
         cfs_alpha: float = 0,
         ) -> dict[str, Any]:
+    """
+    Construct the parameter block for one CPML face.
+
+    Args:
+        axis: Which Cartesian axis the CPML is normal to (`0`, `1`, or `2`).
+        polarity: Which face along that axis (`-1` for the low-index face,
+            `+1` for the high-index face).
+        dt: Timestep used by the Yee update.
+        thickness: Number of Yee cells occupied by the CPML region.
+        ln_R_per_layer: Logarithmic attenuation target per layer.
+        epsilon_eff: Effective permittivity used when choosing the CPML scaling.
+        mu_eff: Effective permeability used when choosing the CPML scaling.
+        m: Polynomial grading exponent for `sigma` and `kappa`.
+        ma: Polynomial grading exponent for the complex-frequency shift `alpha`.
+        cfs_alpha: Maximum complex-frequency shift parameter.
+
+    Returns:
+        Dictionary with:
+
+        - `param_e`: `(p0, p1, p2)` arrays for the E update
+        - `param_h`: `(p0, p1, p2)` arrays for the H update
+        - `region`: slice tuple selecting the CPML cells on that face
+    """
 
     if axis not in range(3):
         raise Exception(f'Invalid axis: {axis}')
@@ -57,8 +90,6 @@ def cpml_params(
         xh -= 0.5
         xe = xe[::-1]
         xh = xh[::-1]
-    else:
-        raise Exception('Bad polarity!')
 
     expand_slice_l: list[Any] = [None, None, None]
     expand_slice_l[axis] = slice(None)
@@ -82,8 +113,6 @@ def cpml_params(
         region_list[axis] = slice(None, thickness)
     elif polarity > 0:
         region_list[axis] = slice(-thickness, None)
-    else:
-        raise Exception('Bad polarity!')
     region = tuple(region_list)
 
     return {
@@ -102,6 +131,27 @@ def updates_with_cpml(
         dtype: DTypeLike = numpy.float32,
         ) -> tuple[Callable[[fdfield_t, fdfield_t, fdfield_t], None],
                    Callable[[fdfield_t, fdfield_t, fdfield_t], None]]:
+    """
+    Build Yee-step update closures augmented with CPML terms.
+
+    Args:
+        cpml_params: Three-by-two sequence indexed as `[axis][polarity_index]`.
+            Entries are the dictionaries returned by `cpml_params(...)`; use
+            `None` to disable CPML on one face.
+        dt: Timestep.
+        dxes: Yee-grid spacing lists `[dx_e, dx_h]`.
+        epsilon: Electric material distribution used by the E update.
+        dtype: Storage dtype for the auxiliary CPML state arrays.
+
+    Returns:
+        `(update_E, update_H)` closures with the same call shape as the basic
+        Yee updates:
+
+        - `update_E(e, h, epsilon)`
+        - `update_H(e, h, mu)`
+
+        The closures retain the CPML auxiliary state internally.
+    """
 
     Dfx, Dfy, Dfz = deriv_forward(dxes[1])
     Dbx, Dby, Dbz = deriv_back(dxes[1])

meanas/test/_bloch_case.py

@@ -0,0 +1,37 @@
+import numpy
+
+
+SHAPE = (2, 2, 2)
+G_MATRIX = numpy.eye(3)
+K0_GENERAL = numpy.array([0.1, 0.2, 0.3], dtype=float)
+K0_AXIAL = numpy.array([0.0, 0.0, 0.25], dtype=float)
+K0_X = numpy.array([0.1, 0.0, 0.0], dtype=float)
+EPSILON = numpy.ones((3, *SHAPE), dtype=float)
+MU = numpy.stack([
+    numpy.linspace(2.0, 2.7, numpy.prod(SHAPE)).reshape(SHAPE),
+    numpy.linspace(2.1, 2.8, numpy.prod(SHAPE)).reshape(SHAPE),
+    numpy.linspace(2.2, 2.9, numpy.prod(SHAPE)).reshape(SHAPE),
+])
+H_SIZE = 2 * numpy.prod(SHAPE)
+H_MN = (numpy.arange(H_SIZE) + 0.25j).astype(complex)
+ZERO_H_MN = numpy.zeros_like(H_MN)
+Y0 = (numpy.arange(H_SIZE, dtype=float) + 1j * numpy.linspace(0.1, 0.9, H_SIZE))[None, :]
+Y0_TWO_MODE = numpy.vstack(
+    [
+        numpy.arange(H_SIZE, dtype=float) + 1j * numpy.linspace(0.1, 0.9, H_SIZE),
+        numpy.linspace(2.0, 3.5, H_SIZE) - 0.5j * numpy.arange(H_SIZE, dtype=float),
+    ],
+)
+
+
+def build_overlap_fixture() -> tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray, numpy.ndarray]:
+    e_in = numpy.zeros((3, *SHAPE), dtype=complex)
+    h_in = numpy.zeros_like(e_in)
+    e_out = numpy.zeros_like(e_in)
+    h_out = numpy.zeros_like(e_in)
+
+    e_in[1] = 1.0
+    h_in[2] = 2.0
+    e_out[1] = 3.0
+    h_out[2] = 4.0
+    return e_in, h_in, e_out, h_out

meanas/test/_fdfd_case.py

@@ -0,0 +1,49 @@
+import numpy
+
+from ..fdmath import vec, unvec
+
+
+OMEGA = 1 / 1500
+SHAPE = (2, 3, 2)
+DXES = [
+    [numpy.array([1.0, 1.5]), numpy.array([0.75, 1.25, 1.5]), numpy.array([1.2, 0.8])],
+    [numpy.array([0.9, 1.4]), numpy.array([0.8, 1.1, 1.4]), numpy.array([1.0, 0.7])],
+]
+
+EPSILON = numpy.stack([
+    numpy.linspace(1.0, 2.2, numpy.prod(SHAPE)).reshape(SHAPE),
+    numpy.linspace(1.1, 2.3, numpy.prod(SHAPE)).reshape(SHAPE),
+    numpy.linspace(1.2, 2.4, numpy.prod(SHAPE)).reshape(SHAPE),
+])
+MU = numpy.stack([
+    numpy.linspace(2.0, 3.2, numpy.prod(SHAPE)).reshape(SHAPE),
+    numpy.linspace(2.1, 3.3, numpy.prod(SHAPE)).reshape(SHAPE),
+    numpy.linspace(2.2, 3.4, numpy.prod(SHAPE)).reshape(SHAPE),
+])
+
+E_FIELD = (numpy.arange(3 * numpy.prod(SHAPE)).reshape((3, *SHAPE)) + 0.5j).astype(complex)
+H_FIELD = (numpy.arange(3 * numpy.prod(SHAPE)).reshape((3, *SHAPE)) * 0.25 - 0.75j).astype(complex)
+
+PEC = numpy.zeros((3, *SHAPE), dtype=float)
+PEC[1, 0, 1, 0] = 1.0
+PMC = numpy.zeros((3, *SHAPE), dtype=float)
+PMC[2, 1, 2, 1] = 1.0
+
+TF_REGION = numpy.zeros((3, *SHAPE), dtype=float)
+TF_REGION[:, 0, 1, 0] = 1.0
+
+BOUNDARY_SHAPE = (3, 4, 3)
+BOUNDARY_DXES = [
+    [numpy.array([1.0, 1.5, 0.8]), numpy.array([0.75, 1.25, 1.5, 0.9]), numpy.array([1.2, 0.8, 1.1])],
+    [numpy.array([0.9, 1.4, 1.0]), numpy.array([0.8, 1.1, 1.4, 1.0]), numpy.array([1.0, 0.7, 1.3])],
+]
+BOUNDARY_EPSILON = numpy.stack([
+    numpy.linspace(1.0, 2.2, numpy.prod(BOUNDARY_SHAPE)).reshape(BOUNDARY_SHAPE),
+    numpy.linspace(1.1, 2.3, numpy.prod(BOUNDARY_SHAPE)).reshape(BOUNDARY_SHAPE),
+    numpy.linspace(1.2, 2.4, numpy.prod(BOUNDARY_SHAPE)).reshape(BOUNDARY_SHAPE),
+])
+BOUNDARY_FIELD = (numpy.arange(3 * numpy.prod(BOUNDARY_SHAPE)).reshape((3, *BOUNDARY_SHAPE)) + 0.5j).astype(complex)
+
+
+def apply_fdfd_matrix(op, field: numpy.ndarray, shape: tuple[int, ...] = SHAPE) -> numpy.ndarray:
+    return unvec(op @ vec(field), shape)

meanas/test/_solver_cases.py

@@ -0,0 +1,56 @@
+import dataclasses
+import numpy
+from scipy import sparse
+import scipy.sparse.linalg as spalg
+
+from ._test_builders import unit_dxes
+
+
+@dataclasses.dataclass(frozen=True)
+class DiagonalEigenCase:
+    operator: sparse.csr_matrix
+    linear_operator: spalg.LinearOperator
+    guess_default: numpy.ndarray
+    guess_sparse: numpy.ndarray
+
+
+@dataclasses.dataclass(frozen=True)
+class SolverPlumbingCase:
+    omega: float
+    a0: sparse.csr_matrix
+    pl: sparse.csr_matrix
+    pr: sparse.csr_matrix
+    j: numpy.ndarray
+    guess: numpy.ndarray
+    solver_result: numpy.ndarray
+    dxes: tuple[tuple[numpy.ndarray, ...], tuple[numpy.ndarray, ...]]
+    epsilon: numpy.ndarray
+
+
+def diagonal_eigen_case() -> DiagonalEigenCase:
+    operator = sparse.diags([5.0, 3.0, 1.0, -2.0]).tocsr()
+    linear_operator = spalg.LinearOperator(
+        shape=operator.shape,
+        dtype=complex,
+        matvec=lambda vv: operator @ vv,
+    )
+    return DiagonalEigenCase(
+        operator=operator,
+        linear_operator=linear_operator,
+        guess_default=numpy.array([0.0, 1.0, 1e-6, 0.0], dtype=complex),
+        guess_sparse=numpy.array([1.0, 0.1, 0.0, 0.0], dtype=complex),
+    )
+
+
+def solver_plumbing_case() -> SolverPlumbingCase:
+    return SolverPlumbingCase(
+        omega=2.0,
+        a0=sparse.csr_matrix(numpy.array([[1.0 + 2.0j, 2.0], [3.0 - 1.0j, 4.0]])),
+        pl=sparse.csr_matrix(numpy.array([[2.0, 0.0], [0.0, 3.0j]])),
+        pr=sparse.csr_matrix(numpy.array([[0.5, 0.0], [0.0, -2.0j]])),
+        j=numpy.array([1.0 + 0.5j, -2.0]),
+        guess=numpy.array([0.25 - 0.75j, 1.5 + 0.5j]),
+        solver_result=numpy.array([3.0 - 1.0j, -4.0 + 2.0j]),
+        dxes=unit_dxes((1, 1, 1)),
+        epsilon=numpy.ones(2),
+    )

meanas/test/_test_builders.py

@@ -0,0 +1,22 @@
+import numpy
+
+
+def real_ramp(shape: tuple[int, ...], *, scale: float = 1.0, offset: float = 0.0) -> numpy.ndarray:
+    return numpy.arange(numpy.prod(shape), dtype=float).reshape(shape, order='C') * scale + offset
+
+
+def complex_ramp(
+        shape: tuple[int, ...],
+        *,
+        scale: float = 1.0,
+        offset: float = 0.0,
+        imag_scale: float = 0.0,
+        imag_offset: float = 0.0,
+        ) -> numpy.ndarray:
+    real = real_ramp(shape, scale=scale, offset=offset)
+    imag = real_ramp(shape, scale=imag_scale, offset=imag_offset)
+    return (real + 1j * imag).astype(complex)
+
+
+def unit_dxes(shape: tuple[int, ...]) -> tuple[tuple[numpy.ndarray, ...], tuple[numpy.ndarray, ...]]:
+    return tuple(tuple(numpy.ones(length) for length in shape) for _ in range(2))  # type: ignore[return-value]

meanas/test/conftest.py

@@ -9,7 +9,7 @@ import numpy
 from numpy.typing import NDArray
 import pytest       # type: ignore
 
-from .utils import PRNG
+from .utils import make_prng
 
 
 FixtureRequest = Any
@@ -42,6 +42,7 @@ def epsilon(
         epsilon_bg: float,
         epsilon_fg: float,
         ) -> NDArray[numpy.float64]:
+    prng = make_prng()
     is3d = (numpy.array(shape) == 1).sum() == 0
     if is3d:
         if request.param == '000':
@@ -57,9 +58,11 @@ def epsilon(
     elif request.param == '000':
         epsilon[:, 0, 0, 0] = epsilon_fg
     elif request.param == 'random':
-        epsilon[:] = PRNG.uniform(low=min(epsilon_bg, epsilon_fg),
+        epsilon[:] = prng.uniform(
-                                  high=max(epsilon_bg, epsilon_fg),
+            low=min(epsilon_bg, epsilon_fg),
-                                  size=shape)
+            high=max(epsilon_bg, epsilon_fg),
+            size=shape,
+        )
 
     return epsilon
 
@@ -80,6 +83,7 @@ def dxes(
         shape: tuple[int, ...],
         dx: float,
         ) -> list[list[NDArray[numpy.float64]]]:
+    prng = make_prng()
     if request.param == 'uniform':
         dxes = [[numpy.full(s, dx) for s in shape[1:]] for _ in range(2)]
     elif request.param == 'centerbig':
@@ -88,8 +92,7 @@ def dxes(
             for ax in (0, 1, 2):
                 dxes[eh][ax][dxes[eh][ax].size // 2] *= 1.1
     elif request.param == 'random':
-        dxe = [PRNG.uniform(low=1.0 * dx, high=1.1 * dx, size=s) for s in shape[1:]]
+        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]
     return dxes
-

meanas/test/test_bloch_foundations.py

@@ -0,0 +1,73 @@
+import numpy
+from numpy.linalg import norm
+
+from ..fdfd import bloch
+from ._bloch_case import EPSILON, G_MATRIX, H_MN, K0_AXIAL, K0_GENERAL, MU, SHAPE, ZERO_H_MN
+from .utils import assert_close
+
+def test_generate_kmn_general_case_returns_orthonormal_basis() -> None:
+    k_mag, m_vecs, n_vecs = bloch.generate_kmn(K0_GENERAL, G_MATRIX, SHAPE)
+
+    assert k_mag.shape == SHAPE + (1,)
+    assert m_vecs.shape == SHAPE + (3,)
+    assert n_vecs.shape == SHAPE + (3,)
+    assert numpy.isfinite(k_mag).all()
+    assert numpy.isfinite(m_vecs).all()
+    assert numpy.isfinite(n_vecs).all()
+
+    assert_close(norm(m_vecs.reshape(-1, 3), axis=1), 1.0)
+    assert_close(norm(n_vecs.reshape(-1, 3), axis=1), 1.0)
+    assert_close(numpy.sum(m_vecs * n_vecs, axis=3), 0.0, atol=1e-12)
+
+
+def test_generate_kmn_z_aligned_uses_default_transverse_basis() -> None:
+    k_mag, m_vecs, n_vecs = bloch.generate_kmn(K0_AXIAL, G_MATRIX, (1, 1, 1))
+
+    assert numpy.isfinite(k_mag).all()
+    assert_close(m_vecs[0, 0, 0], [0.0, 1.0, 0.0])
+    assert_close(numpy.sum(m_vecs * n_vecs, axis=3), 0.0, atol=1e-12)
+    assert_close(norm(n_vecs.reshape(-1, 3), axis=1), 1.0)
+
+
+def test_maxwell_operator_returns_finite_column_vector_without_mu() -> None:
+    operator = bloch.maxwell_operator(K0_GENERAL, G_MATRIX, EPSILON)
+
+    result = operator(H_MN.copy())
+    zero_result = operator(ZERO_H_MN.copy())
+
+    assert result.shape == (2 * numpy.prod(SHAPE), 1)
+    assert numpy.isfinite(result).all()
+    assert_close(zero_result, 0.0)
+
+
+def test_maxwell_operator_returns_finite_column_vector_with_mu() -> None:
+    operator = bloch.maxwell_operator(K0_GENERAL, G_MATRIX, EPSILON, MU)
+
+    result = operator(H_MN.copy())
+    zero_result = operator(ZERO_H_MN.copy())
+
+    assert result.shape == (2 * numpy.prod(SHAPE), 1)
+    assert numpy.isfinite(result).all()
+    assert_close(zero_result, 0.0)
+
+
+def test_inverse_maxwell_operator_returns_finite_column_vector_for_both_mu_branches() -> None:
+    for mu in (None, MU):
+        operator = bloch.inverse_maxwell_operator_approx(K0_GENERAL, G_MATRIX, EPSILON, mu)
+
+        result = operator(H_MN.copy())
+        zero_result = operator(ZERO_H_MN.copy())
+
+        assert result.shape == (2 * numpy.prod(SHAPE), 1)
+        assert numpy.isfinite(result).all()
+        assert_close(zero_result, 0.0)
+
+
+def test_bloch_field_converters_return_finite_fields() -> None:
+    e_field = bloch.hmn_2_exyz(K0_GENERAL, G_MATRIX, EPSILON)(H_MN.copy())
+    h_field = bloch.hmn_2_hxyz(K0_GENERAL, G_MATRIX, EPSILON)(H_MN.copy())
+
+    assert e_field.shape == (3, *SHAPE)
+    assert h_field.shape == (3, *SHAPE)
+    assert numpy.isfinite(e_field).all()
+    assert numpy.isfinite(h_field).all()

meanas/test/test_bloch_interactions.py

@@ -0,0 +1,315 @@
+import numpy
+import pytest
+from numpy.testing import assert_allclose
+from types import SimpleNamespace
+
+from ..fdfd import bloch
+from ._bloch_case import EPSILON, G_MATRIX, H_SIZE, K0_X, SHAPE, Y0, Y0_TWO_MODE, build_overlap_fixture
+from .utils import assert_close
+
+
+def test_rtrace_atb_matches_real_frobenius_inner_product() -> None:
+    a_mat = numpy.array([[1.0 + 2.0j, 3.0 - 1.0j], [2.0j, 4.0]], dtype=complex)
+    b_mat = numpy.array([[5.0 - 1.0j, 1.0 + 1.0j], [2.0, 3.0j]], dtype=complex)
+    expected = numpy.real(numpy.sum(a_mat.conj() * b_mat))
+
+    assert bloch._rtrace_AtB(a_mat, b_mat) == expected
+
+
+def test_symmetrize_returns_hermitian_average() -> None:
+    matrix = numpy.array([[1.0 + 2.0j, 3.0 - 1.0j], [2.0j, 4.0]], dtype=complex)
+    result = bloch._symmetrize(matrix)
+
+    assert_close(result, 0.5 * (matrix + matrix.conj().T))
+    assert_close(result, result.conj().T)
+
+
+def test_inner_product_is_nonmutating_and_obeys_sign_symmetry() -> None:
+    e_in, h_in, e_out, h_out = build_overlap_fixture()
+    originals = (e_in.copy(), h_in.copy(), e_out.copy(), h_out.copy())
+
+    pp = bloch.inner_product(e_out, h_out, e_in, h_in)
+    pn = bloch.inner_product(e_out, h_out, e_in, -h_in)
+    np_term = bloch.inner_product(e_out, -h_out, e_in, h_in)
+    nn = bloch.inner_product(e_out, -h_out, e_in, -h_in)
+
+    assert_close(pp, 0.8164965809277263 + 0.0j)
+    assert_close(pp, -nn, atol=1e-12, rtol=1e-12)
+    assert_close(pn, -np_term, atol=1e-12, rtol=1e-12)
+    assert numpy.array_equal(e_in, originals[0])
+    assert numpy.array_equal(h_in, originals[1])
+    assert numpy.array_equal(e_out, originals[2])
+    assert numpy.array_equal(h_out, originals[3])
+
+
+def test_trq_returns_expected_transmission_and_reflection() -> None:
+    e_in, h_in, e_out, h_out = build_overlap_fixture()
+
+    transmission, reflection = bloch.trq(e_in, h_in, e_out, h_out)
+
+    assert_close(transmission, 0.9797958971132713 + 0.0j, atol=1e-12, rtol=1e-12)
+    assert_close(reflection, 0.2 + 0.0j, atol=1e-12, rtol=1e-12)
+
+
+def test_eigsolve_returns_finite_modes_with_small_residual() -> None:
+    callback_count = 0
+
+    def callback() -> None:
+        nonlocal callback_count
+        callback_count += 1
+
+    eigvals, eigvecs = bloch.eigsolve(
+        1,
+        K0_X,
+        G_MATRIX,
+        EPSILON,
+        tolerance=1e-6,
+        max_iters=50,
+        y0=Y0,
+        callback=callback,
+    )
+
+    operator = bloch.maxwell_operator(K0_X, G_MATRIX, EPSILON)
+    eigvec = eigvecs[0] / numpy.linalg.norm(eigvecs[0])
+    residual = numpy.linalg.norm(operator(eigvec).reshape(-1) - eigvals[0] * eigvec) / numpy.linalg.norm(eigvec)
+
+    assert eigvals.shape == (1,)
+    assert eigvecs.shape == (1, H_SIZE)
+    assert numpy.isfinite(eigvals).all()
+    assert numpy.isfinite(eigvecs).all()
+    assert residual < 1e-5
+    assert callback_count > 0
+
+
+def test_eigsolve_without_initial_guess_returns_finite_modes() -> None:
+    eigvals, eigvecs = bloch.eigsolve(
+        1,
+        K0_X,
+        G_MATRIX,
+        EPSILON,
+        tolerance=1e-6,
+        max_iters=20,
+        y0=None,
+    )
+
+    operator = bloch.maxwell_operator(K0_X, G_MATRIX, EPSILON)
+    eigvec = eigvecs[0] / numpy.linalg.norm(eigvecs[0])
+    residual = numpy.linalg.norm(operator(eigvec).reshape(-1) - eigvals[0] * eigvec) / numpy.linalg.norm(eigvec)
+
+    assert eigvals.shape == (1,)
+    assert eigvecs.shape == (1, H_SIZE)
+    assert numpy.isfinite(eigvals).all()
+    assert numpy.isfinite(eigvecs).all()
+    assert residual < 1e-5
+
+
+def test_eigsolve_recovers_from_singular_initial_subspace(monkeypatch: pytest.MonkeyPatch) -> None:
+    class FakeRng:
+        def __init__(self) -> None:
+            self.calls = 0
+
+        def random(self, shape: tuple[int, ...]) -> numpy.ndarray:
+            self.calls += 1
+            return numpy.arange(numpy.prod(shape), dtype=float).reshape(shape) + self.calls
+
+    fake_rng = FakeRng()
+    singular_y0 = numpy.vstack([Y0_TWO_MODE[0], Y0_TWO_MODE[0]])
+    monkeypatch.setattr(bloch.numpy.random, 'default_rng', lambda: fake_rng)
+
+    eigvals, eigvecs = bloch.eigsolve(
+        2,
+        K0_X,
+        G_MATRIX,
+        EPSILON,
+        tolerance=1e-6,
+        max_iters=20,
+        y0=singular_y0,
+    )
+
+    assert fake_rng.calls == 2
+    assert eigvals.shape == (2,)
+    assert eigvecs.shape == (2, H_SIZE)
+    assert numpy.isfinite(eigvals).all()
+    assert numpy.isfinite(eigvecs).all()
+
+
+def test_eigsolve_reconditions_large_trace_initial_subspace(monkeypatch: pytest.MonkeyPatch) -> None:
+    original_inv = bloch.numpy.linalg.inv
+    original_sqrtm = bloch.scipy.linalg.sqrtm
+    sqrtm_calls = 0
+    inv_calls = 0
+
+    def inv_with_large_first_trace(matrix: numpy.ndarray) -> numpy.ndarray:
+        nonlocal inv_calls
+        inv_calls += 1
+        if inv_calls == 1:
+            return numpy.eye(matrix.shape[0], dtype=complex) * 1e9
+        return original_inv(matrix)
+
+    def sqrtm_wrapper(matrix: numpy.ndarray) -> numpy.ndarray:
+        nonlocal sqrtm_calls
+        sqrtm_calls += 1
+        return original_sqrtm(matrix)
+
+    monkeypatch.setattr(bloch.numpy.linalg, 'inv', inv_with_large_first_trace)
+    monkeypatch.setattr(bloch.scipy.linalg, 'sqrtm', sqrtm_wrapper)
+
+    eigvals, eigvecs = bloch.eigsolve(
+        2,
+        K0_X,
+        G_MATRIX,
+        EPSILON,
+        tolerance=1e-6,
+        max_iters=20,
+        y0=Y0_TWO_MODE,
+    )
+
+    assert sqrtm_calls >= 2
+    assert eigvals.shape == (2,)
+    assert eigvecs.shape == (2, H_SIZE)
+    assert numpy.isfinite(eigvals).all()
+    assert numpy.isfinite(eigvecs).all()
+
+
+def test_eigsolve_qi_memoization_reuses_cached_theta(monkeypatch: pytest.MonkeyPatch) -> None:
+    def fake_minimize_scalar(func, method: str, bounds: tuple[float, float], options: dict[str, float]) -> SimpleNamespace:
+        theta = 0.3
+        first = func(theta)
+        second = func(theta)
+        assert_allclose(second, first)
+        return SimpleNamespace(fun=second, x=theta)
+
+    monkeypatch.setattr(bloch.scipy.optimize, 'minimize_scalar', fake_minimize_scalar)
+
+    eigvals, eigvecs = bloch.eigsolve(
+        1,
+        K0_X,
+        G_MATRIX,
+        EPSILON,
+        tolerance=1e-6,
+        max_iters=1,
+        y0=Y0,
+    )
+
+    assert eigvals.shape == (1,)
+    assert eigvecs.shape == (1, H_SIZE)
+    assert numpy.isfinite(eigvals).all()
+    assert numpy.isfinite(eigvecs).all()
+
+
+@pytest.mark.parametrize('theta', [numpy.pi / 2 - 1e-8, 1e-8])
+def test_eigsolve_qi_taylor_expansions_return_finite_modes(monkeypatch: pytest.MonkeyPatch, theta: float) -> None:
+    original_inv = bloch.numpy.linalg.inv
+    inv_calls = 0
+
+    def inv_raise_once_for_q(matrix: numpy.ndarray) -> numpy.ndarray:
+        nonlocal inv_calls
+        inv_calls += 1
+        if inv_calls == 3:
+            raise numpy.linalg.LinAlgError('forced singular Q')
+        return original_inv(matrix)
+
+    def fake_minimize_scalar(func, method: str, bounds: tuple[float, float], options: dict[str, float]) -> SimpleNamespace:
+        value = func(theta)
+        return SimpleNamespace(fun=value, x=theta)
+
+    monkeypatch.setattr(bloch.numpy.linalg, 'inv', inv_raise_once_for_q)
+    monkeypatch.setattr(bloch.scipy.optimize, 'minimize_scalar', fake_minimize_scalar)
+
+    eigvals, eigvecs = bloch.eigsolve(
+        1,
+        K0_X,
+        G_MATRIX,
+        EPSILON,
+        tolerance=1e-6,
+        max_iters=1,
+        y0=Y0,
+    )
+
+    assert eigvals.shape == (1,)
+    assert eigvecs.shape == (1, H_SIZE)
+    assert numpy.isfinite(eigvals).all()
+    assert numpy.isfinite(eigvecs).all()
+
+
+def test_eigsolve_qi_inexplicable_singularity_raises(monkeypatch: pytest.MonkeyPatch) -> None:
+    original_inv = bloch.numpy.linalg.inv
+    inv_calls = 0
+
+    def inv_raise_once_for_q(matrix: numpy.ndarray) -> numpy.ndarray:
+        nonlocal inv_calls
+        inv_calls += 1
+        if inv_calls == 3:
+            raise numpy.linalg.LinAlgError('forced singular Q')
+        return original_inv(matrix)
+
+    def fake_minimize_scalar(func, method: str, bounds: tuple[float, float], options: dict[str, float]) -> SimpleNamespace:
+        func(numpy.pi / 4)
+        raise AssertionError('unreachable after trace_func exception')
+
+    monkeypatch.setattr(bloch.numpy.linalg, 'inv', inv_raise_once_for_q)
+    monkeypatch.setattr(bloch.scipy.optimize, 'minimize_scalar', fake_minimize_scalar)
+
+    with pytest.raises(Exception, match='Inexplicable singularity in trace_func'):
+        bloch.eigsolve(
+            1,
+            K0_X,
+            G_MATRIX,
+            EPSILON,
+            tolerance=1e-6,
+            max_iters=1,
+            y0=Y0,
+        )
+
+
+def test_find_k_returns_vector_frequency_and_callbacks() -> None:
+    target_eigvals, _target_eigvecs = bloch.eigsolve(
+        1,
+        K0_X,
+        G_MATRIX,
+        EPSILON,
+        tolerance=1e-6,
+        max_iters=50,
+        y0=Y0,
+    )
+    target_frequency = float(numpy.sqrt(abs(numpy.real(target_eigvals[0]))))
+
+    solve_calls = 0
+    iter_calls = 0
+
+    def solve_callback(k_mag: float, eigvals: numpy.ndarray, eigvecs: numpy.ndarray, frequency: float) -> None:
+        nonlocal solve_calls
+        solve_calls += 1
+        assert eigvals.shape == (1,)
+        assert eigvecs.shape == (1, H_SIZE)
+        assert isinstance(k_mag, float)
+        assert isinstance(frequency, float)
+
+    def iter_callback() -> None:
+        nonlocal iter_calls
+        iter_calls += 1
+
+    found_k, found_frequency, eigvals, eigvecs = bloch.find_k(
+        target_frequency,
+        1e-4,
+        [1, 0, 0],
+        G_MATRIX,
+        EPSILON,
+        band=0,
+        k_bounds=(0.05, 0.15),
+        v0=Y0,
+        solve_callback=solve_callback,
+        iter_callback=iter_callback,
+    )
+
+    assert found_k.shape == (3,)
+    assert numpy.isfinite(found_k).all()
+    assert_close(numpy.cross(found_k, [1.0, 0.0, 0.0]), 0.0, atol=1e-12, rtol=1e-12)
+    assert_close(found_k, K0_X, atol=1e-4, rtol=1e-4)
+    assert abs(found_frequency - target_frequency) <= 1e-4
+    assert eigvals.shape == (1,)
+    assert eigvecs.shape == (1, H_SIZE)
+    assert numpy.isfinite(eigvals).all()
+    assert numpy.isfinite(eigvecs).all()
+    assert solve_calls > 0
+    assert iter_calls > 0

meanas/test/test_eigensolvers_numerics.py

@@ -0,0 +1,101 @@
+import numpy
+from numpy.linalg import norm
+import pytest
+
+from ._solver_cases import diagonal_eigen_case
+from .utils import assert_close
+from ..eigensolvers import power_iteration, rayleigh_quotient_iteration, signed_eigensolve
+
+
+def test_rayleigh_quotient_iteration_with_linear_operator() -> None:
+    case = diagonal_eigen_case()
+
+    def dense_solver(
+            shifted_operator,
+            rhs: numpy.ndarray,
+            ) -> numpy.ndarray:
+        basis = numpy.eye(case.operator.shape[0], dtype=complex)
+        columns = [shifted_operator.matvec(basis[:, ii]) for ii in range(case.operator.shape[0])]
+        dense_matrix = numpy.column_stack(columns)
+        return numpy.linalg.lstsq(dense_matrix, rhs, rcond=None)[0]
+
+    eigval, eigvec = rayleigh_quotient_iteration(
+        case.linear_operator,
+        case.guess_default,
+        iterations=8,
+        solver=dense_solver,
+        )
+
+    residual = norm(case.operator @ eigvec - eigval * eigvec)
+    assert abs(eigval - 3.0) < 1e-12
+    assert residual < 1e-12
+
+
+def test_signed_eigensolve_negative_returns_largest_negative_mode() -> None:
+    case = diagonal_eigen_case()
+
+    eigvals, eigvecs = signed_eigensolve(case.operator, how_many=1, negative=True)
+
+    assert eigvals.shape == (1,)
+    assert eigvecs.shape == (4, 1)
+    assert abs(eigvals[0] + 2.0) < 1e-12
+    assert abs(eigvecs[3, 0]) > 0.99
+
+
+def test_rayleigh_quotient_iteration_uses_default_linear_operator_solver() -> None:
+    case = diagonal_eigen_case()
+
+    eigval, eigvec = rayleigh_quotient_iteration(
+        case.linear_operator,
+        case.guess_default,
+        iterations=8,
+    )
+
+    residual = norm(case.operator @ eigvec - eigval * eigvec)
+    assert abs(eigval - 3.0) < 1e-12
+    assert residual < 1e-12
+
+
+def test_signed_eigensolve_linear_operator_fallback_returns_dominant_positive_mode() -> None:
+    case = diagonal_eigen_case()
+
+    eigvals, eigvecs = signed_eigensolve(case.linear_operator, how_many=1)
+
+    assert eigvals.shape == (1,)
+    assert eigvecs.shape == (4, 1)
+    assert_close(eigvals[0], 5.0, atol=1e-12, rtol=1e-12)
+    assert abs(eigvecs[0, 0]) > 0.99
+
+
+def test_power_iteration_finds_dominant_mode() -> None:
+    case = diagonal_eigen_case()
+    eigval, eigvec = power_iteration(case.operator, guess_vector=numpy.ones(4, dtype=complex), iterations=20)
+
+    assert eigval == pytest.approx(5.0, rel=1e-6)
+    assert abs(eigvec[0]) > abs(eigvec[1])
+
+
+def test_rayleigh_quotient_iteration_refines_known_sparse_mode() -> None:
+    case = diagonal_eigen_case()
+
+    def solver(matrix, rhs: numpy.ndarray) -> numpy.ndarray:
+        return numpy.linalg.lstsq(matrix.toarray(), rhs, rcond=None)[0]
+
+    eigval, eigvec = rayleigh_quotient_iteration(
+        case.operator,
+        case.guess_sparse,
+        iterations=8,
+        solver=solver,
+    )
+
+    residual = numpy.linalg.norm(case.operator @ eigvec - eigval * eigvec)
+    assert eigval == pytest.approx(3.0, rel=1e-6)
+    assert residual < 1e-8
+
+
+def test_signed_eigensolve_returns_largest_positive_modes() -> None:
+    case = diagonal_eigen_case()
+    eigvals, eigvecs = signed_eigensolve(case.operator, how_many=2)
+
+    assert_close(eigvals, [3.0, 5.0], atol=1e-6)
+    assert eigvecs.shape == (4, 2)

meanas/test/test_eme_numerics.py

@@ -0,0 +1,132 @@
+import numpy
+from scipy import sparse
+
+from ..fdmath import vec
+from ..fdfd import eme
+from ._test_builders import complex_ramp, unit_dxes
+from .utils import assert_close
+
+
+SHAPE = (3, 2, 2)
+DXES = unit_dxes((2, 2))
+WAVENUMBERS_L = numpy.array([1.0, 0.8])
+WAVENUMBERS_R = numpy.array([0.9, 0.7])
+
+
+def _mode(scale: float) -> tuple[numpy.ndarray, numpy.ndarray]:
+    e_field = complex_ramp(SHAPE, offset=1.0 + scale)
+    h_field = complex_ramp(SHAPE, scale=0.2, offset=2.0, imag_offset=0.05 * scale)
+    return vec(e_field), vec(h_field)
+
+
+def _mode_sets() -> tuple[list[tuple[numpy.ndarray, numpy.ndarray]], list[tuple[numpy.ndarray, numpy.ndarray]]]:
+    left_modes = [_mode(0.0), _mode(0.7)]
+    right_modes = [_mode(1.4), _mode(2.1)]
+    return left_modes, right_modes
+
+
+def _gain_only_tr(*args, **kwargs) -> tuple[numpy.ndarray, numpy.ndarray]:
+    return numpy.array([[2.0, 0.0], [0.0, 0.5]]), numpy.zeros((2, 2))
+
+
+def _gain_and_reflection_tr(*args, **kwargs) -> tuple[numpy.ndarray, numpy.ndarray]:
+    return numpy.array([[2.0, 0.0], [0.0, 0.5]]), numpy.array([[0.0, 1.0], [2.0, 0.0]])
+
+
+def _nonsymmetric_tr(left_marker: object):
+    def fake_get_tr(_eh_left, wavenumbers_left, _eh_right, _wavenumbers_right, **kwargs):
+        if wavenumbers_left is left_marker:
+            return (
+                numpy.array([[1.0, 2.0], [0.5, 1.0]]),
+                numpy.array([[0.0, 1.0], [2.0, 0.0]]),
+            )
+        return (
+            numpy.array([[1.0, -1.0], [0.0, 1.0]]),
+            numpy.array([[0.0, 0.5], [1.5, 0.0]]),
+        )
+
+    return fake_get_tr
+
+
+def test_get_tr_returns_finite_bounded_transfer_matrices() -> None:
+    left_modes, right_modes = _mode_sets()
+
+    transmission, reflection = eme.get_tr(
+        left_modes,
+        WAVENUMBERS_L,
+        right_modes,
+        WAVENUMBERS_R,
+        dxes=DXES,
+    )
+
+    singular_values = numpy.linalg.svd(transmission, compute_uv=False)
+
+    assert transmission.shape == (2, 2)
+    assert reflection.shape == (2, 2)
+    assert numpy.isfinite(transmission).all()
+    assert numpy.isfinite(reflection).all()
+    assert (singular_values <= 1.0 + 1e-12).all()
+
+
+def test_get_abcd_matches_explicit_block_formula() -> None:
+    left_modes, right_modes = _mode_sets()
+    t12, r12 = eme.get_tr(left_modes, WAVENUMBERS_L, right_modes, WAVENUMBERS_R, dxes=DXES)
+    t21, r21 = eme.get_tr(right_modes, WAVENUMBERS_R, left_modes, WAVENUMBERS_L, dxes=DXES)
+    t21_inv = numpy.linalg.pinv(t21)
+
+    expected = numpy.block([
+        [t12 - r21 @ t21_inv @ r12, r21 @ t21_inv],
+        [-t21_inv @ r12, t21_inv],
+    ])
+    abcd = eme.get_abcd(left_modes, WAVENUMBERS_L, right_modes, WAVENUMBERS_R, dxes=DXES)
+
+    assert sparse.issparse(abcd)
+    assert abcd.shape == (4, 4)
+    assert_close(abcd.toarray(), expected)
+
+
+def test_get_s_plain_matches_block_assembly_from_get_tr() -> None:
+    left_modes, right_modes = _mode_sets()
+    t12, r12 = eme.get_tr(left_modes, WAVENUMBERS_L, right_modes, WAVENUMBERS_R, dxes=DXES)
+    t21, r21 = eme.get_tr(right_modes, WAVENUMBERS_R, left_modes, WAVENUMBERS_L, dxes=DXES)
+    expected = numpy.block([[r12, t12], [t21, r21]])
+
+    ss = eme.get_s(left_modes, WAVENUMBERS_L, right_modes, WAVENUMBERS_R, dxes=DXES)
+
+    assert ss.shape == (4, 4)
+    assert numpy.isfinite(ss).all()
+    assert_close(ss, expected)
+
+
+def test_get_s_force_nogain_caps_singular_values(monkeypatch) -> None:
+    monkeypatch.setattr(eme, 'get_tr', _gain_only_tr)
+
+    plain_s = eme.get_s(None, None, None, None)
+    clipped_s = eme.get_s(None, None, None, None, force_nogain=True)
+
+    plain_singular_values = numpy.linalg.svd(plain_s, compute_uv=False)
+    clipped_singular_values = numpy.linalg.svd(clipped_s, compute_uv=False)
+
+    assert plain_singular_values.max() > 1.0
+    assert (clipped_singular_values <= 1.0 + 1e-12).all()
+    assert numpy.isfinite(clipped_s).all()
+
+
+def test_get_s_force_reciprocal_symmetrizes_output(monkeypatch) -> None:
+    left = object()
+    right = object()
+
+    monkeypatch.setattr(eme, 'get_tr', _nonsymmetric_tr(left))
+    ss = eme.get_s(None, left, None, right, force_reciprocal=True)
+
+    assert_close(ss, ss.T)
+
+
+def test_get_s_force_nogain_and_reciprocal_returns_finite_output(monkeypatch) -> None:
+    monkeypatch.setattr(eme, 'get_tr', _gain_and_reflection_tr)
+    ss = eme.get_s(None, None, None, None, force_nogain=True, force_reciprocal=True)
+
+    assert ss.shape == (4, 4)
+    assert numpy.isfinite(ss).all()
+    assert_close(ss, ss.T)
+    assert (numpy.linalg.svd(ss, compute_uv=False) <= 1.0 + 1e-12).all()

meanas/test/test_fdfd_algebra_helpers.py

@@ -0,0 +1,213 @@
+import numpy
+
+from ..fdmath import vec, unvec
+from ..fdmath import functional as fd_functional
+from ..fdfd import operators, scpml
+from ._fdfd_case import (
+    BOUNDARY_DXES,
+    BOUNDARY_EPSILON,
+    BOUNDARY_FIELD,
+    BOUNDARY_SHAPE,
+    DXES,
+    EPSILON,
+    E_FIELD,
+    MU,
+    H_FIELD,
+    OMEGA,
+    PEC,
+    PMC,
+    SHAPE,
+    apply_fdfd_matrix,
+)
+from .utils import assert_close, assert_fields_close
+
+
+def _dense_e_full(mu: numpy.ndarray | None) -> numpy.ndarray:
+    ce = fd_functional.curl_forward(DXES[0])
+    ch = fd_functional.curl_back(DXES[1])
+    pe = numpy.where(PEC, 0.0, 1.0)
+    pm = numpy.where(PMC, 0.0, 1.0)
+
+    masked_e = pe * E_FIELD
+    curl_term = ce(masked_e)
+    if mu is not None:
+        curl_term = curl_term / mu
+    curl_term = pm * curl_term
+    curl_term = ch(curl_term)
+    return pe * (curl_term - OMEGA**2 * EPSILON * masked_e)
+
+
+def _dense_h_full(mu: numpy.ndarray | None) -> numpy.ndarray:
+    ce = fd_functional.curl_forward(DXES[0])
+    ch = fd_functional.curl_back(DXES[1])
+    pe = numpy.where(PEC, 0.0, 1.0)
+    pm = numpy.where(PMC, 0.0, 1.0)
+    magnetic = numpy.ones_like(EPSILON) if mu is None else mu
+
+    masked_h = pm * H_FIELD
+    curl_term = ch(masked_h)
+    curl_term = pe * (curl_term / EPSILON)
+    curl_term = ce(curl_term)
+    return pm * (curl_term - OMEGA**2 * magnetic * masked_h)
+
+
+def _normalized_distance(u: numpy.ndarray, size: int, thickness: int) -> numpy.ndarray:
+    return ((thickness - u).clip(0) + (u - (size - thickness)).clip(0)) / thickness
+
+
+def test_h_full_matches_dense_reference_with_and_without_mu() -> None:
+    for mu in (None, MU):
+        matrix_result = apply_fdfd_matrix(
+            operators.h_full(OMEGA, DXES, vec(EPSILON), None if mu is None else vec(mu), vec(PEC), vec(PMC)),
+            H_FIELD,
+        )
+        dense_result = _dense_h_full(mu)
+        assert_fields_close(matrix_result, dense_result, atol=1e-10, rtol=1e-10)
+
+
+def test_e_full_matches_dense_reference_with_masks() -> None:
+    for mu in (None, MU):
+        matrix_result = apply_fdfd_matrix(
+            operators.e_full(OMEGA, DXES, vec(EPSILON), None if mu is None else vec(mu), vec(PEC), vec(PMC)),
+            E_FIELD,
+        )
+        dense_result = _dense_e_full(mu)
+        assert_fields_close(matrix_result, dense_result, atol=1e-10, rtol=1e-10)
+
+
+def test_h_full_without_masks_matches_dense_reference() -> None:
+    ce = fd_functional.curl_forward(DXES[0])
+    ch = fd_functional.curl_back(DXES[1])
+    dense_result = ce(ch(H_FIELD) / EPSILON) - OMEGA**2 * MU * H_FIELD
+    matrix_result = apply_fdfd_matrix(
+        operators.h_full(OMEGA, DXES, vec(EPSILON), vec(MU)),
+        H_FIELD,
+    )
+    assert_fields_close(matrix_result, dense_result, atol=1e-10, rtol=1e-10)
+
+
+def test_eh_full_matches_manual_block_operator_with_masks() -> None:
+    pe = numpy.where(PEC, 0.0, 1.0)
+    pm = numpy.where(PMC, 0.0, 1.0)
+    ce = fd_functional.curl_forward(DXES[0])
+    ch = fd_functional.curl_back(DXES[1])
+
+    matrix_result = operators.eh_full(OMEGA, DXES, vec(EPSILON), vec(MU), vec(PEC), vec(PMC)) @ numpy.concatenate(
+        [vec(E_FIELD), vec(H_FIELD)],
+    )
+    matrix_e, matrix_h = (unvec(part, SHAPE) for part in numpy.split(matrix_result, 2))
+
+    dense_e = pe * ch(pm * H_FIELD) - pe * (1j * OMEGA * EPSILON * (pe * E_FIELD))
+    dense_h = pm * ce(pe * E_FIELD) + pm * (1j * OMEGA * MU * (pm * H_FIELD))
+
+    assert_fields_close(matrix_e, dense_e, atol=1e-10, rtol=1e-10)
+    assert_fields_close(matrix_h, dense_h, atol=1e-10, rtol=1e-10)
+
+
+def test_e2h_pmc_mask_matches_masked_unmasked_result() -> None:
+    pmc_complement = numpy.where(PMC, 0.0, 1.0)
+    unmasked = apply_fdfd_matrix(operators.e2h(OMEGA, DXES, vec(MU)), E_FIELD)
+    masked = apply_fdfd_matrix(operators.e2h(OMEGA, DXES, vec(MU), vec(PMC)), E_FIELD)
+
+    assert_fields_close(masked, pmc_complement * unmasked, atol=1e-10, rtol=1e-10)
+
+
+def test_poynting_h_cross_matches_negative_e_cross_relation() -> None:
+    h_cross_e = apply_fdfd_matrix(operators.poynting_h_cross(vec(H_FIELD), DXES), E_FIELD)
+    e_cross_h = apply_fdfd_matrix(operators.poynting_e_cross(vec(E_FIELD), DXES), H_FIELD)
+
+    assert_fields_close(h_cross_e, -e_cross_h, atol=1e-10, rtol=1e-10)
+
+
+def test_e_boundary_source_interior_mask_is_independent_of_periodic_edges() -> None:
+    mask = numpy.zeros((3, *BOUNDARY_SHAPE), dtype=float)
+    mask[:, 1, 1, 1] = 1.0
+
+    periodic = operators.e_boundary_source(vec(mask), OMEGA, BOUNDARY_DXES, vec(BOUNDARY_EPSILON), periodic_mask_edges=True)
+    mirrored = operators.e_boundary_source(vec(mask), OMEGA, BOUNDARY_DXES, vec(BOUNDARY_EPSILON), periodic_mask_edges=False)
+
+    assert_close(periodic.toarray(), mirrored.toarray())
+
+
+def test_e_boundary_source_periodic_edges_add_opposite_face_response() -> None:
+    mask = numpy.zeros((3, *BOUNDARY_SHAPE), dtype=float)
+    mask[:, 0, 1, 1] = 1.0
+
+    periodic = operators.e_boundary_source(vec(mask), OMEGA, BOUNDARY_DXES, vec(BOUNDARY_EPSILON), periodic_mask_edges=True)
+    mirrored = operators.e_boundary_source(vec(mask), OMEGA, BOUNDARY_DXES, vec(BOUNDARY_EPSILON), periodic_mask_edges=False)
+    diff = unvec((periodic - mirrored) @ vec(BOUNDARY_FIELD), BOUNDARY_SHAPE)
+
+    assert numpy.isfinite(diff).all()
+    assert_close(diff[:, 1:-1, :, :], 0.0)
+    assert numpy.linalg.norm(diff[:, -1, :, :]) > 0
+
+
+def test_prepare_s_function_matches_closed_form_polynomial() -> None:
+    ln_r = -12.0
+    order = 3.0
+    distances = numpy.array([0.0, 0.25, 0.5, 1.0])
+    s_function = scpml.prepare_s_function(ln_R=ln_r, m=order)
+    expected = (order + 1) * ln_r / 2 * distances**order
+
+    assert_close(s_function(distances), expected)
+
+
+def test_uniform_grid_scpml_matches_expected_stretch_profile() -> None:
+    s_function = scpml.prepare_s_function(ln_R=-12.0, m=3.0)
+    dxes = scpml.uniform_grid_scpml((6, 4, 3), (2, 0, 1), omega=2.0, epsilon_effective=4.0, s_function=s_function)
+    correction = numpy.sqrt(4.0) * 2.0
+
+    for axis, size, thickness in ((0, 6, 2), (2, 3, 1)):
+        grid = numpy.arange(size, dtype=float)
+        expected_a = 1 + 1j * s_function(_normalized_distance(grid, size, thickness)) / correction
+        expected_b = 1 + 1j * s_function(_normalized_distance(grid + 0.5, size, thickness)) / correction
+        assert_close(dxes[0][axis], expected_a)
+        assert_close(dxes[1][axis], expected_b)
+
+    assert_close(dxes[0][1], 1.0)
+    assert_close(dxes[1][1], 1.0)
+    assert numpy.isfinite(dxes[0][0]).all()
+    assert numpy.isfinite(dxes[1][0]).all()
+
+
+def test_uniform_grid_scpml_default_s_function_matches_explicit_default() -> None:
+    implicit = scpml.uniform_grid_scpml((6, 4, 3), (2, 0, 1), omega=2.0)
+    explicit = scpml.uniform_grid_scpml((6, 4, 3), (2, 0, 1), omega=2.0, s_function=scpml.prepare_s_function())
+
+    for implicit_group, explicit_group in zip(implicit, explicit, strict=True):
+        for implicit_axis, explicit_axis in zip(implicit_group, explicit_group, strict=True):
+            assert_close(implicit_axis, explicit_axis)
+
+
+def test_stretch_with_scpml_only_modifies_requested_front_edge() -> None:
+    s_function = scpml.prepare_s_function(ln_R=-12.0, m=3.0)
+    base = [[numpy.ones(6), numpy.ones(4), numpy.ones(3)] for _ in range(2)]
+    stretched = scpml.stretch_with_scpml(base, axis=0, polarity=1, omega=2.0, epsilon_effective=4.0, thickness=2, s_function=s_function)
+
+    assert_close(stretched[0][0][2:], 1.0)
+    assert_close(stretched[1][0][2:], 1.0)
+    assert_close(stretched[0][0][-2:], 1.0)
+    assert_close(stretched[1][0][-2:], 1.0)
+    assert numpy.linalg.norm(stretched[0][0][:2] - 1.0) > 0
+    assert numpy.linalg.norm(stretched[1][0][:2] - 1.0) > 0
+
+
+def test_stretch_with_scpml_only_modifies_requested_back_edge() -> None:
+    s_function = scpml.prepare_s_function(ln_R=-12.0, m=3.0)
+    base = [[numpy.ones(6), numpy.ones(4), numpy.ones(3)] for _ in range(2)]
+    stretched = scpml.stretch_with_scpml(base, axis=0, polarity=-1, omega=2.0, epsilon_effective=4.0, thickness=2, s_function=s_function)
+
+    assert_close(stretched[0][0][:4], 1.0)
+    assert_close(stretched[1][0][:4], 1.0)
+    assert numpy.linalg.norm(stretched[0][0][-2:] - 1.0) > 0
+    assert numpy.linalg.norm(stretched[1][0][-2:] - 1.0) > 0
+
+
+def test_stretch_with_scpml_thickness_zero_is_noop() -> None:
+    s_function = scpml.prepare_s_function(ln_R=-12.0, m=3.0)
+    base = [[numpy.ones(6), numpy.ones(4), numpy.ones(3)] for _ in range(2)]
+    stretched = scpml.stretch_with_scpml(base, axis=0, polarity=-1, omega=2.0, epsilon_effective=4.0, thickness=0, s_function=s_function)
+
+    for grid_group in stretched:
+        for axis_grid in grid_group:
+            assert_close(axis_grid, 1.0)

meanas/test/test_fdfd_farfield.py

@@ -0,0 +1,100 @@
+import numpy
+import pytest
+
+from ..fdfd import farfield
+
+
+NEAR_SHAPE = (2, 3)
+E_NEAR = [numpy.zeros(NEAR_SHAPE, dtype=complex), numpy.zeros(NEAR_SHAPE, dtype=complex)]
+H_NEAR = [numpy.zeros(NEAR_SHAPE, dtype=complex), numpy.zeros(NEAR_SHAPE, dtype=complex)]
+
+
+def test_near_to_farfield_rejects_wrong_length_inputs() -> None:
+    with pytest.raises(Exception, match='E_near must be a length-2 list'):
+        farfield.near_to_farfield(E_NEAR[:1], H_NEAR, dx=0.2, dy=0.3)
+
+    with pytest.raises(Exception, match='H_near must be a length-2 list'):
+        farfield.near_to_farfield(E_NEAR, H_NEAR[:1], dx=0.2, dy=0.3)
+
+
+def test_near_to_farfield_rejects_mismatched_shapes() -> None:
+    bad_h_near = [H_NEAR[0], numpy.zeros((2, 4), dtype=complex)]
+
+    with pytest.raises(Exception, match='All fields must be the same shape'):
+        farfield.near_to_farfield(E_NEAR, bad_h_near, dx=0.2, dy=0.3)
+
+
+def test_near_to_farfield_uses_default_and_scalar_padding_shapes() -> None:
+    default_result = farfield.near_to_farfield(E_NEAR, H_NEAR, dx=0.2, dy=0.3)
+    scalar_result = farfield.near_to_farfield(E_NEAR, H_NEAR, dx=0.2, dy=0.3, padded_size=8)
+
+    assert default_result['E'][0].shape == (2, 4)
+    assert default_result['H'][0].shape == (2, 4)
+    assert scalar_result['E'][0].shape == (8, 8)
+    assert scalar_result['H'][0].shape == (8, 8)
+
+
+def test_far_to_nearfield_rejects_wrong_length_inputs() -> None:
+    ff = farfield.near_to_farfield(E_NEAR, H_NEAR, dx=0.2, dy=0.3, padded_size=8)
+
+    with pytest.raises(Exception, match='E_far must be a length-2 list'):
+        farfield.far_to_nearfield(ff['E'][:1], ff['H'], ff['dkx'], ff['dky'])
+
+    with pytest.raises(Exception, match='H_far must be a length-2 list'):
+        farfield.far_to_nearfield(ff['E'], ff['H'][:1], ff['dkx'], ff['dky'])
+
+
+def test_far_to_nearfield_rejects_mismatched_shapes() -> None:
+    ff = farfield.near_to_farfield(E_NEAR, H_NEAR, dx=0.2, dy=0.3, padded_size=8)
+    bad_h_far = [ff['H'][0], numpy.zeros((8, 4), dtype=complex)]
+
+    with pytest.raises(Exception, match='All fields must be the same shape'):
+        farfield.far_to_nearfield(ff['E'], bad_h_far, ff['dkx'], ff['dky'])
+
+
+def test_far_to_nearfield_uses_default_and_scalar_padding_shapes() -> None:
+    ff = farfield.near_to_farfield(E_NEAR, H_NEAR, dx=0.2, dy=0.3, padded_size=8)
+    default_result = farfield.far_to_nearfield(
+        [field.copy() for field in ff['E']],
+        [field.copy() for field in ff['H']],
+        ff['dkx'],
+        ff['dky'],
+    )
+    scalar_result = farfield.far_to_nearfield(
+        [field.copy() for field in ff['E']],
+        [field.copy() for field in ff['H']],
+        ff['dkx'],
+        ff['dky'],
+        padded_size=4,
+    )
+
+    assert default_result['E'][0].shape == (8, 8)
+    assert default_result['H'][0].shape == (8, 8)
+    assert scalar_result['E'][0].shape == (4, 4)
+    assert scalar_result['H'][0].shape == (4, 4)
+
+
+def test_farfield_roundtrip_supports_rectangular_arrays() -> None:
+    e_near = [numpy.zeros((4, 8), dtype=complex), numpy.zeros((4, 8), dtype=complex)]
+    h_near = [numpy.zeros((4, 8), dtype=complex), numpy.zeros((4, 8), dtype=complex)]
+    e_near[0][1, 3] = 1.0 + 0.25j
+    h_near[1][2, 5] = -0.5j
+
+    ff = farfield.near_to_farfield(e_near, h_near, dx=0.2, dy=0.3, padded_size=(4, 8))
+    restored = farfield.far_to_nearfield(
+        [field.copy() for field in ff['E']],
+        [field.copy() for field in ff['H']],
+        ff['dkx'],
+        ff['dky'],
+        padded_size=(4, 8),
+    )
+
+    assert isinstance(ff['dkx'], float)
+    assert isinstance(ff['dky'], float)
+    assert ff['E'][0].shape == (4, 8)
+    assert restored['E'][0].shape == (4, 8)
+    assert restored['H'][0].shape == (4, 8)
+    assert restored['dx'] == pytest.approx(0.2)
+    assert restored['dy'] == pytest.approx(0.3)
+    assert numpy.isfinite(restored['E'][0]).all()
+    assert numpy.isfinite(restored['H'][0]).all()

meanas/test/test_fdfd_functional.py

@@ -0,0 +1,122 @@
+import numpy
+
+from ..fdmath import unvec, vec
+from ..fdfd import functional, operators
+from ._fdfd_case import DXES, EPSILON, E_FIELD, H_FIELD, MU, OMEGA, SHAPE, TF_REGION, apply_fdfd_matrix
+from .utils import assert_fields_close
+
+
+ATOL = 1e-9
+RTOL = 1e-9
+
+
+def assert_fields_match(actual: numpy.ndarray, expected: numpy.ndarray) -> None:
+    assert_fields_close(actual, expected, atol=ATOL, rtol=RTOL)
+
+
+def test_e_full_matches_sparse_operator_without_mu() -> None:
+    matrix_result = apply_fdfd_matrix(
+        operators.e_full(OMEGA, DXES, vec(EPSILON)),
+        E_FIELD,
+    )
+    functional_result = functional.e_full(OMEGA, DXES, EPSILON)(E_FIELD)
+
+    assert_fields_match(functional_result, matrix_result)
+
+
+def test_e_full_matches_sparse_operator_with_mu() -> None:
+    matrix_result = apply_fdfd_matrix(
+        operators.e_full(OMEGA, DXES, vec(EPSILON), vec(MU)),
+        E_FIELD,
+    )
+    functional_result = functional.e_full(OMEGA, DXES, EPSILON, MU)(E_FIELD)
+
+    assert_fields_match(functional_result, matrix_result)
+
+
+def test_eh_full_matches_sparse_operator_with_mu() -> None:
+    matrix_result = operators.eh_full(OMEGA, DXES, vec(EPSILON), vec(MU)) @ numpy.concatenate([vec(E_FIELD), vec(H_FIELD)])
+    matrix_e, matrix_h = (unvec(part, SHAPE) for part in numpy.split(matrix_result, 2))
+    functional_e, functional_h = functional.eh_full(OMEGA, DXES, EPSILON, MU)(E_FIELD, H_FIELD)
+
+    assert_fields_match(functional_e, matrix_e)
+    assert_fields_match(functional_h, matrix_h)
+
+
+def test_eh_full_matches_sparse_operator_without_mu() -> None:
+    matrix_result = operators.eh_full(OMEGA, DXES, vec(EPSILON)) @ numpy.concatenate([vec(E_FIELD), vec(H_FIELD)])
+    matrix_e, matrix_h = (unvec(part, SHAPE) for part in numpy.split(matrix_result, 2))
+    functional_e, functional_h = functional.eh_full(OMEGA, DXES, EPSILON)(E_FIELD, H_FIELD)
+
+    assert_fields_match(functional_e, matrix_e)
+    assert_fields_match(functional_h, matrix_h)
+
+
+def test_e2h_matches_sparse_operator_with_mu() -> None:
+    matrix_result = apply_fdfd_matrix(
+        operators.e2h(OMEGA, DXES, vec(MU)),
+        E_FIELD,
+    )
+    functional_result = functional.e2h(OMEGA, DXES, MU)(E_FIELD)
+
+    assert_fields_match(functional_result, matrix_result)
+
+
+def test_e2h_matches_sparse_operator_without_mu() -> None:
+    matrix_result = apply_fdfd_matrix(
+        operators.e2h(OMEGA, DXES),
+        E_FIELD,
+    )
+    functional_result = functional.e2h(OMEGA, DXES)(E_FIELD)
+
+    assert_fields_match(functional_result, matrix_result)
+
+
+def test_m2j_matches_sparse_operator_without_mu() -> None:
+    matrix_result = apply_fdfd_matrix(
+        operators.m2j(OMEGA, DXES),
+        H_FIELD,
+    )
+    functional_result = functional.m2j(OMEGA, DXES)(H_FIELD)
+
+    assert_fields_match(functional_result, matrix_result)
+
+
+def test_m2j_matches_sparse_operator_with_mu() -> None:
+    matrix_result = apply_fdfd_matrix(
+        operators.m2j(OMEGA, DXES, vec(MU)),
+        H_FIELD,
+    )
+    functional_result = functional.m2j(OMEGA, DXES, MU)(H_FIELD)
+
+    assert_fields_match(functional_result, matrix_result)
+
+
+def test_e_tfsf_source_matches_sparse_operator_without_mu() -> None:
+    matrix_result = apply_fdfd_matrix(
+        operators.e_tfsf_source(vec(TF_REGION), OMEGA, DXES, vec(EPSILON)),
+        E_FIELD,
+    )
+    functional_result = functional.e_tfsf_source(TF_REGION, OMEGA, DXES, EPSILON)(E_FIELD)
+
+    assert_fields_match(functional_result, matrix_result)
+
+
+def test_e_tfsf_source_matches_sparse_operator_with_mu() -> None:
+    matrix_result = apply_fdfd_matrix(
+        operators.e_tfsf_source(vec(TF_REGION), OMEGA, DXES, vec(EPSILON), vec(MU)),
+        E_FIELD,
+    )
+    functional_result = functional.e_tfsf_source(TF_REGION, OMEGA, DXES, EPSILON, MU)(E_FIELD)
+
+    assert_fields_match(functional_result, matrix_result)
+
+
+def test_poynting_e_cross_h_matches_sparse_operator() -> None:
+    matrix_result = apply_fdfd_matrix(
+        operators.poynting_e_cross(vec(E_FIELD), DXES),
+        H_FIELD,
+    )
+    functional_result = functional.poynting_e_cross_h(DXES)(E_FIELD, H_FIELD)
+
+    assert_fields_match(functional_result, matrix_result)

meanas/test/test_fdfd_solvers.py

@@ -0,0 +1,126 @@
+import numpy
+
+from ..fdfd import solvers
+from ._solver_cases import solver_plumbing_case
+from .utils import assert_close
+
+
+def test_scipy_qmr_wraps_user_callback_without_recursion(monkeypatch) -> None:
+    seen: list[tuple[float, ...]] = []
+
+    def fake_qmr(a, b: numpy.ndarray, **kwargs):
+        kwargs['callback'](numpy.array([1.0, 2.0]))
+        return numpy.array([3.0, 4.0]), 0
+
+    monkeypatch.setattr(solvers.scipy.sparse.linalg, 'qmr', fake_qmr)
+    result = solvers._scipy_qmr(
+        solver_plumbing_case().a0,
+        numpy.array([1.0, 0.0]),
+        callback=lambda xk: seen.append(tuple(xk)),
+    )
+
+    assert_close(result, [3.0, 4.0])
+    assert seen == [(1.0, 2.0)]
+
+
+def test_scipy_qmr_installs_logging_callback_when_missing(monkeypatch) -> None:
+    callback_seen: list[numpy.ndarray] = []
+
+    def fake_qmr(a, b: numpy.ndarray, **kwargs):
+        callback = kwargs['callback']
+        callback(numpy.array([5.0, 6.0]))
+        callback_seen.append(b.copy())
+        return numpy.array([7.0, 8.0]), 0
+
+    monkeypatch.setattr(solvers.scipy.sparse.linalg, 'qmr', fake_qmr)
+    result = solvers._scipy_qmr(solver_plumbing_case().a0, numpy.array([1.0, 0.0]))
+
+    assert_close(result, [7.0, 8.0])
+    assert len(callback_seen) == 1
+
+
+def test_generic_forward_preconditions_system_and_guess(monkeypatch) -> None:
+    case = solver_plumbing_case()
+    captured: dict[str, object] = {}
+
+    monkeypatch.setattr(solvers.operators, 'e_full', lambda *args, **kwargs: case.a0)
+    monkeypatch.setattr(solvers.operators, 'e_full_preconditioners', lambda dxes: (case.pl, case.pr))
+
+    def fake_solver(a, b: numpy.ndarray, **kwargs):
+        captured['a'] = a
+        captured['b'] = b
+        captured['x0'] = kwargs['x0']
+        captured['atol'] = kwargs['atol']
+        return case.solver_result
+
+    result = solvers.generic(
+        omega=case.omega,
+        dxes=case.dxes,
+        J=case.j,
+        epsilon=case.epsilon,
+        matrix_solver=fake_solver,
+        matrix_solver_opts={'atol': 1e-12},
+        E_guess=case.guess,
+    )
+
+    assert_close(captured['a'].toarray(), (case.pl @ case.a0 @ case.pr).toarray())
+    assert_close(captured['b'], case.pl @ (-1j * case.omega * case.j))
+    assert_close(captured['x0'], case.pl @ case.guess)
+    assert captured['atol'] == 1e-12
+    assert_close(result, case.pr @ case.solver_result)
+
+
+def test_generic_adjoint_preconditions_system_and_guess(monkeypatch) -> None:
+    case = solver_plumbing_case()
+    captured: dict[str, object] = {}
+
+    monkeypatch.setattr(solvers.operators, 'e_full', lambda *args, **kwargs: case.a0)
+    monkeypatch.setattr(solvers.operators, 'e_full_preconditioners', lambda dxes: (case.pl, case.pr))
+
+    def fake_solver(a, b: numpy.ndarray, **kwargs):
+        captured['a'] = a
+        captured['b'] = b
+        captured['x0'] = kwargs['x0']
+        captured['rtol'] = kwargs['rtol']
+        return case.solver_result
+
+    result = solvers.generic(
+        omega=case.omega,
+        dxes=case.dxes,
+        J=case.j,
+        epsilon=case.epsilon,
+        matrix_solver=fake_solver,
+        matrix_solver_opts={'rtol': 1e-9},
+        E_guess=case.guess,
+        adjoint=True,
+    )
+
+    expected_matrix = (case.pl @ case.a0 @ case.pr).T.conjugate()
+    assert_close(captured['a'].toarray(), expected_matrix.toarray())
+    assert_close(captured['b'], case.pr.T.conjugate() @ (-1j * case.omega * case.j))
+    assert_close(captured['x0'], case.pr.T.conjugate() @ case.guess)
+    assert captured['rtol'] == 1e-9
+    assert_close(result, case.pl.T.conjugate() @ case.solver_result)
+
+
+def test_generic_without_guess_does_not_inject_x0(monkeypatch) -> None:
+    case = solver_plumbing_case()
+    captured: dict[str, object] = {}
+
+    monkeypatch.setattr(solvers.operators, 'e_full', lambda *args, **kwargs: case.a0)
+    monkeypatch.setattr(solvers.operators, 'e_full_preconditioners', lambda dxes: (case.pl, case.pr))
+
+    def fake_solver(a, b: numpy.ndarray, **kwargs):
+        captured['kwargs'] = kwargs
+        return numpy.array([1.0, -1.0])
+
+    result = solvers.generic(
+        omega=1.0,
+        dxes=case.dxes,
+        J=numpy.array([2.0, 3.0]),
+        epsilon=case.epsilon,
+        matrix_solver=fake_solver,
+    )
+
+    assert 'x0' not in captured['kwargs']
+    assert_close(result, case.pr @ numpy.array([1.0, -1.0]))

meanas/test/test_fdmath_functional.py

@@ -0,0 +1,60 @@
+import numpy
+
+from ..fdmath import functional as fd_functional
+from ..fdmath import operators as fd_operators
+from ..fdmath import vec, unvec
+from .utils import assert_close, assert_fields_close
+
+
+SHAPE = (2, 3, 2)
+DX_E = [numpy.array([1.0, 1.5]), numpy.array([0.75, 1.25, 1.5]), numpy.array([1.2, 0.8])]
+DX_H = [numpy.array([0.9, 1.4]), numpy.array([0.8, 1.1, 1.4]), numpy.array([1.0, 0.7])]
+
+SCALAR_FIELD = (
+    numpy.arange(numpy.prod(SHAPE)).reshape(SHAPE)
+    + 0.1j * numpy.arange(numpy.prod(SHAPE)).reshape(SHAPE)
+).astype(complex)
+VECTOR_FIELD = (numpy.arange(3 * numpy.prod(SHAPE)).reshape((3, *SHAPE)) + 0.25j).astype(complex)
+
+
+def test_deriv_forward_without_dx_matches_numpy_roll() -> None:
+    for axis, deriv in enumerate(fd_functional.deriv_forward()):
+        expected = numpy.roll(SCALAR_FIELD, -1, axis=axis) - SCALAR_FIELD
+        assert_close(deriv(SCALAR_FIELD), expected)
+
+
+def test_deriv_back_without_dx_matches_numpy_roll() -> None:
+    for axis, deriv in enumerate(fd_functional.deriv_back()):
+        expected = SCALAR_FIELD - numpy.roll(SCALAR_FIELD, 1, axis=axis)
+        assert_close(deriv(SCALAR_FIELD), expected)
+
+
+def test_curl_parts_sum_to_full_curl() -> None:
+    curl_forward = fd_functional.curl_forward(DX_E)(VECTOR_FIELD)
+    curl_back = fd_functional.curl_back(DX_H)(VECTOR_FIELD)
+    forward_parts = fd_functional.curl_forward_parts(DX_E)(VECTOR_FIELD)
+    back_parts = fd_functional.curl_back_parts(DX_H)(VECTOR_FIELD)
+
+    for axis in range(3):
+        assert_close(forward_parts[axis][0] + forward_parts[axis][1], curl_forward[axis])
+        assert_close(back_parts[axis][0] + back_parts[axis][1], curl_back[axis])
+
+
+def test_derivatives_match_sparse_operators_on_nonuniform_grid() -> None:
+    for axis, deriv in enumerate(fd_functional.deriv_forward(DX_E)):
+        matrix_result = (fd_operators.deriv_forward(DX_E)[axis] @ SCALAR_FIELD.ravel(order='C')).reshape(SHAPE, order='C')
+        assert_close(deriv(SCALAR_FIELD), matrix_result, atol=1e-12, rtol=1e-12)
+
+    for axis, deriv in enumerate(fd_functional.deriv_back(DX_H)):
+        matrix_result = (fd_operators.deriv_back(DX_H)[axis] @ SCALAR_FIELD.ravel(order='C')).reshape(SHAPE, order='C')
+        assert_close(deriv(SCALAR_FIELD), matrix_result, atol=1e-12, rtol=1e-12)
+
+
+def test_curls_match_sparse_operators_on_nonuniform_grid() -> None:
+    curl_forward = fd_functional.curl_forward(DX_E)(VECTOR_FIELD)
+    curl_back = fd_functional.curl_back(DX_H)(VECTOR_FIELD)
+    matrix_forward = unvec(fd_operators.curl_forward(DX_E) @ vec(VECTOR_FIELD), SHAPE)
+    matrix_back = unvec(fd_operators.curl_back(DX_H) @ vec(VECTOR_FIELD), SHAPE)
+
+    assert_fields_close(curl_forward, matrix_forward, atol=1e-12, rtol=1e-12)
+    assert_fields_close(curl_back, matrix_back, atol=1e-12, rtol=1e-12)

meanas/test/test_fdmath_operators.py

@@ -0,0 +1,90 @@
+import numpy
+import pytest
+
+from ..fdmath import operators, unvec, vec
+from ._test_builders import real_ramp
+from .utils import assert_close
+
+
+SHAPE = (2, 3, 2)
+SCALAR_FIELD = real_ramp(SHAPE)
+VECTOR_LEFT = real_ramp((3, *SHAPE), offset=0.5)
+VECTOR_RIGHT = real_ramp((3, *SHAPE), scale=1 / 3, offset=2.0)
+
+
+def _apply_scalar_matrix(op: operators.sparse.spmatrix) -> numpy.ndarray:
+    return (op @ SCALAR_FIELD.ravel(order='C')).reshape(SHAPE, order='C')
+
+
+def _mirrored_indices(size: int, shift_distance: int) -> numpy.ndarray:
+    indices = numpy.arange(size) + shift_distance
+    indices = numpy.where(indices >= size, 2 * size - indices - 1, indices)
+    indices = numpy.where(indices < 0, -1 - indices, indices)
+    return indices
+
+
+@pytest.mark.parametrize(('axis', 'shift_distance'), [(0, 1), (1, -1), (2, 1)])
+def test_shift_circ_matches_numpy_roll(axis: int, shift_distance: int) -> None:
+    matrix_result = _apply_scalar_matrix(operators.shift_circ(axis, SHAPE, shift_distance))
+    expected = numpy.roll(SCALAR_FIELD, -shift_distance, axis=axis)
+    assert_close(matrix_result, expected)
+
+
+@pytest.mark.parametrize(('axis', 'shift_distance'), [(0, 1), (1, -1), (2, 1)])
+def test_shift_with_mirror_matches_explicit_mirrored_indices(axis: int, shift_distance: int) -> None:
+    matrix_result = _apply_scalar_matrix(operators.shift_with_mirror(axis, SHAPE, shift_distance))
+    indices = [numpy.arange(length) for length in SHAPE]
+    indices[axis] = _mirrored_indices(SHAPE[axis], shift_distance)
+    expected = SCALAR_FIELD[numpy.ix_(*indices)]
+    assert_close(matrix_result, expected)
+
+
+@pytest.mark.parametrize(
+    ('args', 'message'),
+    [
+        ((0, (2,), 1), 'Invalid shape'),
+        ((3, SHAPE, 1), 'Invalid direction'),
+    ],
+    )
+def test_shift_circ_rejects_invalid_arguments(args: tuple[int, tuple[int, ...], int], message: str) -> None:
+    with pytest.raises(Exception, match=message):
+        operators.shift_circ(*args)
+
+
+@pytest.mark.parametrize(
+    ('args', 'message'),
+    [
+        ((0, (2,), 1), 'Invalid shape'),
+        ((3, SHAPE, 1), 'Invalid direction'),
+        ((0, SHAPE, SHAPE[0]), 'too large'),
+    ],
+    )
+def test_shift_with_mirror_rejects_invalid_arguments(args: tuple[int, tuple[int, ...], int], message: str) -> None:
+    with pytest.raises(Exception, match=message):
+        operators.shift_with_mirror(*args)
+
+
+def test_vec_cross_matches_pointwise_cross_product() -> None:
+    matrix_result = unvec(operators.vec_cross(vec(VECTOR_LEFT)) @ vec(VECTOR_RIGHT), SHAPE)
+    expected = numpy.empty_like(VECTOR_LEFT)
+    expected[0] = VECTOR_LEFT[1] * VECTOR_RIGHT[2] - VECTOR_LEFT[2] * VECTOR_RIGHT[1]
+    expected[1] = VECTOR_LEFT[2] * VECTOR_RIGHT[0] - VECTOR_LEFT[0] * VECTOR_RIGHT[2]
+    expected[2] = VECTOR_LEFT[0] * VECTOR_RIGHT[1] - VECTOR_LEFT[1] * VECTOR_RIGHT[0]
+    assert_close(matrix_result, expected)
+
+
+def test_avg_forward_matches_half_sum_with_forward_neighbor() -> None:
+    matrix_result = _apply_scalar_matrix(operators.avg_forward(1, SHAPE))
+    expected = 0.5 * (SCALAR_FIELD + numpy.roll(SCALAR_FIELD, -1, axis=1))
+    assert_close(matrix_result, expected)
+
+
+def test_avg_back_matches_half_sum_with_backward_neighbor() -> None:
+    matrix_result = _apply_scalar_matrix(operators.avg_back(1, SHAPE))
+    expected = 0.5 * (SCALAR_FIELD + numpy.roll(SCALAR_FIELD, 1, axis=1))
+    assert_close(matrix_result, expected)
+
+
+def test_avg_forward_rejects_invalid_shape() -> None:
+    with pytest.raises(Exception, match='Invalid shape'):
+        operators.avg_forward(0, (2,))

meanas/test/test_fdmath_vectorization.py

@@ -0,0 +1,46 @@
+import numpy
+
+from ..fdmath import unvec, vec
+from ._test_builders import complex_ramp, real_ramp
+from .utils import assert_close
+
+
+SHAPE = (2, 3, 2)
+FIELD = real_ramp((3, *SHAPE))
+COMPLEX_FIELD = complex_ramp((3, *SHAPE), imag_scale=0.5)
+
+
+def test_vec_and_unvec_return_none_for_none_input() -> None:
+    assert vec(None) is None
+    assert unvec(None, SHAPE) is None
+
+
+def test_real_field_round_trip_preserves_shape_and_values() -> None:
+    vector = vec(FIELD)
+    assert vector is not None
+    restored = unvec(vector, SHAPE)
+    assert restored is not None
+    assert restored.shape == (3, *SHAPE)
+    assert_close(restored, FIELD)
+
+
+def test_complex_field_round_trip_preserves_shape_and_values() -> None:
+    vector = vec(COMPLEX_FIELD)
+    assert vector is not None
+    restored = unvec(vector, SHAPE)
+    assert restored is not None
+    assert restored.shape == (3, *SHAPE)
+    assert_close(restored, COMPLEX_FIELD)
+
+
+def test_unvec_with_two_components_round_trips_vector() -> None:
+    vector = numpy.arange(2 * numpy.prod(SHAPE), dtype=float)
+    field = unvec(vector, SHAPE, nvdim=2)
+    assert field is not None
+    assert field.shape == (2, *SHAPE)
+    assert_close(vec(field), vector)
+
+
+def test_vec_accepts_arraylike_input() -> None:
+    arraylike = [[[1, 2], [3, 4]], [[5, 6], [7, 8]]]
+    assert_close(vec(arraylike), numpy.ravel(arraylike, order='C'))

meanas/test/test_fdtd.py

@@ -7,7 +7,7 @@ from numpy.typing import NDArray
 #from numpy.testing import assert_allclose, assert_array_equal
 
 from .. import fdtd
-from .utils import assert_close, assert_fields_close, PRNG
+from .utils import assert_close, assert_fields_close, make_prng
 from .conftest import FixtureRequest
 
 
@@ -179,13 +179,14 @@ def j_distribution(
         shape: tuple[int, ...],
         j_mag: float,
         ) -> NDArray[numpy.float64]:
+    prng = make_prng()
     j = numpy.zeros(shape)
     if request.param == 'center':
         j[:, shape[1] // 2, shape[2] // 2, shape[3] // 2] = j_mag
     elif request.param == '000':
         j[:, 0, 0, 0] = j_mag
     elif request.param == 'random':
-        j[:] = PRNG.uniform(low=-j_mag, high=j_mag, size=shape)
+        j[:] = prng.uniform(low=-j_mag, high=j_mag, size=shape)
     return j
 
 

meanas/test/test_fdtd_base.py

@@ -0,0 +1,45 @@
+import numpy
+
+from ..fdmath import functional as fd_functional
+from ..fdtd import base
+from ._test_builders import real_ramp
+from .utils import assert_close
+
+
+DT = 0.25
+SHAPE = (3, 2, 2, 2)
+E_FIELD = real_ramp(SHAPE, scale=1 / 5)
+H_FIELD = real_ramp(SHAPE, scale=1 / 7, offset=1 / 7)
+EPSILON = 1.5 + E_FIELD / 10.0
+MU_FIELD = 2.0 + H_FIELD / 8.0
+MU_SCALAR = 3.0
+
+
+def test_maxwell_e_without_dxes_matches_unit_spacing_update() -> None:
+    updater = base.maxwell_e(dt=DT)
+    expected = E_FIELD + DT * fd_functional.curl_back()(H_FIELD) / EPSILON
+
+    updated = updater(E_FIELD.copy(), H_FIELD.copy(), EPSILON)
+
+    assert_close(updated, expected)
+
+
+def test_maxwell_h_without_dxes_and_without_mu_matches_unit_spacing_update() -> None:
+    updater = base.maxwell_h(dt=DT)
+    expected = H_FIELD - DT * fd_functional.curl_forward()(E_FIELD)
+
+    updated = updater(E_FIELD.copy(), H_FIELD.copy())
+
+    assert_close(updated, expected)
+
+
+def test_maxwell_h_without_dxes_accepts_scalar_and_field_mu() -> None:
+    updater = base.maxwell_h(dt=DT)
+
+    updated_scalar = updater(E_FIELD.copy(), H_FIELD.copy(), MU_SCALAR)
+    expected_scalar = H_FIELD - DT * fd_functional.curl_forward()(E_FIELD) / MU_SCALAR
+    assert_close(updated_scalar, expected_scalar)
+
+    updated_field = updater(E_FIELD.copy(), H_FIELD.copy(), MU_FIELD)
+    expected_field = H_FIELD - DT * fd_functional.curl_forward()(E_FIELD) / MU_FIELD
+    assert_close(updated_field, expected_field)

meanas/test/test_fdtd_boundaries.py

@@ -0,0 +1,62 @@
+import numpy
+import pytest
+from numpy.testing import assert_allclose
+
+from ..fdtd.boundaries import conducting_boundary
+
+
+def _axis_index(axis: int, index: int) -> tuple[slice | int, ...]:
+    coords: list[slice | int] = [slice(None), slice(None), slice(None)]
+    coords[axis] = index
+    return tuple(coords)
+
+
+@pytest.mark.parametrize('direction', [0, 1, 2])
+@pytest.mark.parametrize('polarity', [-1, 1])
+def test_conducting_boundary_updates_expected_faces(direction: int, polarity: int) -> None:
+    e = numpy.arange(3 * 4 * 4 * 4, dtype=float).reshape(3, 4, 4, 4)
+    h = e.copy()
+    e0 = e.copy()
+    h0 = h.copy()
+
+    update_e, update_h = conducting_boundary(direction, polarity)
+    update_e(e)
+    update_h(h)
+
+    dirs = [0, 1, 2]
+    dirs.remove(direction)
+    u, v = dirs
+
+    if polarity < 0:
+        boundary = _axis_index(direction, 0)
+        shifted1 = _axis_index(direction, 1)
+
+        assert_allclose(e[direction][boundary], 0)
+        assert_allclose(e[u][boundary], e0[u][shifted1])
+        assert_allclose(e[v][boundary], e0[v][shifted1])
+        assert_allclose(h[direction][boundary], h0[direction][shifted1])
+        assert_allclose(h[u][boundary], 0)
+        assert_allclose(h[v][boundary], 0)
+    else:
+        boundary = _axis_index(direction, -1)
+        shifted1 = _axis_index(direction, -2)
+        shifted2 = _axis_index(direction, -3)
+
+        assert_allclose(e[direction][boundary], -e0[direction][shifted2])
+        assert_allclose(e[direction][shifted1], 0)
+        assert_allclose(e[u][boundary], e0[u][shifted1])
+        assert_allclose(e[v][boundary], e0[v][shifted1])
+        assert_allclose(h[direction][boundary], h0[direction][shifted1])
+        assert_allclose(h[u][boundary], -h0[u][shifted2])
+        assert_allclose(h[u][shifted1], 0)
+        assert_allclose(h[v][boundary], -h0[v][shifted2])
+        assert_allclose(h[v][shifted1], 0)
+
+
+@pytest.mark.parametrize(
+    ('direction', 'polarity'),
+    [(-1, 1), (3, 1), (0, 0)],
+)
+def test_conducting_boundary_rejects_invalid_arguments(direction: int, polarity: int) -> None:
+    with pytest.raises(Exception):
+        conducting_boundary(direction, polarity)

meanas/test/test_fdtd_energy.py

@@ -0,0 +1,99 @@
+import numpy
+
+from .. import fdtd
+from ..fdtd import energy as fdtd_energy
+from ._test_builders import real_ramp, unit_dxes
+from .utils import assert_close
+
+
+SHAPE = (2, 2, 2)
+DT = 0.25
+UNIT_DXES = unit_dxes(SHAPE)
+DXES = (
+    (
+        numpy.array([1.0, 1.5]),
+        numpy.array([0.75, 1.25]),
+        numpy.array([1.1, 0.9]),
+    ),
+    (
+        numpy.array([0.8, 1.2]),
+        numpy.array([1.4, 0.6]),
+        numpy.array([0.7, 1.3]),
+    ),
+)
+E0 = real_ramp((3, *SHAPE))
+E1 = E0 + 0.5
+E2 = E0 + 1.0
+E3 = E0 + 1.5
+H0 = real_ramp((3, *SHAPE), scale=1 / 3, offset=2 / 3)
+H1 = H0 + 0.25
+H2 = H0 + 0.5
+H3 = H0 + 0.75
+J0 = (E0 + 2.0) / 5.0
+EPSILON = 1.0 + E0 / 20.0
+MU = 1.5 + H0 / 10.0
+
+
+def test_poynting_default_spacing_matches_explicit_unit_spacing() -> None:
+    default_spacing = fdtd.poynting(E1, H1)
+    explicit_spacing = fdtd.poynting(E1, H1, dxes=UNIT_DXES)
+    assert_close(default_spacing, explicit_spacing)
+
+
+def test_poynting_divergence_matches_precomputed_poynting_vector() -> None:
+    s = fdtd.poynting(E2, H2, dxes=DXES)
+    from_fields = fdtd.poynting_divergence(e=E2, h=H2, dxes=DXES)
+    from_vector = fdtd.poynting_divergence(s=s)
+    assert_close(from_fields, from_vector)
+
+
+def test_delta_energy_h2e_matches_direct_dxmul_formula() -> None:
+    expected = fdtd_energy.dxmul(
+        E2 * (E2 - E0) / DT,
+        H1 * (H3 - H1) / DT,
+        EPSILON,
+        MU,
+        DXES,
+    )
+    actual = fdtd.delta_energy_h2e(
+        dt=DT,
+        e0=E0,
+        h1=H1,
+        e2=E2,
+        h3=H3,
+        epsilon=EPSILON,
+        mu=MU,
+        dxes=DXES,
+    )
+    assert_close(actual, expected)
+
+
+def test_delta_energy_e2h_matches_direct_dxmul_formula() -> None:
+    expected = fdtd_energy.dxmul(
+        E1 * (E3 - E1) / DT,
+        H2 * (H2 - H0) / DT,
+        EPSILON,
+        MU,
+        DXES,
+    )
+    actual = fdtd_energy.delta_energy_e2h(
+        dt=DT,
+        h0=H0,
+        e1=E1,
+        h2=H2,
+        e3=E3,
+        epsilon=EPSILON,
+        mu=MU,
+        dxes=DXES,
+    )
+    assert_close(actual, expected)
+
+
+def test_delta_energy_j_defaults_to_unit_cell_volume() -> None:
+    expected = (J0 * E1).sum(axis=0)
+    assert_close(fdtd.delta_energy_j(j0=J0, e1=E1), expected)
+
+
+def test_dxmul_defaults_to_unit_materials_and_spacing() -> None:
+    expected = E1.sum(axis=0) + H1.sum(axis=0)
+    assert_close(fdtd_energy.dxmul(E1, H1), expected)

meanas/test/test_fdtd_misc.py

@@ -0,0 +1,42 @@
+import numpy
+import pytest
+
+from ..fdtd.misc import gaussian_beam, gaussian_packet, ricker_pulse
+
+
+@pytest.mark.parametrize('one_sided', [False, True])
+def test_gaussian_packet_accepts_array_input(one_sided: bool) -> None:
+    dt = 0.01
+    source, delay = gaussian_packet(1.55, 0.1, dt, one_sided=one_sided)
+    steps = numpy.array([0, int(numpy.ceil(delay / dt)) + 5])
+    envelope, cc, ss = source(steps)
+
+    assert envelope.shape == (2,)
+    assert numpy.isfinite(envelope).all()
+    assert numpy.isfinite(cc).all()
+    assert numpy.isfinite(ss).all()
+    if one_sided:
+        assert envelope[-1] == pytest.approx(1.0)
+
+
+def test_ricker_pulse_returns_finite_values() -> None:
+    source, delay = ricker_pulse(1.55, 0.01)
+    envelope, cc, ss = source(numpy.array([0, 1, 2]))
+
+    assert numpy.isfinite(delay)
+    assert numpy.isfinite(envelope).all()
+    assert numpy.isfinite(cc).all()
+    assert numpy.isfinite(ss).all()
+
+
+def test_gaussian_beam_centered_grid_is_finite_and_normalized() -> None:
+    beam = gaussian_beam(
+        xyz=[numpy.linspace(-1, 1, 3), numpy.linspace(-1, 1, 3), numpy.linspace(-1, 1, 3)],
+        center=[0, 0, 0],
+        waist_radius=1.0,
+        wl=1.55,
+    )
+
+    row = beam[:, :, beam.shape[2] // 2]
+    assert numpy.isfinite(beam).all()
+    assert numpy.linalg.norm(row) == pytest.approx(1.0)

meanas/test/test_fdtd_phasor.py

@@ -0,0 +1,208 @@
+import dataclasses
+from functools import lru_cache
+
+import numpy
+import pytest
+import scipy.sparse.linalg
+
+from .. import fdfd, fdtd
+from ..fdmath import unvec, vec
+from ._test_builders import unit_dxes
+from .utils import assert_close, assert_fields_close
+
+
+@dataclasses.dataclass(frozen=True)
+class ContinuousWaveCase:
+    omega: float
+    dt: float
+    dxes: tuple[tuple[numpy.ndarray, ...], tuple[numpy.ndarray, ...]]
+    epsilon: numpy.ndarray
+    e_ph: numpy.ndarray
+    h_ph: numpy.ndarray
+    j_ph: numpy.ndarray
+
+
+def test_phasor_accumulator_matches_direct_sum_for_multi_frequency_weights() -> None:
+    omegas = numpy.array([0.25, 0.5])
+    dt = 0.2
+    sample_0 = numpy.array([[1.0, 2.0], [3.0, 4.0]])
+    sample_1 = numpy.array([[0.5, 1.5], [2.5, 3.5]])
+    weight_0 = numpy.array([1.0, 2.0])
+    weight_1 = 0.75
+    accumulator = numpy.zeros((omegas.size, *sample_0.shape), dtype=complex)
+
+    fdtd.accumulate_phasor(accumulator, omegas, dt, sample_0, 0, weight=weight_0)
+    fdtd.accumulate_phasor(accumulator, omegas, dt, sample_1, 3, offset_steps=0.5, weight=weight_1)
+
+    expected = numpy.zeros((2, *sample_0.shape), dtype=complex)
+    for idx, omega in enumerate(omegas):
+        expected[idx] += dt * weight_0[idx] * numpy.exp(-1j * omega * 0.0) * sample_0
+        expected[idx] += dt * weight_1 * numpy.exp(-1j * omega * ((3 + 0.5) * dt)) * sample_1
+
+    assert_close(accumulator, expected)
+
+
+def test_phasor_accumulator_convenience_methods_apply_yee_offsets() -> None:
+    omega = 1.25
+    dt = 0.1
+    sample = numpy.arange(6, dtype=float).reshape(2, 3)
+    e_acc = numpy.zeros((1, *sample.shape), dtype=complex)
+    h_acc = numpy.zeros((1, *sample.shape), dtype=complex)
+    j_acc = numpy.zeros((1, *sample.shape), dtype=complex)
+
+    fdtd.accumulate_phasor_e(e_acc, omega, dt, sample, 4)
+    fdtd.accumulate_phasor_h(h_acc, omega, dt, sample, 4)
+    fdtd.accumulate_phasor_j(j_acc, omega, dt, sample, 4)
+
+    expected_e = dt * numpy.exp(-1j * omega * (4 * dt)) * sample
+    expected_h = dt * numpy.exp(-1j * omega * ((4.5) * dt)) * sample
+
+    assert_close(e_acc[0], expected_e)
+    assert_close(h_acc[0], expected_h)
+    assert_close(j_acc[0], expected_h)
+
+
+def test_phasor_accumulator_matches_delayed_weighted_example_pattern() -> None:
+    omega = 0.75
+    dt = 0.2
+    delay = 0.6
+    phasor_norm = 0.5
+    steps = numpy.arange(5)
+    samples = numpy.arange(20, dtype=float).reshape(5, 2, 2) + 1.0
+    accumulator = numpy.zeros((1, 2, 2), dtype=complex)
+
+    for step, sample in zip(steps, samples, strict=True):
+        fdtd.accumulate_phasor(
+            accumulator,
+            omega,
+            dt,
+            sample,
+            int(step),
+            offset_steps=0.5 - delay / dt,
+            weight=phasor_norm / dt,
+        )
+
+    expected = numpy.zeros((2, 2), dtype=complex)
+    for step, sample in zip(steps, samples, strict=True):
+        time = (step + 0.5 - delay / dt) * dt
+        expected += dt * (phasor_norm / dt) * numpy.exp(-1j * omega * time) * sample
+
+    assert_close(accumulator[0], expected)
+
+
+def test_phasor_accumulator_validation_reset_and_squeeze() -> None:
+    with pytest.raises(ValueError, match='dt must be positive'):
+        fdtd.accumulate_phasor(numpy.zeros((1, 2, 2), dtype=complex), [1.0], 0.0, numpy.ones((2, 2)), 0)
+
+    with pytest.raises(ValueError, match='omegas must be a scalar or non-empty 1D sequence'):
+        fdtd.accumulate_phasor(numpy.zeros((1, 2, 2), dtype=complex), numpy.ones((2, 2)), 0.2, numpy.ones((2, 2)), 0)
+
+    accumulator = numpy.zeros((2, 2, 2), dtype=complex)
+
+    with pytest.raises(ValueError, match='accumulator must have shape'):
+        fdtd.accumulate_phasor(accumulator, [1.0], 0.2, numpy.ones((2, 2)), 0)
+
+    with pytest.raises(ValueError, match='weight must be scalar'):
+        fdtd.accumulate_phasor(accumulator, [1.0, 2.0], 0.2, numpy.ones((2, 2)), 0, weight=numpy.ones((2, 2)))
+
+    fdtd.accumulate_phasor(accumulator, [1.0, 2.0], 0.2, numpy.ones((2, 2)), 0)
+    accumulator.fill(0)
+    assert_close(accumulator, 0.0)
+
+
+@lru_cache(maxsize=1)
+def _continuous_wave_case() -> ContinuousWaveCase:
+    spatial_shape = (5, 1, 5)
+    full_shape = (3, *spatial_shape)
+    dt = 0.25
+    period_steps = 64
+    warmup_periods = 8
+    accumulation_periods = 8
+    omega = 2 * numpy.pi / (period_steps * dt)
+    total_steps = period_steps * (warmup_periods + accumulation_periods)
+    warmup_steps = period_steps * warmup_periods
+    accumulation_steps = period_steps * accumulation_periods
+    source_amplitude = 1.0
+    source_index = (1, spatial_shape[0] // 2, spatial_shape[1] // 2, spatial_shape[2] // 2)
+
+    dxes = unit_dxes(spatial_shape)
+    epsilon = numpy.ones(full_shape, dtype=float)
+    e_field = numpy.zeros(full_shape, dtype=float)
+    h_field = numpy.zeros(full_shape, dtype=float)
+    update_e = fdtd.maxwell_e(dt=dt, dxes=dxes)
+    update_h = fdtd.maxwell_h(dt=dt, dxes=dxes)
+
+    e_accumulator = numpy.zeros((1, *full_shape), dtype=complex)
+    h_accumulator = numpy.zeros((1, *full_shape), dtype=complex)
+    j_accumulator = numpy.zeros((1, *full_shape), dtype=complex)
+
+    for step in range(total_steps):
+        update_e(e_field, h_field, epsilon)
+
+        j_step = numpy.zeros_like(e_field)
+        current_density = source_amplitude * numpy.cos(omega * (step + 0.5) * dt)
+        j_step[source_index] = current_density
+        e_field -= dt * j_step / epsilon
+
+        if step >= warmup_steps:
+            fdtd.accumulate_phasor_j(j_accumulator, omega, dt, j_step, step)
+            fdtd.accumulate_phasor_e(e_accumulator, omega, dt, e_field, step + 1)
+
+        update_h(e_field, h_field)
+
+        if step >= warmup_steps:
+            fdtd.accumulate_phasor_h(h_accumulator, omega, dt, h_field, step + 1)
+
+    return ContinuousWaveCase(
+        omega=omega,
+        dt=dt,
+        dxes=dxes,
+        epsilon=epsilon,
+        e_ph=e_accumulator[0],
+        h_ph=h_accumulator[0],
+        j_ph=j_accumulator[0],
+    )
+
+
+def test_continuous_wave_current_phasor_matches_analytic_discrete_sum() -> None:
+    case = _continuous_wave_case()
+
+    accumulation_indices = numpy.arange(64 * 8, 64 * 16)
+    times = (accumulation_indices + 0.5) * case.dt
+    expected = numpy.zeros_like(case.j_ph)
+    expected[1, 2, 0, 2] = case.dt * numpy.sum(
+        numpy.exp(-1j * case.omega * times) * numpy.cos(case.omega * times),
+    )
+
+    assert_fields_close(case.j_ph, expected, atol=1e-12, rtol=1e-12)
+
+
+def test_continuous_wave_electric_phasor_matches_fdfd_solution() -> None:
+    case = _continuous_wave_case()
+    operator = fdfd.operators.e_full(case.omega, case.dxes, vec(case.epsilon)).tocsr()
+    rhs = -1j * case.omega * vec(case.j_ph)
+    e_fdfd = unvec(scipy.sparse.linalg.spsolve(operator, rhs), case.epsilon.shape[1:])
+
+    rel_err = numpy.linalg.norm(vec(case.e_ph - e_fdfd)) / numpy.linalg.norm(vec(e_fdfd))
+    assert rel_err < 5e-2
+
+
+def test_continuous_wave_magnetic_phasor_matches_fdfd_conversion() -> None:
+    case = _continuous_wave_case()
+    operator = fdfd.operators.e_full(case.omega, case.dxes, vec(case.epsilon)).tocsr()
+    rhs = -1j * case.omega * vec(case.j_ph)
+    e_fdfd = unvec(scipy.sparse.linalg.spsolve(operator, rhs), case.epsilon.shape[1:])
+    h_fdfd = fdfd.functional.e2h(case.omega, case.dxes)(e_fdfd)
+
+    rel_err = numpy.linalg.norm(vec(case.h_ph - h_fdfd)) / numpy.linalg.norm(vec(h_fdfd))
+    assert rel_err < 5e-2
+
+
+def test_continuous_wave_extracted_electric_phasor_has_small_fdfd_residual() -> None:
+    case = _continuous_wave_case()
+    operator = fdfd.operators.e_full(case.omega, case.dxes, vec(case.epsilon)).tocsr()
+    rhs = -1j * case.omega * vec(case.j_ph)
+    residual = operator @ vec(case.e_ph) - rhs
+    rel_residual = numpy.linalg.norm(residual) / numpy.linalg.norm(rhs)
+
+    assert rel_residual < 5e-2

meanas/test/test_fdtd_pml.py

@@ -0,0 +1,44 @@
+import numpy
+import pytest
+
+from ..fdtd.pml import cpml_params, updates_with_cpml
+
+
+@pytest.mark.parametrize(
+    ('axis', 'polarity', 'thickness', 'epsilon_eff'),
+    [(3, 1, 4, 1.0), (0, 0, 4, 1.0), (0, 1, 2, 1.0), (0, 1, 4, 0.0)],
+)
+def test_cpml_params_reject_invalid_arguments(axis: int, polarity: int, thickness: int, epsilon_eff: float) -> None:
+    with pytest.raises(Exception):
+        cpml_params(axis=axis, polarity=polarity, dt=0.1, thickness=thickness, epsilon_eff=epsilon_eff)
+
+
+def test_cpml_params_shapes_and_region() -> None:
+    params = cpml_params(axis=1, polarity=1, dt=0.1, thickness=3)
+    p0e, p1e, p2e = params['param_e']
+    p0h, p1h, p2h = params['param_h']
+
+    assert p0e.shape == (1, 3, 1)
+    assert p1e.shape == (1, 3, 1)
+    assert p2e.shape == (1, 3, 1)
+    assert p0h.shape == (1, 3, 1)
+    assert p1h.shape == (1, 3, 1)
+    assert p2h.shape == (1, 3, 1)
+    assert params['region'][1] == slice(-3, None)
+
+
+def test_updates_with_cpml_keeps_zero_fields_zero() -> None:
+    shape = (3, 4, 4, 4)
+    epsilon = numpy.ones(shape, dtype=float)
+    e = numpy.zeros(shape, dtype=float)
+    h = numpy.zeros(shape, dtype=float)
+    dxes = [[numpy.ones(4), numpy.ones(4), numpy.ones(4)] for _ in range(2)]
+    params = [[None, None] for _ in range(3)]
+    params[0][0] = cpml_params(axis=0, polarity=-1, dt=0.1, thickness=3)
+
+    update_e, update_h = updates_with_cpml(params, dt=0.1, dxes=dxes, epsilon=epsilon)
+    update_e(e, h, epsilon)
+    update_h(e, h)
+
+    assert not e.any()
+    assert not h.any()

meanas/test/test_import_fallbacks.py

@@ -0,0 +1,51 @@
+import builtins
+import importlib
+import pathlib
+
+import meanas
+from ..fdfd import bloch
+from .utils import assert_close
+
+
+def _reload(module):
+    return importlib.reload(module)
+
+
+def _restore_reloaded(monkeypatch, module):
+    monkeypatch.undo()
+    return _reload(module)
+
+
+def test_meanas_import_survives_readme_open_failure(monkeypatch) -> None:  # type: ignore[no-untyped-def]
+    original_open = pathlib.Path.open
+
+    def failing_open(self: pathlib.Path, *args, **kwargs):  # type: ignore[no-untyped-def]
+        if self.name == 'README.md':
+            raise FileNotFoundError('forced README failure')
+        return original_open(self, *args, **kwargs)
+
+    monkeypatch.setattr(pathlib.Path, 'open', failing_open)
+    reloaded = _reload(meanas)
+
+    assert reloaded.__version__ == '0.10'
+    assert reloaded.__author__ == 'Jan Petykiewicz'
+    assert reloaded.__doc__ is not None
+
+    _restore_reloaded(monkeypatch, meanas)
+
+
+def test_bloch_reloads_with_numpy_fft_when_pyfftw_is_unavailable(monkeypatch) -> None:  # type: ignore[no-untyped-def]
+    original_import = builtins.__import__
+
+    def fake_import(name: str, globals=None, locals=None, fromlist=(), level: int = 0):  # type: ignore[no-untyped-def]
+        if name.startswith('pyfftw'):
+            raise ImportError('forced pyfftw failure')
+        return original_import(name, globals, locals, fromlist, level)
+
+    monkeypatch.setattr(builtins, '__import__', fake_import)
+    reloaded = _reload(bloch)
+
+    assert reloaded.fftn.__module__ == 'numpy.fft'
+    assert reloaded.ifftn.__module__ == 'numpy.fft'
+
+    _restore_reloaded(monkeypatch, bloch)

meanas/test/test_waveguide_2d_numerics.py

@@ -0,0 +1,284 @@
+import numpy
+from numpy.linalg import norm
+from numpy.testing import assert_allclose
+
+from ..fdmath import vec
+from ..fdfd import waveguide_2d
+
+
+OMEGA = 1 / 1500
+GRID_SHAPE = (5, 5)
+DXES_2D = [[numpy.ones(GRID_SHAPE[0]), numpy.ones(GRID_SHAPE[1])] for _ in range(2)]
+DXES_2D_NONUNIFORM = [[
+    numpy.array([1.0, 1.2, 0.9, 1.1, 1.3]),
+    numpy.array([0.8, 1.1, 1.0, 1.2, 0.9]),
+] for _ in range(2)]
+
+
+def build_asymmetric_epsilon() -> numpy.ndarray:
+    epsilon = numpy.ones((3, *GRID_SHAPE), dtype=float)
+    epsilon[:, 2, 1] = 2.0
+    return vec(epsilon)
+
+
+def build_mu_profile() -> numpy.ndarray:
+    return numpy.linspace(1.5, 2.2, 3 * GRID_SHAPE[0] * GRID_SHAPE[1])
+
+
+def test_waveguide_2d_solved_modes_are_ordered_and_low_residual() -> None:
+    epsilon = build_asymmetric_epsilon()
+    operator_e = waveguide_2d.operator_e(OMEGA, DXES_2D, epsilon)
+
+    e_xys, wavenumbers = waveguide_2d.solve_modes(
+        [0, 1],
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        )
+
+    assert numpy.all(numpy.diff(numpy.real(wavenumbers)) <= 0)
+
+    for e_xy, wavenumber in zip(e_xys, wavenumbers, strict=True):
+        residual = norm(operator_e @ e_xy - (wavenumber ** 2) * e_xy) / norm(e_xy)
+        assert residual < 1e-6
+
+
+def test_waveguide_2d_normalized_fields_are_consistent() -> None:
+    epsilon = build_asymmetric_epsilon()
+    operator_h = waveguide_2d.operator_h(OMEGA, DXES_2D, epsilon)
+
+    e_xy, wavenumber = waveguide_2d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        )
+    e_field, h_field = waveguide_2d.normalized_fields_e(
+        e_xy,
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        )
+    h_xy = numpy.concatenate(numpy.split(h_field, 3)[:2])
+
+    overlap = waveguide_2d.inner_product(e_field, h_field, DXES_2D, conj_h=True)
+    h_residual = norm(operator_h @ h_xy - (wavenumber ** 2) * h_xy) / norm(h_xy)
+
+    assert abs(overlap.real - 1.0) < 1e-10
+    assert abs(overlap.imag) < 1e-10
+    assert waveguide_2d.e_err(e_field, wavenumber, OMEGA, DXES_2D, epsilon) < 1e-6
+    assert waveguide_2d.h_err(h_field, wavenumber, OMEGA, DXES_2D, epsilon) < 1e-6
+    assert h_residual < 1e-6
+
+
+def test_waveguide_2d_sensitivity_matches_finite_difference() -> None:
+    epsilon = build_asymmetric_epsilon()
+    e_xy, wavenumber = waveguide_2d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        )
+    e_field, h_field = waveguide_2d.normalized_fields_e(
+        e_xy,
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        )
+    sensitivity = waveguide_2d.sensitivity(
+        e_field,
+        h_field,
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        )
+
+    target_index = int(numpy.argmax(numpy.abs(sensitivity)))
+    delta = 1e-4
+    epsilon_perturbed = epsilon.copy()
+    epsilon_perturbed[target_index] += delta
+
+    _, perturbed_wavenumber = waveguide_2d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon_perturbed,
+        )
+    finite_difference = (perturbed_wavenumber - wavenumber) / delta
+
+    assert numpy.isfinite(sensitivity[target_index])
+    assert numpy.isfinite(finite_difference)
+    assert abs(sensitivity[target_index].imag) < 1e-10
+    assert abs(finite_difference.imag) < 1e-10
+
+    ratio = abs(sensitivity[target_index] / finite_difference)
+    assert sensitivity[target_index].real > 0
+    assert finite_difference.real > 0
+    assert 0.4 < ratio < 1.8
+
+
+def test_waveguide_2d_normalized_fields_h_are_finite_and_unit_normalized_with_mu() -> None:
+    epsilon = build_asymmetric_epsilon()
+    mu = build_mu_profile()
+
+    e_xy, wavenumber = waveguide_2d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        )
+    _e_ref, h_ref = waveguide_2d.normalized_fields_e(
+        e_xy,
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        mu=mu,
+        )
+    h_xy = numpy.concatenate(numpy.split(h_ref, 3)[:2])
+
+    e_field, h_field = waveguide_2d.normalized_fields_h(
+        h_xy,
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        mu=mu,
+        )
+    overlap = waveguide_2d.inner_product(e_field, h_field, DXES_2D, conj_h=True)
+
+    assert e_field.shape == (3 * GRID_SHAPE[0] * GRID_SHAPE[1],)
+    assert h_field.shape == (3 * GRID_SHAPE[0] * GRID_SHAPE[1],)
+    assert numpy.isfinite(e_field).all()
+    assert numpy.isfinite(h_field).all()
+    assert abs(overlap.real - 1.0) < 1e-10
+    assert abs(overlap.imag) < 1e-10
+
+
+def test_waveguide_2d_helper_operators_with_mu_return_finite_outputs() -> None:
+    epsilon = build_asymmetric_epsilon()
+    mu = build_mu_profile()
+
+    e_xy, wavenumber = waveguide_2d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        )
+    _e_ref, h_ref = waveguide_2d.normalized_fields_e(
+        e_xy,
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        mu=mu,
+        )
+    h_xy = numpy.concatenate(numpy.split(h_ref, 3)[:2])
+    n_pts = GRID_SHAPE[0] * GRID_SHAPE[1]
+
+    operators = [
+        ('exy2h', waveguide_2d.exy2h(wavenumber, OMEGA, DXES_2D, epsilon, mu), e_xy),
+        ('hxy2e', waveguide_2d.hxy2e(wavenumber, OMEGA, DXES_2D, epsilon, mu), h_xy),
+        ('hxy2h', waveguide_2d.hxy2h(wavenumber, DXES_2D, mu), h_xy),
+    ]
+
+    for _name, operator, vector in operators:
+        result = operator @ vector
+        assert operator.shape == (3 * n_pts, 2 * n_pts)
+        assert numpy.isfinite(operator.data).all()
+        assert result.shape == (3 * n_pts,)
+        assert numpy.isfinite(result).all()
+
+
+def test_waveguide_2d_error_helpers_with_mu_return_finite_values() -> None:
+    epsilon = build_asymmetric_epsilon()
+    mu = build_mu_profile()
+
+    e_xy, wavenumber = waveguide_2d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        )
+    e_field, h_field = waveguide_2d.normalized_fields_e(
+        e_xy,
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        mu=mu,
+        )
+
+    h_error = waveguide_2d.h_err(h_field, wavenumber, OMEGA, DXES_2D, epsilon, mu)
+    e_error = waveguide_2d.e_err(e_field, wavenumber, OMEGA, DXES_2D, epsilon, mu)
+
+    assert numpy.isfinite(h_error)
+    assert numpy.isfinite(e_error)
+    assert 0 < h_error < 0.1
+    assert 0 < e_error < 2.0
+
+
+def test_waveguide_2d_inner_product_phase_and_conjugation_branches() -> None:
+    epsilon = build_asymmetric_epsilon()
+    mu = build_mu_profile()
+
+    e_xy, wavenumber = waveguide_2d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        )
+    e_field, h_field = waveguide_2d.normalized_fields_e(
+        e_xy,
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=DXES_2D,
+        epsilon=epsilon,
+        mu=mu,
+        )
+
+    overlap_no_conj = waveguide_2d.inner_product(e_field, h_field, DXES_2D, conj_h=False)
+    overlap_conj = waveguide_2d.inner_product(e_field, h_field, DXES_2D, conj_h=True)
+    overlap_phase = waveguide_2d.inner_product(e_field, h_field, DXES_2D, conj_h=True, prop_phase=0.3)
+
+    assert numpy.isfinite(overlap_no_conj)
+    assert numpy.isfinite(overlap_phase)
+    assert abs(overlap_no_conj.real - 1.0) < 1e-10
+    assert abs(overlap_no_conj.imag) < 1e-10
+    assert_allclose(overlap_phase / overlap_conj, numpy.exp(-0.15j), atol=1e-12, rtol=1e-12)
+
+
+def test_waveguide_2d_inner_product_trapezoid_branch_on_nonuniform_grid() -> None:
+    epsilon = build_asymmetric_epsilon()
+
+    e_xy, wavenumber = waveguide_2d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=DXES_2D_NONUNIFORM,
+        epsilon=epsilon,
+        )
+    e_field, h_field = waveguide_2d.normalized_fields_e(
+        e_xy,
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=DXES_2D_NONUNIFORM,
+        epsilon=epsilon,
+        )
+
+    overlap_rect = waveguide_2d.inner_product(e_field, h_field, DXES_2D_NONUNIFORM, conj_h=True)
+    overlap_trap = waveguide_2d.inner_product(
+        e_field,
+        h_field,
+        DXES_2D_NONUNIFORM,
+        conj_h=True,
+        trapezoid=True,
+        )
+
+    assert numpy.isfinite(overlap_rect)
+    assert numpy.isfinite(overlap_trap)
+    assert abs(overlap_rect.imag) < 1e-10
+    assert abs(overlap_trap.imag) < 1e-10
+    assert abs(overlap_rect - overlap_trap) > 0.1

meanas/test/test_waveguide_fdtd_fdfd.py

@@ -0,0 +1,436 @@
+import dataclasses
+from functools import lru_cache
+
+import numpy
+
+from .. import fdfd, fdtd
+from ..fdmath import vec, unvec
+from ..fdfd import functional, scpml, waveguide_3d
+
+
+DT = 0.25
+PERIOD_STEPS = 64
+OMEGA = 2 * numpy.pi / (PERIOD_STEPS * DT)
+CPML_THICKNESS = 3
+WARMUP_PERIODS = 9
+ACCUMULATION_PERIODS = 9
+SHAPE = (3, 25, 13, 13)
+SOURCE_SLICES = (slice(4, 5), slice(None), slice(None))
+MONITOR_SLICES = (slice(18, 19), slice(None), slice(None))
+CHOSEN_VARIANT = 'base'
+SCATTERING_SHAPE = (3, 35, 15, 15)
+SCATTERING_SOURCE_SLICES = (slice(4, 5), slice(None), slice(None))
+SCATTERING_REFLECT_SLICES = (slice(10, 11), slice(None), slice(None))
+SCATTERING_TRANSMIT_SLICES = (slice(29, 30), slice(None), slice(None))
+SCATTERING_STEP_X = 18
+SCATTERING_WARMUP_PERIODS = 10
+SCATTERING_ACCUMULATION_PERIODS = 10
+
+
+@dataclasses.dataclass(frozen=True)
+class WaveguideCalibrationResult:
+    variant: str
+    e_ph: numpy.ndarray
+    h_ph: numpy.ndarray
+    j_ph: numpy.ndarray
+    e_fdfd: numpy.ndarray
+    h_fdfd: numpy.ndarray
+    overlap_td: complex
+    overlap_fd: complex
+    flux_td: float
+    flux_fd: float
+
+    @property
+    def overlap_rel_err(self) -> float:
+        return float(abs(self.overlap_td - self.overlap_fd) / abs(self.overlap_fd))
+
+    @property
+    def overlap_mag_rel_err(self) -> float:
+        return float(abs(abs(self.overlap_td) - abs(self.overlap_fd)) / abs(self.overlap_fd))
+
+    @property
+    def overlap_phase_deg(self) -> float:
+        return float(abs(numpy.degrees(numpy.angle(self.overlap_td / self.overlap_fd))))
+
+    @property
+    def flux_rel_err(self) -> float:
+        return float(abs(self.flux_td - self.flux_fd) / abs(self.flux_fd))
+
+    @property
+    def combined_error(self) -> float:
+        return self.overlap_mag_rel_err + self.flux_rel_err
+
+
+@dataclasses.dataclass(frozen=True)
+class WaveguideScatteringResult:
+    e_ph: numpy.ndarray
+    h_ph: numpy.ndarray
+    j_ph: numpy.ndarray
+    e_fdfd: numpy.ndarray
+    h_fdfd: numpy.ndarray
+    reflected_td: complex
+    reflected_fd: complex
+    transmitted_td: complex
+    transmitted_fd: complex
+    reflected_flux_td: float
+    reflected_flux_fd: float
+    transmitted_flux_td: float
+    transmitted_flux_fd: float
+
+    @property
+    def reflected_overlap_mag_rel_err(self) -> float:
+        return float(abs(abs(self.reflected_td) - abs(self.reflected_fd)) / abs(self.reflected_fd))
+
+    @property
+    def transmitted_overlap_mag_rel_err(self) -> float:
+        return float(abs(abs(self.transmitted_td) - abs(self.transmitted_fd)) / abs(self.transmitted_fd))
+
+    @property
+    def reflected_flux_rel_err(self) -> float:
+        return float(abs(self.reflected_flux_td - self.reflected_flux_fd) / abs(self.reflected_flux_fd))
+
+    @property
+    def transmitted_flux_rel_err(self) -> float:
+        return float(abs(self.transmitted_flux_td - self.transmitted_flux_fd) / abs(self.transmitted_flux_fd))
+
+
+def _build_uniform_dxes(shape: tuple[int, int, int, int]) -> list[list[numpy.ndarray]]:
+    return [[numpy.ones(shape[axis + 1]) for axis in range(3)] for _ in range(2)]
+
+
+def _build_base_dxes() -> list[list[numpy.ndarray]]:
+    return _build_uniform_dxes(SHAPE)
+
+
+def _build_stretched_dxes(base_dxes: list[list[numpy.ndarray]]) -> list[list[numpy.ndarray]]:
+    stretched_dxes = [[dx.copy() for dx in group] for group in base_dxes]
+    for axis in (0, 1, 2):
+        for polarity in (-1, 1):
+            stretched_dxes = scpml.stretch_with_scpml(
+                stretched_dxes,
+                axis=axis,
+                polarity=polarity,
+                omega=OMEGA,
+                epsilon_effective=1.0,
+                thickness=CPML_THICKNESS,
+            )
+    return stretched_dxes
+
+
+def _build_epsilon() -> numpy.ndarray:
+    epsilon = numpy.ones(SHAPE, dtype=float)
+    y0 = (SHAPE[2] - 3) // 2
+    z0 = (SHAPE[3] - 3) // 2
+    epsilon[:, :, y0:y0 + 3, z0:z0 + 3] = 12.0
+    return epsilon
+
+
+def _build_scattering_epsilon() -> numpy.ndarray:
+    epsilon = numpy.ones(SCATTERING_SHAPE, dtype=float)
+    y0 = SCATTERING_SHAPE[2] // 2
+    z0 = SCATTERING_SHAPE[3] // 2
+    epsilon[:, :SCATTERING_STEP_X, y0 - 1:y0 + 2, z0 - 1:z0 + 2] = 12.0
+    epsilon[:, SCATTERING_STEP_X:, y0 - 2:y0 + 3, z0 - 2:z0 + 3] = 12.0
+    return epsilon
+
+
+def _build_cpml_params() -> list[list[dict[str, numpy.ndarray | float]]]:
+    return [
+        [fdtd.cpml_params(axis=axis, polarity=polarity, dt=DT, thickness=CPML_THICKNESS, epsilon_eff=1.0)
+         for polarity in (-1, 1)]
+        for axis in range(3)
+    ]
+
+
+@lru_cache(maxsize=2)
+def _run_straight_waveguide_case(variant: str) -> WaveguideCalibrationResult:
+    assert variant in ('stretched', 'base')
+
+    epsilon = _build_epsilon()
+    base_dxes = _build_base_dxes()
+    stretched_dxes = _build_stretched_dxes(base_dxes)
+    mode_dxes = stretched_dxes if variant == 'stretched' else base_dxes
+
+    source_mode = waveguide_3d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=mode_dxes,
+        axis=0,
+        polarity=1,
+        slices=SOURCE_SLICES,
+        epsilon=epsilon,
+    )
+    j_mode = waveguide_3d.compute_source(
+        E=source_mode['E'],
+        wavenumber=source_mode['wavenumber'],
+        omega=OMEGA,
+        dxes=mode_dxes,
+        axis=0,
+        polarity=1,
+        slices=SOURCE_SLICES,
+        epsilon=epsilon,
+    )
+    monitor_mode = waveguide_3d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=mode_dxes,
+        axis=0,
+        polarity=1,
+        slices=MONITOR_SLICES,
+        epsilon=epsilon,
+    )
+    overlap_e = waveguide_3d.compute_overlap_e(
+        E=monitor_mode['E'],
+        wavenumber=monitor_mode['wavenumber'],
+        dxes=mode_dxes,
+        axis=0,
+        polarity=1,
+        slices=MONITOR_SLICES,
+        omega=OMEGA,
+    )
+
+    update_e, update_h = fdtd.updates_with_cpml(cpml_params=_build_cpml_params(), dt=DT, dxes=base_dxes, epsilon=epsilon)
+
+    e_field = numpy.zeros_like(epsilon)
+    h_field = numpy.zeros_like(epsilon)
+    e_accumulator = numpy.zeros((1, *SHAPE), dtype=complex)
+    h_accumulator = numpy.zeros((1, *SHAPE), dtype=complex)
+    j_accumulator = numpy.zeros((1, *SHAPE), dtype=complex)
+
+    warmup_steps = WARMUP_PERIODS * PERIOD_STEPS
+    accumulation_steps = ACCUMULATION_PERIODS * PERIOD_STEPS
+    for step in range(warmup_steps + accumulation_steps):
+        update_e(e_field, h_field, epsilon)
+
+        t_half = (step + 0.5) * DT
+        j_real = (j_mode.real * numpy.cos(OMEGA * t_half) - j_mode.imag * numpy.sin(OMEGA * t_half)).real
+        e_field -= DT * j_real / epsilon
+
+        if step >= warmup_steps:
+            fdtd.accumulate_phasor_j(j_accumulator, OMEGA, DT, j_real, step)
+            fdtd.accumulate_phasor_e(e_accumulator, OMEGA, DT, e_field, step + 1)
+
+        update_h(e_field, h_field)
+
+        if step >= warmup_steps:
+            fdtd.accumulate_phasor_h(h_accumulator, OMEGA, DT, h_field, step + 1)
+
+    e_ph = e_accumulator[0]
+    h_ph = h_accumulator[0]
+    j_ph = j_accumulator[0]
+
+    e_fdfd = unvec(
+        fdfd.solvers.generic(
+            J=vec(j_ph),
+            omega=OMEGA,
+            dxes=stretched_dxes,
+            epsilon=vec(epsilon),
+            matrix_solver_opts={'atol': 1e-10, 'rtol': 1e-7},
+        ),
+        SHAPE[1:],
+    )
+    h_fdfd = functional.e2h(OMEGA, stretched_dxes)(e_fdfd)
+
+    overlap_td = vec(e_ph) @ vec(overlap_e).conj()
+    overlap_fd = vec(e_fdfd) @ vec(overlap_e).conj()
+
+    poynting_td = functional.poynting_e_cross_h(stretched_dxes)(e_ph, h_ph.conj())
+    poynting_fd = functional.poynting_e_cross_h(stretched_dxes)(e_fdfd, h_fdfd.conj())
+    flux_td = float(0.5 * poynting_td[0, MONITOR_SLICES[0], :, :].real.sum())
+    flux_fd = float(0.5 * poynting_fd[0, MONITOR_SLICES[0], :, :].real.sum())
+
+    return WaveguideCalibrationResult(
+        variant=variant,
+        e_ph=e_ph,
+        h_ph=h_ph,
+        j_ph=j_ph,
+        e_fdfd=e_fdfd,
+        h_fdfd=h_fdfd,
+        overlap_td=overlap_td,
+        overlap_fd=overlap_fd,
+        flux_td=flux_td,
+        flux_fd=flux_fd,
+    )
+
+
+@lru_cache(maxsize=1)
+def _run_width_step_scattering_case() -> WaveguideScatteringResult:
+    epsilon = _build_scattering_epsilon()
+    base_dxes = _build_uniform_dxes(SCATTERING_SHAPE)
+    stretched_dxes = _build_stretched_dxes(base_dxes)
+
+    source_mode = waveguide_3d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=base_dxes,
+        axis=0,
+        polarity=1,
+        slices=SCATTERING_SOURCE_SLICES,
+        epsilon=epsilon,
+    )
+    j_mode = waveguide_3d.compute_source(
+        E=source_mode['E'],
+        wavenumber=source_mode['wavenumber'],
+        omega=OMEGA,
+        dxes=base_dxes,
+        axis=0,
+        polarity=1,
+        slices=SCATTERING_SOURCE_SLICES,
+        epsilon=epsilon,
+    )
+    reflected_mode = waveguide_3d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=base_dxes,
+        axis=0,
+        polarity=-1,
+        slices=SCATTERING_REFLECT_SLICES,
+        epsilon=epsilon,
+    )
+    reflected_overlap = waveguide_3d.compute_overlap_e(
+        E=reflected_mode['E'],
+        wavenumber=reflected_mode['wavenumber'],
+        dxes=base_dxes,
+        axis=0,
+        polarity=-1,
+        slices=SCATTERING_REFLECT_SLICES,
+        omega=OMEGA,
+    )
+    transmitted_mode = waveguide_3d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=base_dxes,
+        axis=0,
+        polarity=1,
+        slices=SCATTERING_TRANSMIT_SLICES,
+        epsilon=epsilon,
+    )
+    transmitted_overlap = waveguide_3d.compute_overlap_e(
+        E=transmitted_mode['E'],
+        wavenumber=transmitted_mode['wavenumber'],
+        dxes=base_dxes,
+        axis=0,
+        polarity=1,
+        slices=SCATTERING_TRANSMIT_SLICES,
+        omega=OMEGA,
+    )
+
+    update_e, update_h = fdtd.updates_with_cpml(cpml_params=_build_cpml_params(), dt=DT, dxes=base_dxes, epsilon=epsilon)
+
+    e_field = numpy.zeros_like(epsilon)
+    h_field = numpy.zeros_like(epsilon)
+    e_accumulator = numpy.zeros((1, *SCATTERING_SHAPE), dtype=complex)
+    h_accumulator = numpy.zeros((1, *SCATTERING_SHAPE), dtype=complex)
+    j_accumulator = numpy.zeros((1, *SCATTERING_SHAPE), dtype=complex)
+
+    warmup_steps = SCATTERING_WARMUP_PERIODS * PERIOD_STEPS
+    accumulation_steps = SCATTERING_ACCUMULATION_PERIODS * PERIOD_STEPS
+    for step in range(warmup_steps + accumulation_steps):
+        update_e(e_field, h_field, epsilon)
+
+        t_half = (step + 0.5) * DT
+        j_real = (j_mode.real * numpy.cos(OMEGA * t_half) - j_mode.imag * numpy.sin(OMEGA * t_half)).real
+        e_field -= DT * j_real / epsilon
+
+        if step >= warmup_steps:
+            fdtd.accumulate_phasor_j(j_accumulator, OMEGA, DT, j_real, step)
+            fdtd.accumulate_phasor_e(e_accumulator, OMEGA, DT, e_field, step + 1)
+
+        update_h(e_field, h_field)
+
+        if step >= warmup_steps:
+            fdtd.accumulate_phasor_h(h_accumulator, OMEGA, DT, h_field, step + 1)
+
+    e_ph = e_accumulator[0]
+    h_ph = h_accumulator[0]
+    j_ph = j_accumulator[0]
+
+    e_fdfd = unvec(
+        fdfd.solvers.generic(
+            J=vec(j_ph),
+            omega=OMEGA,
+            dxes=stretched_dxes,
+            epsilon=vec(epsilon),
+            matrix_solver_opts={'atol': 1e-10, 'rtol': 1e-7},
+        ),
+        SCATTERING_SHAPE[1:],
+    )
+    h_fdfd = functional.e2h(OMEGA, stretched_dxes)(e_fdfd)
+
+    reflected_td = vec(e_ph) @ vec(reflected_overlap).conj()
+    reflected_fd = vec(e_fdfd) @ vec(reflected_overlap).conj()
+    transmitted_td = vec(e_ph) @ vec(transmitted_overlap).conj()
+    transmitted_fd = vec(e_fdfd) @ vec(transmitted_overlap).conj()
+
+    poynting_td = functional.poynting_e_cross_h(stretched_dxes)(e_ph, h_ph.conj())
+    poynting_fd = functional.poynting_e_cross_h(stretched_dxes)(e_fdfd, h_fdfd.conj())
+    reflected_flux_td = float(0.5 * poynting_td[0, SCATTERING_REFLECT_SLICES[0], :, :].real.sum())
+    reflected_flux_fd = float(0.5 * poynting_fd[0, SCATTERING_REFLECT_SLICES[0], :, :].real.sum())
+    transmitted_flux_td = float(0.5 * poynting_td[0, SCATTERING_TRANSMIT_SLICES[0], :, :].real.sum())
+    transmitted_flux_fd = float(0.5 * poynting_fd[0, SCATTERING_TRANSMIT_SLICES[0], :, :].real.sum())
+
+    return WaveguideScatteringResult(
+        e_ph=e_ph,
+        h_ph=h_ph,
+        j_ph=j_ph,
+        e_fdfd=e_fdfd,
+        h_fdfd=h_fdfd,
+        reflected_td=reflected_td,
+        reflected_fd=reflected_fd,
+        transmitted_td=transmitted_td,
+        transmitted_fd=transmitted_fd,
+        reflected_flux_td=reflected_flux_td,
+        reflected_flux_fd=reflected_flux_fd,
+        transmitted_flux_td=transmitted_flux_td,
+        transmitted_flux_fd=transmitted_flux_fd,
+    )
+
+
+def test_straight_waveguide_base_variant_outperforms_stretched_variant() -> None:
+    base_result = _run_straight_waveguide_case('base')
+    stretched_result = _run_straight_waveguide_case('stretched')
+
+    assert base_result.variant == CHOSEN_VARIANT
+    assert base_result.combined_error < stretched_result.combined_error
+
+
+def test_straight_waveguide_fdtd_fdfd_overlap_and_flux_agree() -> None:
+    result = _run_straight_waveguide_case(CHOSEN_VARIANT)
+
+    assert numpy.isfinite(result.e_ph).all()
+    assert numpy.isfinite(result.h_ph).all()
+    assert numpy.isfinite(result.j_ph).all()
+    assert numpy.isfinite(result.e_fdfd).all()
+    assert numpy.isfinite(result.h_fdfd).all()
+    assert abs(result.overlap_td) > 0
+    assert abs(result.overlap_fd) > 0
+    assert abs(result.flux_td) > 0
+    assert abs(result.flux_fd) > 0
+
+    assert result.overlap_mag_rel_err < 0.01
+    assert result.flux_rel_err < 0.01
+    assert result.overlap_rel_err < 0.01
+    assert result.overlap_phase_deg < 0.5
+
+
+def test_width_step_waveguide_fdtd_fdfd_modal_powers_and_flux_agree() -> None:
+    result = _run_width_step_scattering_case()
+
+    assert numpy.isfinite(result.e_ph).all()
+    assert numpy.isfinite(result.h_ph).all()
+    assert numpy.isfinite(result.j_ph).all()
+    assert numpy.isfinite(result.e_fdfd).all()
+    assert numpy.isfinite(result.h_fdfd).all()
+    assert abs(result.reflected_td) > 0
+    assert abs(result.reflected_fd) > 0
+    assert abs(result.transmitted_td) > 0
+    assert abs(result.transmitted_fd) > 0
+    assert abs(result.reflected_flux_td) > 0
+    assert abs(result.reflected_flux_fd) > 0
+    assert abs(result.transmitted_flux_td) > 0
+    assert abs(result.transmitted_flux_fd) > 0
+
+    assert result.transmitted_overlap_mag_rel_err < 0.03
+    assert result.reflected_overlap_mag_rel_err < 0.03
+    assert result.transmitted_flux_rel_err < 0.01
+    assert result.reflected_flux_rel_err < 0.01

meanas/test/test_waveguide_mode_helpers.py

@@ -0,0 +1,299 @@
+import contextlib
+import io
+import numpy
+from numpy.linalg import norm
+import pytest
+import warnings
+
+from ..fdmath import vec, unvec
+from ..fdfd import waveguide_2d, waveguide_3d, waveguide_cyl
+
+
+OMEGA = 1 / 1500
+
+
+def build_waveguide_3d_mode(
+        *,
+        slice_start: int,
+        polarity: int,
+        ) -> tuple[numpy.ndarray, list[list[numpy.ndarray]], tuple[slice, slice, slice], dict[str, complex | numpy.ndarray]]:
+    epsilon = numpy.ones((3, 5, 5, 1), dtype=float)
+    dxes = [[numpy.ones(5), numpy.ones(5), numpy.ones(1)] for _ in range(2)]
+    slices = (slice(slice_start, slice_start + 1), slice(None), slice(None))
+    result = waveguide_3d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=dxes,
+        axis=0,
+        polarity=polarity,
+        slices=slices,
+        epsilon=epsilon,
+        )
+    return epsilon, dxes, slices, result
+
+
+def build_waveguide_cyl_fixture(
+        *,
+        nonuniform: bool = False,
+        ) -> tuple[list[list[numpy.ndarray]], numpy.ndarray, float]:
+    if nonuniform:
+        dxes = [
+            [numpy.array([1.0, 1.5, 1.2, 0.8, 1.1]), numpy.ones(5)],
+            [numpy.array([0.9, 1.4, 1.0, 0.7, 1.2]), numpy.ones(5)],
+        ]
+    else:
+        dxes = [[numpy.ones(5), numpy.ones(5)] for _ in range(2)]
+    epsilon = vec(numpy.ones((3, 5, 5), dtype=float))
+    return dxes, epsilon, 10.0
+
+
+def test_waveguide_3d_solve_mode_and_expand_e_are_phase_consistent() -> None:
+    epsilon, dxes, slices, result = build_waveguide_3d_mode(slice_start=0, polarity=1)
+    axis = 0
+    polarity = 1
+    expanded = waveguide_3d.expand_e(
+        E=result['E'],
+        wavenumber=result['wavenumber'],
+        dxes=dxes,
+        axis=axis,
+        polarity=polarity,
+        slices=slices,
+        )
+
+    dx_prop = 0.5 * numpy.array([dx[2][slices[2]] for dx in dxes]).sum()
+    expected_wavenumber = 2 / dx_prop * numpy.arcsin(result['wavenumber_2d'] * dx_prop / 2)
+    solved_slice = (slice(None), *slices)
+
+    assert result['E'].shape == epsilon.shape
+    assert result['H'].shape == epsilon.shape
+    assert numpy.isfinite(result['E']).all()
+    assert numpy.isfinite(result['H']).all()
+    assert abs(result['wavenumber'] - expected_wavenumber) < 1e-12
+    assert numpy.allclose(expanded[solved_slice], result['E'][solved_slice])
+
+    component, _x, y_index, z_index = numpy.unravel_index(
+        numpy.abs(result['E']).argmax(),
+        result['E'].shape,
+        )
+    values = expanded[component, :, y_index, z_index]
+    ratios = values[1:] / values[:-1]
+    expected_ratio = numpy.exp(-1j * result['wavenumber'])
+
+    numpy.testing.assert_allclose(ratios, expected_ratio, rtol=1e-6, atol=1e-9)
+
+
+@pytest.mark.parametrize(
+    ('polarity', 'expected_range'),
+    [(1, (0, 1)), (-1, (3, 4))],
+    )
+def test_waveguide_3d_compute_overlap_e_uses_adjacent_window(
+        polarity: int,
+        expected_range: tuple[int, int],
+        ) -> None:
+    _epsilon, dxes, slices, result = build_waveguide_3d_mode(slice_start=2, polarity=polarity)
+
+    with warnings.catch_warnings(record=True) as caught:
+        overlap = waveguide_3d.compute_overlap_e(
+            E=result['E'],
+            wavenumber=result['wavenumber'],
+            dxes=dxes,
+            axis=0,
+            polarity=polarity,
+            slices=slices,
+            omega=OMEGA,
+            )
+
+    nonzero = numpy.argwhere(numpy.abs(overlap) > 0)
+
+    assert not caught
+    assert numpy.isfinite(overlap).all()
+    assert nonzero[:, 1].min() == expected_range[0]
+    assert nonzero[:, 1].max() == expected_range[1]
+
+
+@pytest.mark.parametrize(
+    ('polarity', 'slice_start', 'expected_index'),
+    [(1, 1, 0), (-1, 3, 4)],
+    )
+def test_waveguide_3d_compute_overlap_e_warns_when_window_is_clipped(
+        polarity: int,
+        slice_start: int,
+        expected_index: int,
+        ) -> None:
+    _epsilon, dxes, slices, result = build_waveguide_3d_mode(slice_start=slice_start, polarity=polarity)
+
+    with pytest.warns(RuntimeWarning, match='clipped'):
+        overlap = waveguide_3d.compute_overlap_e(
+            E=result['E'],
+            wavenumber=result['wavenumber'],
+            dxes=dxes,
+            axis=0,
+            polarity=polarity,
+            slices=slices,
+            omega=OMEGA,
+            )
+
+    nonzero = numpy.argwhere(numpy.abs(overlap) > 0)
+
+    assert numpy.isfinite(overlap).all()
+    assert nonzero[:, 1].min() == expected_index
+    assert nonzero[:, 1].max() == expected_index
+
+
+@pytest.mark.parametrize(
+    ('polarity', 'slice_start'),
+    [(1, 0), (-1, 4)],
+    )
+def test_waveguide_3d_compute_overlap_e_rejects_empty_overlap_window(
+        polarity: int,
+        slice_start: int,
+        ) -> None:
+    _epsilon, dxes, slices, result = build_waveguide_3d_mode(slice_start=slice_start, polarity=polarity)
+
+    with pytest.raises(ValueError, match='outside the domain'):
+        waveguide_3d.compute_overlap_e(
+            E=result['E'],
+            wavenumber=result['wavenumber'],
+            dxes=dxes,
+            axis=0,
+            polarity=polarity,
+            slices=slices,
+            omega=OMEGA,
+            )
+
+
+def test_waveguide_3d_compute_overlap_e_rejects_zero_support_window() -> None:
+    _epsilon, dxes, slices, result = build_waveguide_3d_mode(slice_start=2, polarity=1)
+
+    with pytest.raises(ValueError, match='no overlap field support'):
+        waveguide_3d.compute_overlap_e(
+            E=numpy.zeros_like(result['E']),
+            wavenumber=result['wavenumber'],
+            dxes=dxes,
+            axis=0,
+            polarity=1,
+            slices=slices,
+            omega=OMEGA,
+            )
+
+
+def test_waveguide_cyl_solved_modes_are_ordered_and_low_residual() -> None:
+    dxes, epsilon, rmin = build_waveguide_cyl_fixture()
+
+    e_xys, angular_wavenumbers = waveguide_cyl.solve_modes(
+        [0, 1],
+        omega=OMEGA,
+        dxes=dxes,
+        epsilon=epsilon,
+        rmin=rmin,
+        )
+    operator = waveguide_cyl.cylindrical_operator(OMEGA, dxes, epsilon, rmin=rmin)
+
+    assert numpy.all(numpy.diff(numpy.real(angular_wavenumbers)) <= 0)
+
+    for e_xy, angular_wavenumber in zip(e_xys, angular_wavenumbers, strict=True):
+        eigenvalue = (angular_wavenumber / rmin) ** 2
+        residual = norm(operator @ e_xy - eigenvalue * e_xy) / norm(e_xy)
+        assert residual < 1e-6
+
+
+def test_waveguide_cyl_linear_wavenumbers_are_finite_and_ordered() -> None:
+    dxes, epsilon, rmin = build_waveguide_cyl_fixture()
+
+    e_xys, angular_wavenumbers = waveguide_cyl.solve_modes(
+        [0, 1],
+        omega=OMEGA,
+        dxes=dxes,
+        epsilon=epsilon,
+        rmin=10.0,
+        )
+    linear_wavenumbers = waveguide_cyl.linear_wavenumbers(
+        e_xys,
+        angular_wavenumbers,
+        epsilon=epsilon,
+        dxes=dxes,
+        rmin=rmin,
+        )
+
+    assert numpy.isfinite(linear_wavenumbers).all()
+    assert numpy.all(numpy.real(linear_wavenumbers) > 0)
+    assert numpy.all(numpy.diff(numpy.real(linear_wavenumbers)) <= 0)
+
+
+def test_waveguide_cyl_dxes2t_matches_expected_radius_scaling() -> None:
+    dxes, _epsilon, rmin = build_waveguide_cyl_fixture(nonuniform=True)
+    Ta, Tb = waveguide_cyl.dxes2T(dxes, rmin)
+
+    ta = (rmin + numpy.cumsum(dxes[0][0])) / rmin
+    tb = (rmin + dxes[0][0] / 2 + numpy.cumsum(dxes[1][0])) / rmin
+
+    numpy.testing.assert_allclose(Ta.diagonal(), numpy.repeat(ta, dxes[0][1].size))
+    numpy.testing.assert_allclose(Tb.diagonal(), numpy.repeat(tb, dxes[1][1].size))
+
+
+def test_waveguide_cyl_exy2e_and_exy2h_return_finite_full_fields() -> None:
+    dxes, epsilon, rmin = build_waveguide_cyl_fixture()
+    mu = vec(2 * numpy.ones((3, 5, 5), dtype=float))
+    e_xy, angular_wavenumber = waveguide_cyl.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=dxes,
+        epsilon=epsilon,
+        rmin=rmin,
+        )
+
+    e_field = waveguide_cyl.exy2e(
+        angular_wavenumber=angular_wavenumber,
+        omega=OMEGA,
+        dxes=dxes,
+        rmin=rmin,
+        epsilon=epsilon,
+        ) @ e_xy
+    h_field = waveguide_cyl.exy2h(
+        angular_wavenumber=angular_wavenumber,
+        omega=OMEGA,
+        dxes=dxes,
+        rmin=rmin,
+        epsilon=epsilon,
+        mu=mu,
+        ) @ e_xy
+
+    assert e_field.shape == (3 * 25,)
+    assert h_field.shape == (3 * 25,)
+    assert numpy.isfinite(e_field).all()
+    assert numpy.isfinite(h_field).all()
+    assert unvec(e_field, (5, 5)).shape == (3, 5, 5)
+    assert unvec(h_field, (5, 5)).shape == (3, 5, 5)
+
+
+@pytest.mark.parametrize('use_mu', [False, True])
+def test_waveguide_cyl_normalized_fields_are_unit_norm_and_silent(use_mu: bool) -> None:
+    dxes, epsilon, rmin = build_waveguide_cyl_fixture()
+    mu = vec(2 * numpy.ones((3, 5, 5), dtype=float)) if use_mu else None
+    e_xy, angular_wavenumber = waveguide_cyl.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=dxes,
+        epsilon=epsilon,
+        rmin=rmin,
+        )
+
+    output = io.StringIO()
+    with contextlib.redirect_stdout(output):
+        e_field, h_field = waveguide_cyl.normalized_fields_e(
+            e_xy,
+            angular_wavenumber=angular_wavenumber,
+            omega=OMEGA,
+            dxes=dxes,
+            rmin=rmin,
+            epsilon=epsilon,
+            mu=mu,
+            )
+
+    overlap = waveguide_2d.inner_product(e_field, h_field, dxes, conj_h=True)
+
+    assert output.getvalue() == ''
+    assert numpy.isfinite(e_field).all()
+    assert numpy.isfinite(h_field).all()
+    assert abs(overlap.real - 1.0) < 1e-10
+    assert abs(overlap.imag) < 1e-10

meanas/test/utils.py

@@ -2,7 +2,8 @@ import numpy
 from numpy.typing import NDArray
 
 
-PRNG = numpy.random.RandomState(12345)
+def make_prng(seed: int = 12345) -> numpy.random.RandomState:
+    return numpy.random.RandomState(seed)
 
 
 def assert_fields_close(
@@ -29,4 +30,3 @@ def assert_close(
         **kwargs,
         ) -> None:
     numpy.testing.assert_allclose(x, y, *args, **kwargs)
-

mkdocs.yml

@@ -0,0 +1,76 @@
+site_name: meanas
+site_description: Electromagnetic simulation tools
+site_url: !ENV [DOCS_SITE_URL, ""]
+repo_url: https://mpxd.net/code/jan/meanas
+repo_name: meanas
+docs_dir: docs
+site_dir: site
+strict: false
+
+theme:
+  name: material
+  font: false
+  features:
+    - navigation.indexes
+    - navigation.sections
+    - navigation.top
+    - content.code.copy
+    - toc.follow
+
+nav:
+  - Home: index.md
+  - API:
+      - Overview: api/index.md
+      - meanas: api/meanas.md
+      - eigensolvers: api/eigensolvers.md
+      - fdfd: api/fdfd.md
+      - waveguides: api/waveguides.md
+      - fdtd: api/fdtd.md
+      - fdmath: api/fdmath.md
+
+plugins:
+  - search
+  - mkdocstrings:
+      handlers:
+        python:
+          paths:
+            - .
+          options:
+            show_root_heading: true
+            show_root_toc_entry: false
+            show_source: false
+            show_signature_annotations: true
+            show_symbol_type_heading: true
+            show_symbol_type_toc: true
+            members_order: source
+            separate_signature: true
+            merge_init_into_class: true
+            docstring_style: google
+  - print-site
+
+markdown_extensions:
+  - admonition
+  - attr_list
+  - md_in_html
+  - tables
+  - toc:
+      permalink: true
+  - pymdownx.arithmatex:
+      generic: true
+  - pymdownx.highlight:
+      anchor_linenums: true
+  - pymdownx.inlinehilite
+  - pymdownx.snippets
+  - pymdownx.superfences
+  - pymdownx.tabbed:
+      alternate_style: true
+
+extra_css:
+  - stylesheets/extra.css
+
+extra_javascript:
+  - javascripts/mathjax.js
+  - assets/vendor/mathjax/startup.js
+
+watch:
+  - meanas

pdoc_templates/config.mako

@@ -1,47 +0,0 @@
-<%!
-    # Template configuration. Copy over in your template directory
-    # (used with --template-dir) and adapt as required.
-    html_lang = 'en'
-    show_inherited_members = False
-    extract_module_toc_into_sidebar = True
-    list_class_variables_in_index = True
-    sort_identifiers = True
-    show_type_annotations = True
-
-    # Show collapsed source code block next to each item.
-    # Disabling this can improve rendering speed of large modules.
-    show_source_code = True
-
-    # If set, format links to objects in online source code repository
-    # according to this template. Supported keywords for interpolation
-    # are: commit, path, start_line, end_line.
-    #git_link_template = 'https://github.com/USER/PROJECT/blob/{commit}/{path}#L{start_line}-L{end_line}'
-    #git_link_template = 'https://gitlab.com/USER/PROJECT/blob/{commit}/{path}#L{start_line}-L{end_line}'
-    #git_link_template = 'https://bitbucket.org/USER/PROJECT/src/{commit}/{path}#lines-{start_line}:{end_line}'
-    #git_link_template = 'https://CGIT_HOSTNAME/PROJECT/tree/{path}?id={commit}#n{start_line}'
-    #git_link_template = None
-    git_link_template = 'https://mpxd.net/code/jan/meanas/src/commit/{commit}/{path}#L{start_line}-L{end_line}'
-
-    # A prefix to use for every HTML hyperlink in the generated documentation.
-    # No prefix results in all links being relative.
-    link_prefix = ''
-
-    # Enable syntax highlighting for code/source blocks by including Highlight.js
-    syntax_highlighting = True
-
-    # Set the style keyword such as 'atom-one-light' or 'github-gist'
-    #     Options: https://github.com/highlightjs/highlight.js/tree/master/src/styles
-    #     Demo: https://highlightjs.org/static/demo/
-    hljs_style = 'github'
-
-    # If set, insert Google Analytics tracking code. Value is GA
-    # tracking id (UA-XXXXXX-Y).
-    google_analytics = ''
-
-    # If set, render LaTeX math syntax within \(...\) (inline equations),
-    # or within \[...\] or $$...$$ or `.. math::` (block equations)
-    # as nicely-formatted math formulas using MathJax.
-    # Note: in Python docstrings, either all backslashes need to be escaped (\\)
-    # or you need to use raw r-strings.
-    latex_math = True
-%>

pdoc_templates/css.mako

@@ -1,389 +0,0 @@
-<%!
-    from pdoc.html_helpers import minify_css
-%>
-
-<%def name="mobile()" filter="minify_css">
-  .flex {
-    display: flex !important;
-  }
-
-  body {
-    line-height: 1.5em;
-    background: black;
-    color: #DDD;
-  }
-
-  #content {
-    padding: 20px;
-  }
-
-  #sidebar {
-    padding: 30px;
-    overflow: hidden;
-  }
-
-  .http-server-breadcrumbs {
-    font-size: 130%;
-    margin: 0 0 15px 0;
-  }
-
-  #footer {
-    font-size: .75em;
-    padding: 5px 30px;
-    border-top: 1px solid #ddd;
-    text-align: right;
-  }
-    #footer p {
-      margin: 0 0 0 1em;
-      display: inline-block;
-    }
-    #footer p:last-child {
-      margin-right: 30px;
-    }
-
-  h1, h2, h3, h4, h5 {
-    font-weight: 300;
-  }
-  h1 {
-    font-size: 2.5em;
-    line-height: 1.1em;
-  }
-  h2 {
-    font-size: 1.75em;
-    margin: 1em 0 .50em 0;
-  }
-  h3 {
-    font-size: 1.4em;
-    margin: 25px 0 10px 0;
-  }
-  h4 {
-    margin: 0;
-    font-size: 105%;
-  }
-
-  a {
-    color: #999;
-    text-decoration: none;
-    transition: color .3s ease-in-out;
-  }
-  a:hover {
-    color: #18d;
-  }
-
-  .title code {
-    font-weight: bold;
-  }
-  h2[id^="header-"] {
-    margin-top: 2em;
-  }
-  .ident {
-    color: #7ff;
-  }
-
-  pre code {
-    background: transparent;
-    font-size: .8em;
-    line-height: 1.4em;
-  }
-  code {
-    background: #0d0d0e;
-    padding: 1px 4px;
-    overflow-wrap: break-word;
-  }
-  h1 code { background: transparent }
-
-  pre {
-    background: #111;
-    border: 0;
-    border-top: 1px solid #ccc;
-    border-bottom: 1px solid #ccc;
-    margin: 1em 0;
-    padding: 1ex;
-  }
-
-  #http-server-module-list {
-    display: flex;
-    flex-flow: column;
-  }
-    #http-server-module-list div {
-      display: flex;
-    }
-    #http-server-module-list dt {
-      min-width: 10%;
-    }
-    #http-server-module-list p {
-      margin-top: 0;
-    }
-
-  .toc ul,
-  #index {
-    list-style-type: none;
-    margin: 0;
-    padding: 0;
-  }
-    #index code {
-      background: transparent;
-    }
-    #index h3 {
-      border-bottom: 1px solid #ddd;
-    }
-    #index ul {
-      padding: 0;
-    }
-    #index h4 {
-      font-weight: bold;
-    }
-    #index h4 + ul {
-      margin-bottom:.6em;
-    }
-    /* Make TOC lists have 2+ columns when viewport is wide enough.
-       Assuming ~20-character identifiers and ~30% wide sidebar. */
-    @media (min-width: 200ex) { #index .two-column { column-count: 2 } }
-    @media (min-width: 300ex) { #index .two-column { column-count: 3 } }
-
-  dl {
-    margin-bottom: 2em;
-  }
-    dl dl:last-child {
-      margin-bottom: 4em;
-    }
-  dd {
-    margin: 0 0 1em 3em;
-  }
-    #header-classes + dl > dd {
-      margin-bottom: 3em;
-    }
-    dd dd {
-      margin-left: 2em;
-    }
-    dd p {
-      margin: 10px 0;
-    }
-    .name {
-      background: #111;
-      font-weight: bold;
-      font-size: .85em;
-      padding: 5px 10px;
-      display: inline-block;
-      min-width: 40%;
-    }
-      .name:hover {
-        background: #101010;
-      }
-      .name > span:first-child {
-        white-space: nowrap;
-      }
-      .name.class > span:nth-child(2) {
-        margin-left: .4em;
-      }
-    .inherited {
-      color: #777;
-      border-left: 5px solid #eee;
-      padding-left: 1em;
-    }
-    .inheritance em {
-      font-style: normal;
-      font-weight: bold;
-    }
-
-    /* Docstrings titles, e.g. in numpydoc format */
-    .desc h2 {
-      font-weight: 400;
-      font-size: 1.25em;
-    }
-    .desc h3 {
-      font-size: 1em;
-    }
-    .desc dt code {
-      background: inherit;  /* Don't grey-back parameters */
-    }
-
-    .source summary,
-    .git-link-div {
-      color: #aaa;
-      text-align: right;
-      font-weight: 400;
-      font-size: .8em;
-      text-transform: uppercase;
-    }
-      .source summary > * {
-        white-space: nowrap;
-        cursor: pointer;
-      }
-      .git-link {
-        color: inherit;
-        margin-left: 1em;
-      }
-    .source pre {
-      max-height: 500px;
-      overflow: auto;
-      margin: 0;
-    }
-    .source pre code {
-      font-size: 12px;
-      overflow: visible;
-    }
-  .hlist {
-    list-style: none;
-  }
-    .hlist li {
-      display: inline;
-    }
-    .hlist li:after {
-      content: ',\2002';
-    }
-    .hlist li:last-child:after {
-      content: none;
-    }
-    .hlist .hlist {
-      display: inline;
-      padding-left: 1em;
-    }
-
-  img {
-    max-width: 100%;
-  }
-
-  .admonition {
-    padding: .1em .5em;
-    margin-bottom: 1em;
-  }
-    .admonition-title {
-      font-weight: bold;
-    }
-    .admonition.note,
-    .admonition.info,
-    .admonition.important {
-      background: #610;
-    }
-    .admonition.todo,
-    .admonition.versionadded,
-    .admonition.tip,
-    .admonition.hint {
-      background: #202;
-    }
-    .admonition.warning,
-    .admonition.versionchanged,
-    .admonition.deprecated {
-      background: #02b;
-    }
-    .admonition.error,
-    .admonition.danger,
-    .admonition.caution {
-      background: darkpink;
-    }
-</%def>
-
-<%def name="desktop()" filter="minify_css">
-  @media screen and (min-width: 700px) {
-    #sidebar {
-      width: 30%;
-    }
-    #content {
-      width: 70%;
-      max-width: 100ch;
-      padding: 3em 4em;
-      border-left: 1px solid #ddd;
-    }
-    pre code {
-      font-size: 1em;
-    }
-    .item .name {
-      font-size: 1em;
-    }
-    main {
-      display: flex;
-      flex-direction: row-reverse;
-      justify-content: flex-end;
-    }
-    .toc ul ul,
-    #index ul {
-      padding-left: 1.5em;
-    }
-    .toc > ul > li {
-      margin-top: .5em;
-    }
-  }
-</%def>
-
-<%def name="print()" filter="minify_css">
-@media print {
-  #sidebar h1 {
-    page-break-before: always;
-  }
-  .source {
-    display: none;
-  }
-}
-@media print {
-    * {
-        background: transparent !important;
-        color: #000 !important; /* Black prints faster: h5bp.com/s */
-        box-shadow: none !important;
-        text-shadow: none !important;
-    }
-
-    a[href]:after {
-        content: " (" attr(href) ")";
-        font-size: 90%;
-    }
-    /* Internal, documentation links, recognized by having a title,
-       don't need the URL explicity stated. */
-    a[href][title]:after {
-        content: none;
-    }
-
-    abbr[title]:after {
-        content: " (" attr(title) ")";
-    }
-
-    /*
-     * Don't show links for images, or javascript/internal links
-     */
-
-    .ir a:after,
-    a[href^="javascript:"]:after,
-    a[href^="#"]:after {
-        content: "";
-    }
-
-    pre,
-    blockquote {
-        border: 1px solid #999;
-        page-break-inside: avoid;
-    }
-
-    thead {
-        display: table-header-group; /* h5bp.com/t */
-    }
-
-    tr,
-    img {
-        page-break-inside: avoid;
-    }
-
-    img {
-        max-width: 100% !important;
-    }
-
-    @page {
-        margin: 0.5cm;
-    }
-
-    p,
-    h2,
-    h3 {
-        orphans: 3;
-        widows: 3;
-    }
-
-    h1,
-    h2,
-    h3,
-    h4,
-    h5,
-    h6 {
-        page-break-after: avoid;
-    }
-}
-</%def>

pdoc_templates/html.mako

@@ -1,445 +0,0 @@
-<%
-  import os
-
-  import pdoc
-  from pdoc.html_helpers import extract_toc, glimpse, to_html as _to_html, format_git_link, _md, to_markdown
-
-  from markdown.inlinepatterns import InlineProcessor
-  from markdown.util import AtomicString, etree
-
-
-  def link(d, name=None, fmt='{}'):
-    name = fmt.format(name or d.qualname + ('()' if isinstance(d, pdoc.Function) else ''))
-    if not isinstance(d, pdoc.Doc) or isinstance(d, pdoc.External) and not external_links:
-        return name
-    url = d.url(relative_to=module, link_prefix=link_prefix,
-                top_ancestor=not show_inherited_members)
-    return '<a title="{}" href="{}">{}</a>'.format(d.refname, url, name)
-
-
-  # Altered latex delimeters (allow inline $...$, wrap in <eq></eq>)
-  class _MathPattern(InlineProcessor):
-      NAME = 'pdoc-math'
-      PATTERN = r'(?<!\S|\\)(?:\\\((.+?)\\\)|\\\[(.+?)\\\]|\$\$(.+?)\$\$|\$(\S.*?)\$)'
-      PRIORITY = 181  # Larger than that of 'escape' pattern
-
-      def handleMatch(self, m, data):
-          for value, is_block in zip(m.groups(), (False, True, True, False)):
-              if value:
-                  break
-          wrapper = etree.Element('eq')
-          wrapper.text = AtomicString(value)
-          return wrapper, m.start(0), m.end(0)
-
-  def to_html(text: str):
-      if not latex_math and _MathPattern.NAME in _md.inlinePatterns:
-          _md.inlinePatterns.deregister(_MathPattern.NAME)
-      elif latex_math and _MathPattern.NAME not in _md.inlinePatterns:
-          _md.inlinePatterns.register(_MathPattern(_MathPattern.PATTERN),
-                                      _MathPattern.NAME,
-                                      _MathPattern.PRIORITY)
-      md = to_markdown(text, docformat='numpy,google', module=module, link=link)
-      return _md.reset().convert(md)
-
-
-#  def to_html(text):
-#    return _to_html(text, module=module, link=link, latex_math=latex_math)
-%>
-
-<%def name="ident(name)"><span class="ident">${name}</span></%def>
-
-<%def name="show_source(d)">
-  % if (show_source_code or git_link_template) and d.source and d.obj is not getattr(d.inherits, 'obj', None):
-    <% git_link = format_git_link(git_link_template, d) %>
-    % if show_source_code:
-      <details class="source">
-        <summary>
-            <span>Expand source code</span>
-            % if git_link:
-              <a href="${git_link}" class="git-link">Browse git</a>
-            %endif
-        </summary>
-        <pre><code class="python">${d.source | h}</code></pre>
-      </details>
-    % elif git_link:
-      <div class="git-link-div"><a href="${git_link}" class="git-link">Browse git</a></div>
-    %endif
-  %endif
-</%def>
-
-<%def name="show_desc(d, short=False)">
-  <%
-  inherits = ' inherited' if d.inherits else ''
-  docstring = glimpse(d.docstring) if short or inherits else d.docstring
-  %>
-  % if d.inherits:
-      <p class="inheritance">
-          <em>Inherited from:</em>
-          % if hasattr(d.inherits, 'cls'):
-              <code>${link(d.inherits.cls)}</code>.<code>${link(d.inherits, d.name)}</code>
-          % else:
-              <code>${link(d.inherits)}</code>
-          % endif
-      </p>
-  % endif
-  <section class="desc${inherits}">${docstring | to_html}</section>
-  % if not isinstance(d, pdoc.Module):
-  ${show_source(d)}
-  % endif
-</%def>
-
-<%def name="show_module_list(modules)">
-<h1>Python module list</h1>
-
-% if not modules:
-  <p>No modules found.</p>
-% else:
-  <dl id="http-server-module-list">
-  % for name, desc in modules:
-      <div class="flex">
-      <dt><a href="${link_prefix}${name}">${name}</a></dt>
-      <dd>${desc | glimpse, to_html}</dd>
-      </div>
-  % endfor
-  </dl>
-% endif
-</%def>
-
-<%def name="show_column_list(items)">
-  <%
-      two_column = len(items) >= 6 and all(len(i.name) < 20 for i in items)
-  %>
-  <ul class="${'two-column' if two_column else ''}">
-  % for item in items:
-    <li><code>${link(item, item.name)}</code></li>
-  % endfor
-  </ul>
-</%def>
-
-<%def name="show_module(module)">
-  <%
-  variables = module.variables(sort=sort_identifiers)
-  classes = module.classes(sort=sort_identifiers)
-  functions = module.functions(sort=sort_identifiers)
-  submodules = module.submodules()
-  %>
-
-  <%def name="show_func(f)">
-    <dt id="${f.refname}"><code class="name flex">
-        <%
-            params = ', '.join(f.params(annotate=show_type_annotations, link=link))
-            returns = show_type_annotations and f.return_annotation(link=link) or ''
-            if returns:
-                returns = ' ->\N{NBSP}' + returns
-        %>
-        <span>${f.funcdef()} ${ident(f.name)}</span>(<span>${params})${returns}</span>
-    </code></dt>
-    <dd>${show_desc(f)}</dd>
-  </%def>
-
-  <header>
-  % if http_server:
-    <nav class="http-server-breadcrumbs">
-      <a href="/">All packages</a>
-      <% parts = module.name.split('.')[:-1] %>
-      % for i, m in enumerate(parts):
-        <% parent = '.'.join(parts[:i+1]) %>
-        :: <a href="/${parent.replace('.', '/')}/">${parent}</a>
-      % endfor
-    </nav>
-  % endif
-  <h1 class="title">${'Namespace' if module.is_namespace else 'Module'} <code>${module.name}</code></h1>
-  </header>
-
-  <section id="section-intro">
-  ${module.docstring | to_html}
-  ${show_source(module)}
-  </section>
-
-  <section>
-    % if submodules:
-    <h2 class="section-title" id="header-submodules">Sub-modules</h2>
-    <dl>
-    % for m in submodules:
-      <dt><code class="name">${link(m)}</code></dt>
-      <dd>${show_desc(m, short=True)}</dd>
-    % endfor
-    </dl>
-    % endif
-  </section>
-
-  <section>
-    % if variables:
-    <h2 class="section-title" id="header-variables">Global variables</h2>
-    <dl>
-    % for v in variables:
-      <dt id="${v.refname}"><code class="name">var ${ident(v.name)}</code></dt>
-      <dd>${show_desc(v)}</dd>
-    % endfor
-    </dl>
-    % endif
-  </section>
-
-  <section>
-    % if functions:
-    <h2 class="section-title" id="header-functions">Functions</h2>
-    <dl>
-    % for f in functions:
-      ${show_func(f)}
-    % endfor
-    </dl>
-    % endif
-  </section>
-
-  <section>
-    % if classes:
-    <h2 class="section-title" id="header-classes">Classes</h2>
-    <dl>
-    % for c in classes:
-      <%
-      class_vars = c.class_variables(show_inherited_members, sort=sort_identifiers)
-      smethods = c.functions(show_inherited_members, sort=sort_identifiers)
-      inst_vars = c.instance_variables(show_inherited_members, sort=sort_identifiers)
-      methods = c.methods(show_inherited_members, sort=sort_identifiers)
-      mro = c.mro()
-      subclasses = c.subclasses()
-      params = ', '.join(c.params(annotate=show_type_annotations, link=link))
-      %>
-      <dt id="${c.refname}"><code class="flex name class">
-          <span>class ${ident(c.name)}</span>
-          % if params:
-              <span>(</span><span>${params})</span>
-          % endif
-      </code></dt>
-
-      <dd>${show_desc(c)}
-
-      % if mro:
-          <h3>Ancestors</h3>
-          <ul class="hlist">
-          % for cls in mro:
-              <li>${link(cls)}</li>
-          % endfor
-          </ul>
-      %endif
-
-      % if subclasses:
-          <h3>Subclasses</h3>
-          <ul class="hlist">
-          % for sub in subclasses:
-              <li>${link(sub)}</li>
-          % endfor
-          </ul>
-      % endif
-      % if class_vars:
-          <h3>Class variables</h3>
-          <dl>
-          % for v in class_vars:
-              <dt id="${v.refname}"><code class="name">var ${ident(v.name)}</code></dt>
-              <dd>${show_desc(v)}</dd>
-          % endfor
-          </dl>
-      % endif
-      % if smethods:
-          <h3>Static methods</h3>
-          <dl>
-          % for f in smethods:
-              ${show_func(f)}
-          % endfor
-          </dl>
-      % endif
-      % if inst_vars:
-          <h3>Instance variables</h3>
-          <dl>
-          % for v in inst_vars:
-              <dt id="${v.refname}"><code class="name">var ${ident(v.name)}</code></dt>
-              <dd>${show_desc(v)}</dd>
-          % endfor
-          </dl>
-      % endif
-      % if methods:
-          <h3>Methods</h3>
-          <dl>
-          % for f in methods:
-              ${show_func(f)}
-          % endfor
-          </dl>
-      % endif
-
-      % if not show_inherited_members:
-          <%
-              members = c.inherited_members()
-          %>
-          % if members:
-              <h3>Inherited members</h3>
-              <ul class="hlist">
-              % for cls, mems in members:
-                  <li><code><b>${link(cls)}</b></code>:
-                      <ul class="hlist">
-                          % for m in mems:
-                              <li><code>${link(m, name=m.name)}</code></li>
-                          % endfor
-                      </ul>
-
-                  </li>
-              % endfor
-              </ul>
-          % endif
-      % endif
-
-      </dd>
-    % endfor
-    </dl>
-    % endif
-  </section>
-</%def>
-
-<%def name="module_index(module)">
-  <%
-  variables = module.variables(sort=sort_identifiers)
-  classes = module.classes(sort=sort_identifiers)
-  functions = module.functions(sort=sort_identifiers)
-  submodules = module.submodules()
-  supermodule = module.supermodule
-  %>
-  <nav id="sidebar">
-
-    <%include file="logo.mako"/>
-
-    <h1>Index</h1>
-    ${extract_toc(module.docstring) if extract_module_toc_into_sidebar else ''}
-    <ul id="index">
-    % if supermodule:
-    <li><h3>Super-module</h3>
-      <ul>
-        <li><code>${link(supermodule)}</code></li>
-      </ul>
-    </li>
-    % endif
-
-    % if submodules:
-    <li><h3><a href="#header-submodules">Sub-modules</a></h3>
-      <ul>
-      % for m in submodules:
-        <li><code>${link(m)}</code></li>
-      % endfor
-      </ul>
-    </li>
-    % endif
-
-    % if variables:
-    <li><h3><a href="#header-variables">Global variables</a></h3>
-      ${show_column_list(variables)}
-    </li>
-    % endif
-
-    % if functions:
-    <li><h3><a href="#header-functions">Functions</a></h3>
-      ${show_column_list(functions)}
-    </li>
-    % endif
-
-    % if classes:
-    <li><h3><a href="#header-classes">Classes</a></h3>
-      <ul>
-      % for c in classes:
-        <li>
-        <h4><code>${link(c)}</code></h4>
-        <%
-            members = c.functions(sort=sort_identifiers) + c.methods(sort=sort_identifiers)
-            if list_class_variables_in_index:
-                members += (c.instance_variables(sort=sort_identifiers) +
-                            c.class_variables(sort=sort_identifiers))
-            if not show_inherited_members:
-                members = [i for i in members if not i.inherits]
-            if sort_identifiers:
-              members = sorted(members)
-        %>
-        % if members:
-          ${show_column_list(members)}
-        % endif
-        </li>
-      % endfor
-      </ul>
-    </li>
-    % endif
-
-    </ul>
-  </nav>
-</%def>
-
-<!doctype html>
-<html lang="${html_lang}">
-<head>
-  <meta charset="utf-8">
-  <meta name="viewport" content="width=device-width, initial-scale=1, minimum-scale=1" />
-  <meta name="generator" content="pdoc ${pdoc.__version__}" />
-
-<%
-    module_list = 'modules' in context.keys()  # Whether we're showing module list in server mode
-%>
-
-  % if module_list:
-    <title>Python module list</title>
-    <meta name="description" content="A list of documented Python modules." />
-  % else:
-    <title>${module.name} API documentation</title>
-    <meta name="description" content="${module.docstring | glimpse, trim, h}" />
-  % endif
-
-    <link href='https://mpxd.net/scripts/normalize.css/normalize.css' rel='stylesheet'>
-    <link href='https://mpxd.net/scripts/sanitize.css/sanitize.css' rel='stylesheet'>
-  % if syntax_highlighting:
-    <link href="https://mpxd.net/scripts/highlightjs/styles/${hljs_style}.min.css" rel="stylesheet">
-  %endif
-
-  <%namespace name="css" file="css.mako" />
-  <style>${css.mobile()}</style>
-  <style media="screen and (min-width: 700px)">${css.desktop()}</style>
-  <style media="print">${css.print()}</style>
-
-  % if google_analytics:
-    <script>
-    window.ga=window.ga||function(){(ga.q=ga.q||[]).push(arguments)};ga.l=+new Date;
-    ga('create', '${google_analytics}', 'auto'); ga('send', 'pageview');
-    </script><script async src='https://www.google-analytics.com/analytics.js'></script>
-  % endif
-
-  <%include file="head.mako"/>
-</head>
-<body>
-<main>
-  % if module_list:
-    <article id="content">
-      ${show_module_list(modules)}
-    </article>
-  % else:
-    <article id="content">
-      ${show_module(module)}
-    </article>
-    ${module_index(module)}
-  % endif
-</main>
-
-<footer id="footer">
-    <%include file="credits.mako"/>
-    <p>Generated by <a href="https://pdoc3.github.io/pdoc"><cite>pdoc</cite> ${pdoc.__version__}</a>.</p>
-</footer>
-
-% if syntax_highlighting:
-    <script src="https://mpxd.net/scripts/highlightjs/highlight.pack.js"></script>
-    <script>hljs.initHighlightingOnLoad()</script>
-% endif
-
-% if http_server and module:  ## Auto-reload on file change in dev mode
-    <script>
-    setInterval(() =>
-        fetch(window.location.href, {
-            method: "HEAD",
-            cache: "no-store",
-            headers: {"If-None-Match": "${os.stat(module.obj.__file__).st_mtime}"},
-        }).then(response => response.ok && window.location.reload()), 700);
-    </script>
-% endif
-</body>
-</html>