[docs] expand API and derivation docs

README.md

@@ -94,6 +94,59 @@ python3 -m pytest -rsxX | tee test_results.txt
 
 ## Use
 
-See `examples/` for some simple examples; you may need additional
+`meanas` is organized around a few core workflows:
-packages such as [gridlock](https://mpxd.net/code/jan/gridlock)
+
-to run the examples.
+- `meanas.fdfd`: frequency-domain wave equations, sparse operators, SCPML, and
+  iterative solves for driven problems.
+- `meanas.fdfd.waveguide_2d` / `meanas.fdfd.waveguide_3d`: waveguide mode
+  solvers, mode-source construction, and overlap windows for port-based
+  excitation and analysis.
+- `meanas.fdtd`: Yee-step updates, CPML boundaries, flux/energy accounting, and
+  on-the-fly phasor extraction for comparing time-domain runs against FDFD.
+- `meanas.fdmath`: low-level finite-difference operators, vectorization helpers,
+  and derivations shared by the FDTD and FDFD layers.
+
+The most mature user-facing workflows are:
+
+1. Build an FDFD operator or waveguide port source, then solve a driven
+   frequency-domain problem.
+2. Run an FDTD simulation, extract one or more frequency-domain phasors with
+   `meanas.fdtd.accumulate_phasor(...)`, and compare those phasors against an
+   FDFD reference on the same Yee grid.
+
+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.

examples/fdtd.py

@@ -1,7 +1,12 @@
 """
-Example code for running an FDTD simulation
+Example code for a broadband FDTD run with phasor extraction.
 
-See main() for simulation setup.
+This script shows the intended low-level workflow for:
+
+1. building a Yee-grid simulation with CPML on all faces,
+2. driving it with an electric-current pulse,
+3. extracting a single-frequency phasor on the fly, and
+4. checking that phasor against the matching stretched-grid FDFD operator.
 """
 
 import sys
@@ -150,7 +155,8 @@ def main():
     # print(f'Save time interval would be {sample_interval} * dt = {sample_interval * dt:3g}')
 
 
-    # Source parameters and function
+    # 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
@@ -170,7 +176,8 @@ def main():
         update_E(ee, hh, epsilon)
 
         if tt < src_maxt:
-            # This codebase uses E -= dt * J / epsilon for electric-current injection.
+            # 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)
 
@@ -193,9 +200,11 @@ def main():
             dt,
             ee,
             tt,
-            # The pulse is delayed relative to t=0, so the readout needs the same phase shift.
+            # 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 includes dt, so undo the dt in phasor_norm here.
+            # accumulate_phasor() already multiplies by dt, so pass the
+            # discrete-sum normalization without its extra dt factor.
             weight=phasor_norm / dt,
         )
 
@@ -205,6 +214,8 @@ def main():
     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}')

examples/waveguide.py

@@ -1,7 +1,11 @@
 """
-Example code for running an OpenCL FDTD simulation
+Example code for guided-wave FDFD and FDTD comparison.
 
-See main() for simulation setup.
+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
@@ -78,7 +82,7 @@ def get_waveguide_mode(
         omega: float,
         epsilon: fdfield_t,
         ) -> tuple[vcfdfield_t, vcfdfield_t]:
-    """ Create a mode source and overlap window """
+    """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),
@@ -94,6 +98,8 @@ def get_waveguide_mode(
     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
 
@@ -111,7 +117,8 @@ def main(
 
     grid, epsilon = draw_grid(dx=dx, pml_thickness=pml_thickness)
 
-    # Add PML
+    # 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):
@@ -275,7 +282,8 @@ def main2():
     # print(f'Save time interval would be {sample_interval} * dt = {sample_interval * dt:3g}')
 
 
-    # Source parameters and function
+    # 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
@@ -295,7 +303,8 @@ def main2():
         update_E(ee, hh, epsilon)
 
         if tt < src_maxt:
-            # This codebase uses E -= dt * J / epsilon for electric-current injection.
+            # 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)
 
@@ -317,9 +326,11 @@ def main2():
             dt,
             ee,
             tt,
-            # The pulse is delayed relative to t=0, so the readout needs the same phase shift.
+            # 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 includes dt, so undo the dt in phasor_norm here.
+            # accumulate_phasor() already contributes dt, so remove the extra dt
+            # from the externally computed normalization.
             weight=phasor_norm / dt,
         )
 
@@ -329,6 +340,8 @@ def main2():
     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}')

meanas/fdfd/__init__.py

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

meanas/fdfd/functional.py

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

meanas/fdfd/operators.py

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

meanas/fdfd/waveguide_2d.py

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

meanas/fdfd/waveguide_3d.py

@@ -3,6 +3,21 @@ 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
@@ -26,9 +41,9 @@ def solve_mode(
         epsilon: fdfield,
         mu: fdfield | None = None,
         ) -> dict[str, complex | NDArray[complexfloating]]:
-    """
+    r"""
     Given a 3D grid, selects a slice from the grid and attempts to
-     solve for an eigenmode propagating through that slice.
+    solve for an eigenmode propagating through that slice.
 
     Args:
         mode_number: Number of the mode, 0-indexed
@@ -37,19 +52,22 @@ def solve_mode(
         axis: Propagation axis (0=x, 1=y, 2=z)
         polarity: Propagation direction (+1 for +ve, -1 for -ve)
         slices: `epsilon[tuple(slices)]` is used to select the portion of the grid to use
-                as the waveguide cross-section. `slices[axis]` should select only one item.
+                as the waveguide cross-section. `slices[axis]` must select exactly one item.
         epsilon: Dielectric constant
         mu: Magnetic permeability (default 1 everywhere)
 
     Returns:
-        ```
+        Dictionary containing:
-        {
+
-            'E': NDArray[complexfloating],
+        - `E`: full-grid electric field for the solved mode
-            'H': NDArray[complexfloating],
+        - `H`: full-grid magnetic field for the solved mode
-            'wavenumber': complex,
+        - `wavenumber`: propagation constant corrected for the discretized
-            'wavenumber_2d': complex,
+          propagation axis
-        }
+        - `wavenumber_2d`: propagation constant of the reduced 2D eigenproblem
-        ```
+
+    Notes:
+        The returned fields are normalized through the `waveguide_2d`
+        normalization convention before being expanded back to 3D.
     """
     if mu is None:
         mu = numpy.ones_like(epsilon)
@@ -139,7 +157,14 @@ def compute_source(
         mu: Magnetic permeability (default 1 everywhere)
 
     Returns:
-        J distribution for the unidirectional source
+        `J` distribution for a one-sided electric-current source.
+
+    Notes:
+        The source is built from the expanded mode field and a boundary-source
+        operator. The resulting current is intended to be injected with the
+        same sign convention used elsewhere in the package:
+
+        `E -= dt * J / epsilon`
     """
     E_expanded = expand_e(E=E, dxes=dxes, wavenumber=wavenumber, axis=axis,
                           polarity=polarity, slices=slices)
@@ -167,10 +192,33 @@ def compute_overlap_e(
         slices: Sequence[slice],
         omega: float,
         ) -> cfdfield_t:
-    """
+    r"""
-    Given an eigenmode obtained by `solve_mode`, calculates an overlap_e for the
+    Build an overlap field for projecting another 3D electric field onto a mode.
-    mode orthogonality relation Integrate(((E x H_mode) + (E_mode x H)) dot dn)
+
-    [assumes reflection symmetry].
+    The returned field is intended for the discrete overlap expression
+
+    $$
+    \sum \mathrm{overlap\_e} \; E_\mathrm{other}^*
+    $$
+
+    where the sum is over the full Yee-grid field storage.
+
+    The construction uses a two-cell window immediately upstream of the selected
+    slice:
+
+    - for `polarity=+1`, the two cells just before `slices[axis].start`
+    - for `polarity=-1`, the two cells just after `slices[axis].stop`
+
+    The window is clipped to the simulation domain if necessary. A clipped but
+    non-empty window raises `RuntimeWarning`; an empty window raises
+    `ValueError`.
+
+    The derivation below assumes reflection symmetry and the standard waveguide
+    overlap relation involving
+
+    $$
+    \int ((E \times H_\mathrm{mode}) + (E_\mathrm{mode} \times H)) \cdot dn.
+    $$
 
     E x H_mode + E_mode x H
     -> Ex Hmy - EyHmx + Emx Hy - Emy Hx  (Z-prop)
@@ -183,8 +231,6 @@ def compute_overlap_e(
     B/wu (Ex Emx + Ey Emy) - j/wu (Ex dx Emz + Ey dy Emz)
 
 
-    TODO: add reference
-
     Args:
         E: E-field of the mode
         wavenumber: Wavenumber of the mode
@@ -195,7 +241,8 @@ def compute_overlap_e(
                 as the waveguide cross-section. slices[axis] should select only one item.
 
     Returns:
-        overlap_e such that `numpy.sum(overlap_e * other_e.conj())` computes the overlap integral
+        `overlap_e` normalized so that `numpy.sum(overlap_e * E.conj()) == 1`
+        over the retained overlap window.
     """
     slices = tuple(slices)
 
@@ -240,7 +287,7 @@ def expand_e(
         polarity: int,
         slices: Sequence[slice],
         ) -> cfdfield_t:
-    """
+    r"""
     Given an eigenmode obtained by `solve_mode`, expands the E-field from the 2D
     slice where the mode was calculated to the entire domain (along the propagation
     axis). This assumes the epsilon cross-section remains constant throughout the
@@ -258,6 +305,16 @@ def expand_e(
 
     Returns:
         `E`, with the original field expanded along the specified `axis`.
+
+    Notes:
+        This helper assumes that the waveguide cross-section remains constant
+        along the propagation axis and applies the phase factor
+
+        $$
+        e^{-i \, \mathrm{polarity} \, wavenumber \, \Delta z}
+        $$
+
+        to each copied slice.
     """
     slices = tuple(slices)
 

meanas/fdfd/waveguide_cyl.py

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

meanas/fdtd/__init__.py

@@ -168,7 +168,44 @@ t=0 (assuming the source is off for t<0 this gives $\sim 10^{-3}$ error at t=0).
 
 Boundary conditions
 ===================
-# TODO notes about boundaries / PMLs
+
+`meanas.fdtd` exposes two boundary-related building blocks:
+
+- `conducting_boundary(...)` for simple perfect-electric-conductor style field
+  clamping at one face of the domain.
+- `cpml_params(...)` and `updates_with_cpml(...)` for convolutional perfectly
+  matched layers (CPMLs) attached to one or more faces of the Yee grid.
+
+`updates_with_cpml(...)` accepts a three-by-two table of CPML parameter blocks:
+
+```
+cpml_params[axis][polarity_index]
+```
+
+where `axis` is `0`, `1`, or `2` and `polarity_index` corresponds to `(-1, +1)`.
+Passing `None` for one entry disables CPML on that face while leaving the other
+directions unchanged. This is how mixed boundary setups such as "absorbing in x,
+periodic in y/z" are expressed.
+
+When comparing an FDTD run against an FDFD solve, use the same stretched
+coordinate system in both places:
+
+1. Build the FDTD update with the desired CPML parameters.
+2. Stretch the FDFD `dxes` with the matching SCPML transform.
+3. Compare the extracted phasor against the FDFD field or residual on those
+   stretched `dxes`.
+
+The electric-current sign convention used throughout the examples and tests is
+
+$$
+E \leftarrow E - \Delta_t J / \epsilon
+$$
+
+which matches the FDFD right-hand side
+
+$$
+-i \omega J.
+$$
 """
 
 from .base import (

meanas/fdtd/energy.py

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

meanas/fdtd/pml.py

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