.forgejo/workflows/docs.yml

@@ -1,57 +0,0 @@
-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,10 +54,6 @@ coverage.xml
 
 # documentation
 doc/
-site/
-_doc_mathimg/
-doc.md
-doc.htex
 
 # PyBuilder
 target/

README.md

@@ -94,99 +94,6 @@ python3 -m pytest -rsxX | tee test_results.txt
 
 ## Use
 
-`meanas` is organized around a few core workflows:
-
-- `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.
+See `examples/` for some simple examples; you may need additional
+packages such as [gridlock](https://mpxd.net/code/jan/gridlock)
+to run the examples.

docs/api/eigensolvers.md

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

docs/api/fdfd.md

@@ -1,15 +0,0 @@
-# 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

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

docs/api/fdtd.md

@@ -1,15 +0,0 @@
-# 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

@@ -1,14 +0,0 @@
-# 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

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

docs/api/waveguides.md

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

docs/assets/vendor/mathjax/manifest.json

@@ -1,38 +0,0 @@
-{
-  "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

@@ -1,33 +0,0 @@
-# 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

@@ -1,19 +0,0 @@
-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

@@ -1,13 +0,0 @@
-.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,5 +1,4 @@
 import importlib
-import logging
 import numpy
 from numpy.linalg import norm
 from matplotlib import pyplot, colors
@@ -7,14 +6,12 @@ import logging
 
 import meanas
 from meanas import fdtd
-from meanas.fdmath import vec, unvec, fdfield_t
+from meanas.fdmath import vec, unvec
 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,30 +1,18 @@
 """
-Example code for a broadband FDTD run with phasor extraction.
+Example code for running an OpenCL FDTD simulation
 
-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.
+See main() for simulation setup.
 """
 
 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 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
+from masque import Pattern, shapes
 import gridlock
 import pcgen
 
@@ -53,7 +41,8 @@ def perturbed_l3(a: float, radius: float, **kwargs) -> Pattern:
         `masque.Pattern` object containing the L3 design
     """
 
-    default_args = {
+    default_args = {'hole_dose':    1,
+                    'trench_dose':  1,
                     'hole_layer':   0,
                     'trench_layer': 1,
                     'shifts_a':     (0.15, 0, 0.075),
@@ -64,40 +53,38 @@ def perturbed_l3(a: float, radius: float, **kwargs) -> Pattern:
                     }
     kwargs = {**default_args, **kwargs}
 
-    xyr = pcgen.l3_shift_perturbed_defect(
-        mirror_dims=kwargs['xy_size'],
+    xyr = pcgen.l3_shift_perturbed_defect(mirror_dims=kwargs['xy_size'],
                                           perturbed_radius=kwargs['perturbed_radius'],
                                           shifts_a=kwargs['shifts_a'],
-        shifts_r=kwargs['shifts_r'],
-        )
+                                          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.shapes[(kwargs['hole_layer'], 0)] += [
-        Circle(radius=r, offset=(x, y))
+    pat.name = f'L3p-a{a:g}r{radius:g}rp{kwargs["perturbed_radius"]:g}'
+    pat.shapes += [shapes.Circle(radius=r, offset=(x, y),
+                                 dose=kwargs['hole_dose'],
+                                 layer=kwargs['hole_layer'])
                    for x, y, r in xyr]
 
     maxes = numpy.max(numpy.fabs(xyr), axis=0)
-    pat.shapes[(kwargs['trench_layer'], 0)] += [
-        Polygon.rectangle(
+    pat.shapes += [shapes.Polygon.rectangle(
         lx=(2 * maxes[0]), ly=kwargs['trench_width'],
-            offset=(0, s * (maxes[1] + a + kwargs['trench_width'] / 2))
-            )
+        offset=(0, s * (maxes[1] + a + kwargs['trench_width'] / 2)),
+        dose=kwargs['trench_dose'], layer=kwargs['trench_layer'])
                    for s in (-1, 1)]
     return pat
 
 
 def main():
     dtype = numpy.float32
-    max_t = 3600            # number of timesteps
+    max_t = 8000            # number of timesteps
 
     dx = 40                 # discretization (nm/cell)
     pml_thickness = 8       # (number of cells)
 
     wl = 1550               # Excitation wavelength and fwhm
-    dwl = 100
+    dwl = 200
 
     # Device design parameters
     xy_size = numpy.array([10, 10])
@@ -121,105 +108,68 @@ def main():
     # #### Create the grid, mask, and draw the device ####
     grid = gridlock.Grid(edge_coords)
     epsilon = grid.allocate(n_air**2, dtype=dtype)
-    grid.draw_slab(
-        epsilon,
-        slab = dict(axis='z', center=0, span=th),
-        foreground = n_slab ** 2,
-        )
-
+    grid.draw_slab(epsilon,
+                   surface_normal=2,
+                   center=[0, 0, 0],
+                   thickness=th,
+                   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),
-        )
 
-    print(f'{grid.shape=}')
+    grid.draw_polygons(epsilon,
+                       surface_normal=2,
+                       center=[0, 0, 0],
+                       thickness=2 * th,
+                       eps=n_air**2,
+                       polygons=mask.as_polygons())
 
-    dt = dx * 0.99 / numpy.sqrt(3)
-    ee = numpy.zeros_like(epsilon, dtype=dtype)
-    hh = numpy.zeros_like(epsilon, dtype=dtype)
+    print(grid.shape)
+
+    dt = .99/numpy.sqrt(3)
+    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, thickness=pml_thickness, epsilon_eff=n_air ** 2)
+    pml_params = [[cpml_params(axis=dd, polarity=pp, dt=dt,
+                               thickness=pml_thickness, epsilon_eff=1.0**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)
+    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
+    w = 2 * numpy.pi * dx / wl
+    fwhm = dwl * w * w / (2 * numpy.pi * dx)
+    alpha = (fwhm ** 2) / 8 * numpy.log(2)
+    delay = 7/numpy.sqrt(2 * alpha)
 
-
-    # Source parameters and function. The pulse normalization is kept outside
-    # 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)
+    def field_source(i):
+        t0 = i * dt - delay
+        return numpy.sin(w * t0) * numpy.exp(-alpha * t0**2)
 
     # #### 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)
+    for t in range(max_t):
+        update_E(e, h, epsilon)
 
-        if tt < src_maxt:
-            # 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)
+        e[1][tuple(grid.shape//2)] += field_source(t)
+        update_H(e, h)
 
-        avg_rate = (tt + 1) / (time.perf_counter() - start)
+        avg_rate = (t + 1)/(time.perf_counter() - start))
+        print(f'iteration {t}: average {avg_rate} iterations per sec')
         sys.stdout.flush()
 
-        if 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=}')
+        if t % 20 == 0:
+            r = sum([(f * f * e).sum() for f, e in zip(e, epsilon)])
+            print('E sum', r)
 
         # 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
-            # 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 (t % 20 == 0 and (max_t - t <= 1000 or t <= 2000)) or t == max_t-1:
+            print('saving E-field')
+            for j, f in enumerate(e):
+                output_file['/E{}_t{}'.format('xyz'[j], t)] = f[:, :, round(f.shape[2]/2)]
 
 if __name__ == '__main__':
     main()

examples/waveguide.py

@@ -1,351 +0,0 @@
-"""
-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,12 +2,18 @@
 
 set -Eeuo pipefail
 
-ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
-cd "$ROOT"
+cd ~/projects/meanas
 
-mkdocs build --clean
+# Approach 1: pdf to html?
+#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 -
 
-PRINT_PAGE='site/print_page/index.html'
-if [[ -f "$PRINT_PAGE" ]] && command -v htmlark >/dev/null 2>&1; then
-    htmlark "$PRINT_PAGE" -o site/standalone.html
-fi
+# Approach 2: pdf to html with gladtex
+rm -rf _doc_mathimg
+pdoc --pdf --force --template-dir pdoc_templates -o doc . > doc.md
+pandoc --metadata=title:"meanas" --from=markdown+abbreviations --to=html --output=doc.htex --gladtex -s --css pdoc_templates/pdoc.css doc.md
+gladtex -a -n -d _doc_mathimg -c white -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,12 +9,9 @@ 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,
-  including mode-source and overlap-window construction.
-- `farfield`, `bloch`, `eme`: specialized helper modules for near/far transforms,
-  Bloch-periodic problems, and eigenmode expansion.
+- `waveguide_3d`: Functions for transforming `waveguide_2d` results into 3D.
 
 
 ================================================================
@@ -89,6 +86,10 @@ $$
     -\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,14 +136,6 @@ 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,
@@ -261,8 +253,8 @@ def maxwell_operator(
             h_m, h_n = b_m, b_n
         else:
             # transform from mn to xyz
-            b_xyz = (m * b_m
-                   + n * b_n)    # noqa: E128
+            b_xyz = (m * b_m[:, :, :, None]
+                   + n * b_n[:, :, :, None])
 
             # divide by mu
             temp = ifftn(b_xyz, axes=range(3))
@@ -273,7 +265,10 @@ def maxwell_operator(
             h_m = numpy.sum(h_xyz * m, axis=3)
             h_n = numpy.sum(h_xyz * n, axis=3)
 
-        return _assemble_hmn_vector(h_m, h_n)
+        h.shape = (h.size,)
+        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
 
@@ -408,8 +403,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
-                   + n * hin_n)  # noqa: E128
+            h_xyz = (m * hin_m[:, :, :, None]
+                   + n * hin_n[:, :, :, None])
 
             # multiply by mu
             temp = ifftn(h_xyz, axes=range(3))
@@ -417,8 +412,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, keepdims=True)
-            b_n = numpy.sum(b_xyz * n, axis=3, keepdims=True)
+            b_m = numpy.sum(b_xyz * m, axis=3)
+            b_n = numpy.sum(b_xyz * n, axis=3)
 
         # cross product and transform into xyz basis
         e_xyz = (n * b_m
@@ -433,7 +428,10 @@ 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
 
-        return _assemble_hmn_vector(h_m, h_n)
+        h.shape = (h.size,)
+        h = numpy.concatenate((h_m, h_n), axis=None, out=h)
+        h.shape = (h.size, 1)
+        return h
 
     return operator
 
@@ -451,7 +449,7 @@ def find_k(
         solve_callback: Callable[..., None] | None = None,
         iter_callback: Callable[..., None] | None = None,
         v0: NDArray[numpy.complex128] | None = None,
-        ) -> tuple[NDArray[numpy.float64], float, NDArray[numpy.complex128], NDArray[numpy.complex128]]:
+        ) -> tuple[float, float, NDArray[numpy.complex128], NDArray[numpy.complex128]]:
     """
     Search for a bloch vector that has a given frequency.
 
@@ -496,15 +494,15 @@ def find_k(
 
     res = scipy.optimize.minimize_scalar(
         lambda x: abs(get_f(x, band) - frequency),
-        method='bounded',
+        k_guess,
+        method='Bounded',
         bounds=k_bounds,
         options={'xatol': abs(tolerance)},
         )
 
     assert n is not None
     assert v is not None
-    actual_frequency = get_f(float(res.x), band)
-    return direction * float(res.x), float(actual_frequency), n, v
+    return float(res.x * direction), float(res.fun + frequency), n, v
 
 
 def eigsolve(
@@ -725,12 +723,7 @@ def eigsolve(
                 amax=pi,
                 )
 
-        result = scipy.optimize.minimize_scalar(
-            trace_func,
-            method='bounded',
-            bounds=(0, pi),
-            options={'xatol': tolerance},
-            )
+        result = scipy.optimize.minimize_scalar(trace_func, bounds=(0, pi), tol=tolerance)
         new_E = result.fun
         theta = result.x
 
@@ -759,7 +752,7 @@ def eigsolve(
         v = eigvecs[:, i]
         n = eigvals[i]
         v /= norm(v)
-        Av = numpy.asarray(scipy_op @ v.copy()).reshape(-1)
+        Av = (scipy_op @ v.copy())[:, 0]
         eigness = norm(Av - (v.conj() @ Av) * v)
         f = numpy.sqrt(-numpy.real(n))
         df = numpy.sqrt(-numpy.real(n) + eigness)
@@ -828,18 +821,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_norm = eR / numpy.sqrt(abs(norm2R))
-    hR_norm = hR / numpy.sqrt(abs(norm2R))
-    eL_norm = eL / numpy.sqrt(abs(norm2L))
-    hL_norm = hL / numpy.sqrt(abs(norm2L))
+    eR /= numpy.sqrt(norm2R)
+    hR /= numpy.sqrt(norm2R)
+    eL /= numpy.sqrt(norm2L)
+    hL /= numpy.sqrt(norm2L)
 
     # (eR x hL)[0] and (eL x hR)[0]
-    eRxhL_x = eR_norm[1] * hL_norm[2] - eR_norm[2] * hL_norm[1]
-    eLxhR_x = eL_norm[1] * hR_norm[2] - eL_norm[2] * hR_norm[1]
+    eRxhL_x = eR[1] * hL[2] - eR[2] - hL[1]
+    eLxhR_x = eL[1] * hR[2] - eL[2] - hR[1]
 
     #return 1j *  (eRxhL_x - eLxhR_x).sum() / numpy.sqrt(norm2R * norm2L)
     #return (eRxhL_x.sum() - eLxhR_x.sum()) / numpy.sqrt(norm2R * norm2L)
-    return eLxhR_x.sum() - eRxhL_x.sum()
+    return eRxhL_x.sum() - eLxhR_x.sum()
 
 
 def trq(
@@ -853,8 +846,8 @@ def trq(
     np = inner_product(eO, -hO, eI,  hI)
     nn = inner_product(eO, -hO, eI, -hI)
 
-    assert numpy.allclose(pp, -nn, atol=1e-12, rtol=1e-12)
-    assert numpy.allclose(pn, -np, atol=1e-12, rtol=1e-12)
+    assert pp == -nn
+    assert pn == -np
 
     logger.info(f'''
         {pp=:4g} {pn=:4g}

meanas/fdfd/eme.py

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

meanas/fdfd/farfield.py

@@ -78,12 +78,15 @@ def near_to_farfield(
     kx, ky = numpy.meshgrid(kxx, kyy, indexing='ij')
     kxy2 = kx * kx + ky * ky
     kxy = numpy.sqrt(kxy2)
-    kz = numpy.sqrt(numpy.maximum(0, k * k - kxy2))
+    kz = numpy.sqrt(k * k - kxy2)
 
-    sin_th = numpy.divide(ky, kxy, out=numpy.zeros_like(ky), where=kxy != 0)
-    cos_th = numpy.divide(kx, kxy, out=numpy.ones_like(kx), where=kxy != 0)
+    sin_th = ky / kxy
+    cos_th = kx / kxy
     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
@@ -111,8 +114,8 @@ def near_to_farfield(
     outputs = {
         'E': E_far,
         'H': H_far,
-        'dkx': float(kxx[1] - kxx[0]),
-        'dky': float(kyy[1] - kyy[0]),
+        'dkx': kx[1] - kx[0],
+        'dky': ky[1] - ky[0],
         'kx': kx,
         'ky': ky,
         'theta': theta,
@@ -174,19 +177,22 @@ def far_to_nearfield(
     padded_shape = cast('Sequence[int]', padded_size)
 
     k = 2 * pi
-    kxs = dkx * fftshift(fftfreq(s[0], d=1 / s[0]))
-    kys = dky * fftshift(fftfreq(s[1], d=1 / s[1]))
+    kxs = fftshift(fftfreq(s[0], 1 / (s[0] * dkx)))
+    kys = fftshift(fftfreq(s[0], 1 / (s[1] * dky)))
 
     kx, ky = numpy.meshgrid(kxs, kys, indexing='ij')
     kxy2 = kx * kx + ky * ky
     kxy = numpy.sqrt(kxy2)
 
-    kz = numpy.sqrt(numpy.maximum(0, k * k - kxy2))
+    kz = numpy.sqrt(k * k - kxy2)
 
-    sin_th = numpy.divide(ky, kxy, out=numpy.zeros_like(ky), where=kxy != 0)
-    cos_th = numpy.divide(kx, kxy, out=numpy.ones_like(kx), where=kxy != 0)
+    sin_th = ky / kxy
+    cos_th = kx / kxy
     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
@@ -206,41 +212,21 @@ def far_to_nearfield(
     N = [L[1],
         -L[0]]  # noqa: E128
 
-    En_fft = [
-        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,
-            ),
-        ]
+    En_fft = [-( L[0] * sin_th + L[1] * cos_phi * cos_th) / cos_phi,
+              -(-L[0] * cos_th + L[1] * cos_phi * sin_th) / cos_phi]
 
-    Hn_fft = [
-        numpy.divide(
-            N[0] * sin_th + N[1] * cos_phi * cos_th,
-            cos_phi,
-            out=numpy.zeros_like(N[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,
-            ),
-        ]
+    Hn_fft = [( N[0] * sin_th + N[1] * cos_phi * cos_th) / cos_phi,
+              (-N[0] * cos_th + N[1] * cos_phi * sin_th) / cos_phi]
+
+    for i in range(2):
+        En_fft[i][cos_phi == 0] = 0
+        Hn_fft[i][cos_phi == 0] = 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[1] * dky)
+    dy = 2 * pi / (s[0] * dky)
 
     outputs = {
         'E': E_near,
@@ -250,3 +236,4 @@ 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_t) -> cfdfield_t:
-        curls = ch(ce(e) / mu)          # type: ignore   # mu = None ok because we don't return the function
+    def op_mu(e: cfdfield) -> cfdfield_t:
+        curls = ch(mu * ce(e))          # 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_t) -> cfdfield_t:
-        J = ch(m / mu) / (1j * omega)           # type: ignore  # mu=None ok
+    def m2j_mu(m: cfdfield) -> cfdfield_t:
+        J = ch(m / mu) / (-1j * omega)          # type: ignore  # mu=None ok
         return cfdfield_t(J)
 
-    def m2j_1(m: cfdfield_t) -> cfdfield_t:
-        J = ch(m) / (1j * omega)
+    def m2j_1(m: cfdfield) -> cfdfield_t:
+        J = ch(m) / (-1j * omega)
         return cfdfield_t(J)
 
     if mu is None:
@@ -158,23 +158,10 @@ 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).
@@ -188,6 +175,7 @@ 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:
@@ -200,13 +188,7 @@ 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$.
-
-    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(...)`.
+    for the Poynting vector, $S = E \times H$
 
     Note:
         This function also shifts the input `E` field by one cell as required
@@ -222,8 +204,7 @@ 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 the staggered-grid cross product `E \times H`.
-        For time-average power, call it as `f(E, H.conj())` and take `Re(...) / 2`.
+        Function `f` that returns E x H as required for the poynting vector.
     """
     def exh(e: cfdfield, h: cfdfield) -> cfdfield_t:
         s = numpy.empty_like(e)

meanas/fdfd/operators.py

@@ -236,12 +236,10 @@ def eh_full(
     else:
         pm = sparse.diags_array(numpy.where(pmc, 0, 1))    # set pm to (not PMC)
 
-    iwe = pe @ (1j * omega * sparse.diags(epsilon)) @ pe
-    if mu is None:
-        iwm = 1j * omega * sparse.eye(epsilon.size)
-    else:
-        iwm = 1j * omega * sparse.diags(mu)
-
+    iwe = pe @ (1j * omega * sparse.diags_array(epsilon)) @ pe
+    iwm = 1j * omega
+    if mu is not None:
+        iwm *= sparse.diags_array(mu)
     iwm = pm @ iwm @ pm
 
     A1 = pe @ curl_back(dxes[1]) @ pm
@@ -310,22 +308,16 @@ def m2j(
 
 
 def poynting_e_cross(e: vcfdfield, dxes: dx_lists_t) -> sparse.sparray:
-    r"""
-    Operator for computing the staggered-grid `(E \times)` part of the Poynting vector.
-
-    On the Yee grid the E and H components live on different edges, so the
-    electric field must be shifted by one cell in the appropriate direction
-    before forming the discrete cross product. This sparse operator encodes that
-    shifted cross product directly and is the matrix equivalent of
-    `functional.poynting_e_cross_h(...)`.
+    """
+    Operator for computing the Poynting vector, containing the
+    (E x) portion of the Poynting vector.
 
     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 the `(E \times)` part of the staggered Poynting
-        cross product.
+        Sparse matrix containing (E x) portion of Poynting cross product.
     """
     shape = [len(dx) for dx in dxes[0]]
 
@@ -345,26 +337,15 @@ 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 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.
+    """
+    Operator for computing the Poynting vector, containing the (H x) portion of the Poynting vector.
 
     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 the `(H \times)` part of the staggered Poynting
-        cross product.
+        Sparse matrix containing (H x) portion of Poynting cross product.
     """
     shape = [len(dx) for dx in dxes[0]]
 
@@ -389,23 +370,11 @@ 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.
 
-    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.
+    TODO: Reference Rumpf paper
 
     Args:
         TF_region: Mask, which is set to 1 inside the total-field region and 0 in the
@@ -417,7 +386,9 @@ 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)
@@ -431,17 +402,11 @@ 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,11 +128,6 @@ 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]
@@ -158,6 +153,9 @@ 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(None)
+        else:
             slc = slice(-thickness, None)
 
     dx_ai[slc] *= 1 + 1j * s_function(l_d(pos_a[slc])) / d / s_correction

meanas/fdfd/solvers.py

@@ -48,11 +48,9 @@ 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)
-            callback(xk)
+            kwargs['callback'](xk)
 
         kwargs['callback'] = augmented_callback
     else:
@@ -120,15 +118,15 @@ def generic(
     Pl, Pr = operators.e_full_preconditioners(dxes)
 
     if adjoint:
-        A = (Pl @ A0 @ Pr).T.conjugate()
-        b = Pr.T.conjugate() @ b0
+        A = (Pl @ A0 @ Pr).H
+        b = Pr.H @ b0
     else:
         A = Pl @ A0 @ Pr
         b = Pl @ b0
 
     if E_guess is not None:
         if adjoint:
-            x0 = Pr.T.conjugate() @ E_guess
+            x0 = Pr.H @ E_guess
         else:
             x0 = Pl @ E_guess
         matrix_solver_opts['x0'] = x0
@@ -136,7 +134,7 @@ def generic(
     x = matrix_solver(A.tocsr(), b, **matrix_solver_opts)
 
     if adjoint:
-        x0 = Pl.T.conjugate() @ x
+        x0 = Pl.H @ x
     else:
         x0 = Pr @ x
 

meanas/fdfd/waveguide_2d.py

@@ -175,6 +175,8 @@ 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
@@ -337,7 +339,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.
 
@@ -355,21 +357,6 @@ 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
@@ -387,7 +374,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.
 
@@ -405,13 +392,6 @@ 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
@@ -429,25 +409,7 @@ def _normalized_fields(
         mu: vfdslice | None = None,
         prop_phase: float = 0,
         ) -> tuple[vcfdslice_t, vcfdslice_t]:
-    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.
-    """
+    # TODO documentation
     shape = [s.size for s in dxes[0]]
 
     # Find time-averaged Sz and normalize to it
@@ -959,7 +921,7 @@ def solve_mode(
     return vcfdfield2_t(e_xys[0]), wavenumbers[0]
 
 
-def inner_product(
+def inner_product(    # TODO documentation
         e1: vcfdfield2,
         h2: vcfdfield2,
         dxes: dx_lists2_t,
@@ -967,36 +929,6 @@ def inner_product(
         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]]
 
@@ -1019,3 +951,5 @@ def inner_product(
         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,25 +3,8 @@ 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
@@ -41,7 +24,7 @@ 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.
 
@@ -52,22 +35,19 @@ 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]` must select exactly one item.
+                as the waveguide cross-section. `slices[axis]` should select only one item.
         epsilon: Dielectric constant
         mu: Magnetic permeability (default 1 everywhere)
 
     Returns:
-        Dictionary containing:
-
-        - `E`: full-grid electric field for the solved mode
-        - `H`: full-grid magnetic field for the solved mode
-        - `wavenumber`: propagation constant corrected for the discretized
-          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.
+        ```
+        {
+            'E': NDArray[complexfloating],
+            'H': NDArray[complexfloating],
+            'wavenumber': complex,
+            'wavenumber_2d': complex,
+        }
+        ```
     """
     if mu is None:
         mu = numpy.ones_like(epsilon)
@@ -157,14 +137,7 @@ def compute_source(
         mu: Magnetic permeability (default 1 everywhere)
 
     Returns:
-        `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`
+        J distribution for the unidirectional source
     """
     E_expanded = expand_e(E=E, dxes=dxes, wavenumber=wavenumber, axis=axis,
                           polarity=polarity, slices=slices)
@@ -184,52 +157,19 @@ def compute_source(
 
 
 def compute_overlap_e(
-        E: cfdfield_t,
+        E: cfdfield,
         wavenumber: complex,
         dxes: dx_lists_t,
         axis: int,
         polarity: int,
         slices: Sequence[slice],
-        omega: float,
         ) -> cfdfield_t:
-    r"""
-    Build an overlap field for projecting another 3D electric field onto a mode.
-
-    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)
+    """
+    Given an eigenmode obtained by `solve_mode`, calculates an overlap_e for the
+    mode orthogonality relation Integrate(((E x H_mode) + (E_mode x H)) dot dn)
+    [assumes reflection symmetry].
 
+    TODO: add reference or derivation for compute_overlap_e
 
     Args:
         E: E-field of the mode
@@ -241,41 +181,25 @@ def compute_overlap_e(
                 as the waveguide cross-section. slices[axis] should select only one item.
 
     Returns:
-        `overlap_e` normalized so that `numpy.sum(overlap_e * E.conj()) == 1`
-        over the retained overlap window.
+        overlap_e such that `numpy.sum(overlap_e * other_e.conj())` computes the overlap integral
     """
     slices = tuple(slices)
 
     Ee = expand_e(E=E, wavenumber=wavenumber, dxes=dxes,
                   axis=axis, polarity=polarity, slices=slices)
 
-    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)
+    start, stop = sorted((slices[axis].start, slices[axis].start - 2 * polarity))
 
     slices2_l = list(slices)
-    slices2_l[axis] = slice(clipped_start, clipped_stop)
+    slices2_l[axis] = slice(start, stop)
     slices2 = (slice(None), *slices2_l)
 
     Etgt = numpy.zeros_like(Ee)
     Etgt[slices2] = Ee[slices2]
 
-    # Note: We normalize so that (Etgt @ E.conj()) == 1, so (Etgt @ Etgt.conj) != 1
-    norm = (Etgt.conj() * Etgt).sum()
-    if norm == 0:
-        raise ValueError('Requested overlap window contains no overlap field support')
-    Etgt /= norm
+    # note no sqrt() when normalizing below since we want to get 1.0 after overlapping with the
+    # original field, not the normalized one
+    Etgt /= (Etgt.conj() * Etgt).sum()    # type: ignore
     return cfdfield_t(Etgt)
 
 
@@ -287,7 +211,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
@@ -305,16 +229,6 @@ 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,125 +2,63 @@ r"""
 Operators and helper functions for cylindrical waveguides with unchanging cross-section.
 
 Waveguide operator is derived according to 10.1364/OL.33.001848.
-
-As in `waveguide_2d`, the propagation dependence is separated from the
-transverse solve. Here the propagation coordinate is the bend angle `\theta`,
-and the fields are assumed to have the form
-
-$$
-\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
+The curl equations in the complex coordinate system become
 
 $$
 \begin{aligned}
-r_a(n) &= r_{\min} + \sum_{j \le n} \Delta r_{e, j}, \\
-r_b\!\left(n + \tfrac{1}{2}\right) &= r_{\min} + \tfrac{1}{2}\Delta r_{e, n}
-  + \sum_{j < n} \Delta r_{h, j},
+-\imath \omega \mu_{xx} H_x &= \tilde{\partial}_y E_z + \imath \beta frac{E_y}{\tilde{t}_x} \\
+-\imath \omega \mu_{yy} H_y &= -\imath \beta E_x - \frac{1}{\hat{t}_x} \tilde{\partial}_x \tilde{t}_x E_z \\
+-\imath \omega \mu_{zz} H_z &= \tilde{\partial}_x E_y - \tilde{\partial}_y E_x \\
+\imath \omega \epsilon_{xx} E_x &= \hat{\partial}_y H_z + \imath \beta \frac{H_y}{\hat{T}} \\
+\imath \omega \epsilon_{yy} E_y &= -\imath \beta H_x - \{1}{\tilde{t}_x} \hat{\partial}_x \hat{t}_x} H_z \\
+\imath \omega \epsilon_{zz} E_z &= \hat{\partial}_x H_y - \hat{\partial}_y H_x \\
 \end{aligned}
 $$
 
-and from them the diagonal metric matrices
+where $t_x = 1 + \frac{\Delta_{x, m}}{R_0}$ is the grid spacing adjusted by the nominal radius $R0$.
+
+Rewrite the last three equations as
 
 $$
 \begin{aligned}
-T_a &= \operatorname{diag}(r_a / r_{\min}), \\
-T_b &= \operatorname{diag}(r_b / r_{\min}).
+\imath \beta H_y &=  \imath \omega \hat{t}_x \epsilon_{xx} E_x - \hat{t}_x \hat{\partial}_y H_z \\
+\imath \beta H_x &= -\imath \omega \hat{t}_x \epsilon_{yy} E_y - \hat{t}_x \hat{\partial}_x H_z \\
+\imath \omega E_z &= \frac{1}{\epsilon_{zz}} \hat{\partial}_x H_y - \frac{1}{\epsilon_{zz}} \hat{\partial}_y H_x \\
 \end{aligned}
 $$
 
-With the same forward/backward derivative notation used in `waveguide_2d`, the
-coordinate-transformed discrete curl equations used here are
+The derivation then follows the same steps as the straight waveguide, leading to the eigenvalue problem
+
+$$
+\beta^2 \begin{bmatrix} E_x \\
+                        E_y \end{bmatrix} =
+    (\omega^2 \begin{bmatrix} T_b T_b \mu_{yy} \epsilon_{xx} & 0 \\
+                                                            0 & T_a T_a \mu_{xx} \epsilon_{yy} \end{bmatrix} +
+              \begin{bmatrix} -T_b \mu_{yy} \hat{\partial}_y \\
+                               T_a \mu_{xx} \hat{\partial}_x \end{bmatrix} T_b \mu_{zz}^{-1}
+              \begin{bmatrix} -\tilde{\partial}_y & \tilde{\partial}_x \end{bmatrix} +
+      \begin{bmatrix} \tilde{\partial}_x \\
+                      \tilde{\partial}_y \end{bmatrix} T_a \epsilon_{zz}^{-1}
+                 \begin{bmatrix} \hat{\partial}_x T_b \epsilon_{xx} & \hat{\partial}_y T_a \epsilon_{yy} \end{bmatrix})
+    \begin{bmatrix} E_x \\
+                    E_y \end{bmatrix}
+$$
+
+which resembles the straight waveguide eigenproblem with additonal $T_a$ and $T_b$ terms. These
+are diagonal matrices containing the $t_x$ values:
 
 $$
 \begin{aligned}
--\imath \omega \mu_{rr} H_r &= \tilde{\partial}_y E_z + \imath \beta T_a^{-1} E_y, \\
--\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).
+T_a &=  1 + \frac{\Delta_{x, m               }}{R_0}
+T_b &=  1 + \frac{\Delta_{x, m + \frac{1}{2} }}{R_0}
 \end{aligned}
+
+
+TODO: consider 10.1364/OE.20.021583 for an alternate approach
 $$
 
-The first three equations are the cylindrical analogue of the straight-guide
-relations for `H_r`, `H_y`, and `H_z`. The next two are the metric-weighted
-versions of the straight-guide identities for `\imath \beta H_y` and
-`\imath \beta H_r`, and the last equation plays the same role as the
-longitudinal `E_z` reconstruction in `waveguide_2d`.
-
-Following the same elimination steps as in `waveguide_2d`, apply
-`\imath \beta \tilde{\partial}_r` and `\imath \beta \tilde{\partial}_y` to the
-equation for `E_z`, substitute for `\imath \beta H_r` and `\imath \beta H_y`,
-and then eliminate `H_z` with
-
-$$
-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, ...]]]`.
+As in the straight waveguide case, all the functions in this file assume a 2D grid
+ (i.e. `dxes = [[[dr_e_0, dx_e_1, ...], [dy_e_0, ...]], [[dr_h_0, ...], [dy_h_0, ...]]]`).
 """
 from typing import Any, cast
 from collections.abc import Sequence
@@ -156,18 +94,17 @@ 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_r \\
+        \begin{bmatrix} E_x \\
                         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] = beta**2 * [E_r, E_y]
+        A @ [E_r, E_y] = wavenumber**2 * [E_r, E_y]
 
     which can then be solved for the eigenmodes of the system
-    (an `exp(-i * angular_wavenumber * theta)` theta-dependence is assumed for
-    the fields, with `beta = angular_wavenumber / rmin`).
+    (an `exp(-i * wavenumber * theta)` theta-dependence is assumed for the fields).
 
     (NOTE: See module docs and 10.1364/OL.33.001848)
 
@@ -333,7 +270,7 @@ def exy2h(
         mu: vfdslice | None = None
         ) -> sparse.sparray:
     """
-    Operator which transforms the vector `e_xy` containing the vectorized E_r and E_y fields,
+    Operator which transforms the vector `e_xy` containing the vectorized E_x and E_y fields,
      into a vectorized H containing all three H components
 
     Args:
@@ -361,11 +298,11 @@ def exy2e(
         epsilon: vfdslice,
         ) -> sparse.sparray:
     """
-    Operator which transforms the vector `e_xy` containing the vectorized E_r and E_y fields,
+    Operator which transforms the vector `e_xy` containing the vectorized E_x and E_y fields,
      into a vectorized E containing all three E components
 
     Unlike the straight waveguide case, the H_z components do not cancel and must be calculated
-    from E_r and E_y in order to then calculate E_z.
+    from E_x and E_y in order to then calculate E_z.
 
     Args:
         angular_wavenumber: Wavenumber assuming fields have theta-dependence of
@@ -423,10 +360,9 @@ def e2h(
     This operator is created directly from the initial coordinate-transformed equations:
     $$
     \begin{aligned}
-    -\imath \omega \mu_{rr} H_r &= \tilde{\partial}_y E_z + \imath \beta T_a^{-1} E_y, \\
-    -\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 \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}
     $$
 
@@ -461,18 +397,15 @@ def dxes2T(
         rmin: float,
         ) -> tuple[NDArray[numpy.float64], NDArray[numpy.float64]]:
     r"""
-    Construct the cylindrical metric matrices $T_a$ and $T_b$.
-
-    `T_a` is sampled on the E-grid radial locations, while `T_b` is sampled on
-    the staggered H-grid radial locations. These are the diagonal matrices that
-    convert the straight-waveguide algebra into its cylindrical counterpart.
+    Returns the $T_a$ and $T_b$ diagonal matrices which are used to apply the cylindrical
+      coordinate transformation in various operators.
 
     Args:
         dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
         rmin: Radius at the left edge of the simulation domain (at minimum 'x')
 
     Returns:
-        Sparse diagonal matrices `(T_a, T_b)`.
+        Sparse matrix representations of the operators Ta and Tb
     """
     ra = rmin + numpy.cumsum(dxes[0][0])                      # Radius at Ey points
     rb = rmin + dxes[0][0] / 2.0 + numpy.cumsum(dxes[1][0])   # Radius at Ex points
@@ -494,12 +427,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_r and E_y fields,
+    """
+    Given a vector `e_xy` containing the vectorized E_x and E_y fields,
      returns normalized, vectorized E and H fields for the system.
 
     Args:
-        e_xy: Vector containing E_r and E_y fields
+        e_xy: Vector containing E_x and E_y fields
         angular_wavenumber: Wavenumber assuming fields have theta-dependence of
             `exp(-i * angular_wavenumber * theta)`. It should satisfy
             `operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy`
@@ -514,18 +447,6 @@ 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
@@ -544,30 +465,8 @@ 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  # pragma: no cover
+    pass
 
 @overload
 def vec(f: fdfield_t) -> vfdfield_t:
-    pass  # pragma: no cover
+    pass
 
 @overload
 def vec(f: cfdfield_t) -> vcfdfield_t:
-    pass  # pragma: no cover
+    pass
 
 @overload
 def vec(f: fdfield2_t) -> vfdfield2_t:
-    pass  # pragma: no cover
+    pass
 
 @overload
 def vec(f: cfdfield2_t) -> vcfdfield2_t:
-    pass  # pragma: no cover
+    pass
 
 @overload
 def vec(f: fdslice_t) -> vfdslice_t:
-    pass  # pragma: no cover
+    pass
 
 @overload
 def vec(f: cfdslice_t) -> vcfdslice_t:
-    pass  # pragma: no cover
+    pass
 
 @overload
 def vec(f: ArrayLike) -> NDArray:
-    pass  # pragma: no cover
+    pass
 
 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  # pragma: no cover
+    pass
 
 @overload
 def unvec(v: vfdfield_t, shape: Sequence[int], nvdim: int = 3) -> fdfield_t:
-    pass  # pragma: no cover
+    pass
 
 @overload
 def unvec(v: vcfdfield_t, shape: Sequence[int], nvdim: int = 3) -> cfdfield_t:
-    pass  # pragma: no cover
+    pass
 
 @overload
 def unvec(v: vfdfield2_t, shape: Sequence[int], nvdim: int = 3) -> fdfield2_t:

meanas/fdtd/__init__.py

@@ -144,18 +144,6 @@ 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
 
@@ -168,44 +156,7 @@ t=0 (assuming the source is off for t<0 this gives $\sim 10^{-3}$ error at t=0).
 
 Boundary conditions
 ===================
-
-`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.
-$$
+# TODO notes about boundaries / PMLs
 """
 
 from .base import (
@@ -227,9 +178,3 @@ 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,19 +28,17 @@ 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] = 0
-            e[u][boundary] = e[u][shifted1]
-            e[v][boundary] = e[v][shifted1]
+            e[direction][boundary_slice] = 0
+            e[u][boundary_slice] = e[u][shifted1_slice]
+            e[v][boundary_slice] = e[v][shifted1_slice]
             return e
 
         def hn(h: fdfield_t) -> fdfield_t:
-            h[direction][boundary] = h[direction][shifted1]
-            h[u][boundary] = 0
-            h[v][boundary] = 0
+            h[direction][boundary_slice] = h[direction][shifted1_slice]
+            h[u][boundary_slice] = 0
+            h[v][boundary_slice] = 0
             return h
 
         return en, hn
@@ -52,23 +50,20 @@ 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] = -e[direction][shifted2]
-            e[direction][shifted1] = 0
-            e[u][boundary] = e[u][shifted1]
-            e[v][boundary] = e[v][shifted1]
+            e[direction][boundary_slice] = -e[direction][shifted2_slice]
+            e[direction][shifted1_slice] = 0
+            e[u][boundary_slice] = e[u][shifted1_slice]
+            e[v][boundary_slice] = e[v][shifted1_slice]
             return e
 
         def hp(h: fdfield_t) -> fdfield_t:
-            h[direction][boundary] = h[direction][shifted1]
-            h[u][boundary] = -h[u][shifted2]
-            h[u][shifted1] = 0
-            h[v][boundary] = -h[v][shifted2]
-            h[v][shifted1] = 0
+            h[direction][boundary_slice] = h[direction][shifted1_slice]
+            h[u][boundary_slice] = -h[u][shifted2_slice]
+            h[u][shifted1_slice] = 0
+            h[v][boundary_slice] = -h[v][shifted2_slice]
+            h[v][shifted1_slice] = 0
             return h
 
         return ep, hp

meanas/fdtd/energy.py

@@ -4,21 +4,7 @@ from ..fdmath import dx_lists_t, fdfield_t, fdfield
 from ..fdmath.functional import deriv_back
 
 
-"""
-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.
-"""
+# TODO documentation
 
 
 def poynting(
@@ -266,23 +252,13 @@ def delta_energy_j(
         e1: fdfield,
         dxes: dx_lists_t | None = None,
         ) -> fdfield_t:
-    r"""
-    Calculate the electric-current work term $J \cdot E$ on the Yee grid.
+    """
+    Calculate
 
-    This is the source contribution that appears beside the flux divergence in
-    the discrete Poynting identities documented in `meanas.fdtd`.
+    Note that each value of $J$ contributes to the energy twice (i.e. once per field update)
+    despite only causing the value of $E$ to change once (same for $M$ and $H$).
 
-    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))
@@ -301,20 +277,6 @@ 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:
-            envelope = numpy.where(t0 > 0, 1.0, envelope)
+        if one_sided and t0 > 0:
+            envelope = 1
 
         cc = numpy.cos(omega * t0)
         ss = numpy.sin(omega * t0)
@@ -94,14 +94,10 @@ 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 * tt * tt / 4) - turn_on,
-        x0=0,
-        x1=-2 / omega,
-        )
-
+    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 = delay_results.root
     delay = numpy.ceil(delay * freq) / freq     # force delay to integer number of periods to maintain phase
 
@@ -160,11 +156,11 @@ def gaussian_beam(
     wz = numpy.sqrt(wz2)        # == fwhm(z) / sqrt(2 * ln(2))
 
     kk = 2 * pi / wl
-    inv_Rz = numpy.divide(zz, z2 + zr2, out=numpy.zeros_like(zz), where=(z2 + zr2) != 0)
+    Rz = zz * (1 + zr2 / z2)
     gouy = numpy.arctan(zz / zr)
 
-    gaussian = w0 / wz * numpy.exp(-r2 / wz2) * numpy.exp(1j * (kk * zz + kk * r2 * inv_Rz / 2 - gouy))
+    gaussian = w0 / wz * numpy.exp(-r2 / wz2) * numpy.exp(1j * (kk * zz + kk * r2 / 2 / Rz - gouy))
 
     row = gaussian[:, :, gaussian.shape[2] // 2]
-    norm = numpy.linalg.norm(row)
+    norm = numpy.sqrt((row * row.conj()).sum())
     return gaussian / norm

meanas/fdtd/phasor.py

@@ -1,126 +0,0 @@
-"""
-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,19 +1,9 @@
 """
-Convolutional perfectly matched layer (CPML) support for FDTD updates.
+PML implementations
 
-The helpers in this module construct per-face CPML parameters and then wrap the
-standard Yee updates with the additional auxiliary `psi` fields needed by the
-CPML recurrence.
+#TODO discussion of PMLs
+#TODO cpml documentation
 
-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!
 
@@ -42,29 +32,6 @@ 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}')
@@ -90,6 +57,8 @@ 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)
@@ -113,6 +82,8 @@ 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 {
@@ -131,27 +102,6 @@ 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

@@ -1,37 +0,0 @@
-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

@@ -1,49 +0,0 @@
-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

@@ -1,56 +0,0 @@
-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

@@ -1,22 +0,0 @@
-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 make_prng
+from .utils import PRNG
 
 
 FixtureRequest = Any
@@ -42,7 +42,6 @@ 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':
@@ -58,11 +57,9 @@ 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(low=min(epsilon_bg, epsilon_fg),
                                   high=max(epsilon_bg, epsilon_fg),
-            size=shape,
-        )
+                                  size=shape)
 
     return epsilon
 
@@ -83,7 +80,6 @@ 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':
@@ -92,7 +88,8 @@ 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

@@ -1,73 +0,0 @@
-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

@@ -1,315 +0,0 @@
-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

@@ -1,101 +0,0 @@
-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

@@ -1,132 +0,0 @@
-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

@@ -1,213 +0,0 @@
-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

@@ -1,100 +0,0 @@
-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

@@ -1,122 +0,0 @@
-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

@@ -1,126 +0,0 @@
-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

@@ -1,60 +0,0 @@
-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

@@ -1,90 +0,0 @@
-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

@@ -1,46 +0,0 @@
-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, make_prng
+from .utils import assert_close, assert_fields_close, PRNG
 from .conftest import FixtureRequest
 
 
@@ -179,14 +179,13 @@ 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

@@ -1,45 +0,0 @@
-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

@@ -1,62 +0,0 @@
-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

@@ -1,99 +0,0 @@
-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

@@ -1,42 +0,0 @@
-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

@@ -1,208 +0,0 @@
-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

@@ -1,44 +0,0 @@
-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

@@ -1,51 +0,0 @@
-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

@@ -1,284 +0,0 @@
-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

@@ -1,436 +0,0 @@
-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

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

mkdocs.yml

@@ -1,76 +0,0 @@
-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

@@ -0,0 +1,47 @@
+<%!
+    # 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

@@ -0,0 +1,389 @@
+<%!
+    from pdoc.html_helpers import minify_css
+%>
+
+<%def name="mobile()" filter="minify_css">
+  .flex {
+    display: flex !important;
+  }
+
+  body {
+    line-height: 1.5em;
+    background: black;
+    color: #DDD;
+  }
+
+  #content {
+    padding: 20px;
+  }
+
+  #sidebar {
+    padding: 30px;
+    overflow: hidden;
+  }
+
+  .http-server-breadcrumbs {
+    font-size: 130%;
+    margin: 0 0 15px 0;
+  }
+
+  #footer {
+    font-size: .75em;
+    padding: 5px 30px;
+    border-top: 1px solid #ddd;
+    text-align: right;
+  }
+    #footer p {
+      margin: 0 0 0 1em;
+      display: inline-block;
+    }
+    #footer p:last-child {
+      margin-right: 30px;
+    }
+
+  h1, h2, h3, h4, h5 {
+    font-weight: 300;
+  }
+  h1 {
+    font-size: 2.5em;
+    line-height: 1.1em;
+  }
+  h2 {
+    font-size: 1.75em;
+    margin: 1em 0 .50em 0;
+  }
+  h3 {
+    font-size: 1.4em;
+    margin: 25px 0 10px 0;
+  }
+  h4 {
+    margin: 0;
+    font-size: 105%;
+  }
+
+  a {
+    color: #999;
+    text-decoration: none;
+    transition: color .3s ease-in-out;
+  }
+  a:hover {
+    color: #18d;
+  }
+
+  .title code {
+    font-weight: bold;
+  }
+  h2[id^="header-"] {
+    margin-top: 2em;
+  }
+  .ident {
+    color: #7ff;
+  }
+
+  pre code {
+    background: transparent;
+    font-size: .8em;
+    line-height: 1.4em;
+  }
+  code {
+    background: #0d0d0e;
+    padding: 1px 4px;
+    overflow-wrap: break-word;
+  }
+  h1 code { background: transparent }
+
+  pre {
+    background: #111;
+    border: 0;
+    border-top: 1px solid #ccc;
+    border-bottom: 1px solid #ccc;
+    margin: 1em 0;
+    padding: 1ex;
+  }
+
+  #http-server-module-list {
+    display: flex;
+    flex-flow: column;
+  }
+    #http-server-module-list div {
+      display: flex;
+    }
+    #http-server-module-list dt {
+      min-width: 10%;
+    }
+    #http-server-module-list p {
+      margin-top: 0;
+    }
+
+  .toc ul,
+  #index {
+    list-style-type: none;
+    margin: 0;
+    padding: 0;
+  }
+    #index code {
+      background: transparent;
+    }
+    #index h3 {
+      border-bottom: 1px solid #ddd;
+    }
+    #index ul {
+      padding: 0;
+    }
+    #index h4 {
+      font-weight: bold;
+    }
+    #index h4 + ul {
+      margin-bottom:.6em;
+    }
+    /* Make TOC lists have 2+ columns when viewport is wide enough.
+       Assuming ~20-character identifiers and ~30% wide sidebar. */
+    @media (min-width: 200ex) { #index .two-column { column-count: 2 } }
+    @media (min-width: 300ex) { #index .two-column { column-count: 3 } }
+
+  dl {
+    margin-bottom: 2em;
+  }
+    dl dl:last-child {
+      margin-bottom: 4em;
+    }
+  dd {
+    margin: 0 0 1em 3em;
+  }
+    #header-classes + dl > dd {
+      margin-bottom: 3em;
+    }
+    dd dd {
+      margin-left: 2em;
+    }
+    dd p {
+      margin: 10px 0;
+    }
+    .name {
+      background: #111;
+      font-weight: bold;
+      font-size: .85em;
+      padding: 5px 10px;
+      display: inline-block;
+      min-width: 40%;
+    }
+      .name:hover {
+        background: #101010;
+      }
+      .name > span:first-child {
+        white-space: nowrap;
+      }
+      .name.class > span:nth-child(2) {
+        margin-left: .4em;
+      }
+    .inherited {
+      color: #777;
+      border-left: 5px solid #eee;
+      padding-left: 1em;
+    }
+    .inheritance em {
+      font-style: normal;
+      font-weight: bold;
+    }
+
+    /* Docstrings titles, e.g. in numpydoc format */
+    .desc h2 {
+      font-weight: 400;
+      font-size: 1.25em;
+    }
+    .desc h3 {
+      font-size: 1em;
+    }
+    .desc dt code {
+      background: inherit;  /* Don't grey-back parameters */
+    }
+
+    .source summary,
+    .git-link-div {
+      color: #aaa;
+      text-align: right;
+      font-weight: 400;
+      font-size: .8em;
+      text-transform: uppercase;
+    }
+      .source summary > * {
+        white-space: nowrap;
+        cursor: pointer;
+      }
+      .git-link {
+        color: inherit;
+        margin-left: 1em;
+      }
+    .source pre {
+      max-height: 500px;
+      overflow: auto;
+      margin: 0;
+    }
+    .source pre code {
+      font-size: 12px;
+      overflow: visible;
+    }
+  .hlist {
+    list-style: none;
+  }
+    .hlist li {
+      display: inline;
+    }
+    .hlist li:after {
+      content: ',\2002';
+    }
+    .hlist li:last-child:after {
+      content: none;
+    }
+    .hlist .hlist {
+      display: inline;
+      padding-left: 1em;
+    }
+
+  img {
+    max-width: 100%;
+  }
+
+  .admonition {
+    padding: .1em .5em;
+    margin-bottom: 1em;
+  }
+    .admonition-title {
+      font-weight: bold;
+    }
+    .admonition.note,
+    .admonition.info,
+    .admonition.important {
+      background: #610;
+    }
+    .admonition.todo,
+    .admonition.versionadded,
+    .admonition.tip,
+    .admonition.hint {
+      background: #202;
+    }
+    .admonition.warning,
+    .admonition.versionchanged,
+    .admonition.deprecated {
+      background: #02b;
+    }
+    .admonition.error,
+    .admonition.danger,
+    .admonition.caution {
+      background: darkpink;
+    }
+</%def>
+
+<%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

@@ -0,0 +1,445 @@
+<%
+  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>