[docs] document scattering conventions and frozen-basis workflows

Update examples for sampled-interface tapers, explicit physical CPML spacing, copied SCPML widths, and the consolidated complex mode solver. Document breaking API changes, overlap measures, and PML convergence.

README.md

@@ -49,7 +49,7 @@ linear systems, ideally with double precision.
 
 **Requirements:**
 
-* python >=3.11
+* python >=3.12
 * numpy
 * scipy
 
@@ -118,11 +118,12 @@ python3 -m pytest -q | tee test_results.txt
 - `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.fdfd.waveguide_2d_complex`: targeted signed-`beta` waveguide modes,
-  transpose-left / Lorentz-dual projection fields, PML-mode diagnostics, and
-  overlap-based tracking for open cross-sections.
+  solvers, including targeted complex bases from the default half-sized
+  beta-squared or explicit signed-beta formulations, transpose-left and
+  Lorentz-dual projection fields on a shared field-array axis, separately
+  returned PML-mode diagnostics, overlap-based tracking,
+  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,
@@ -175,7 +176,7 @@ The tracked examples under `examples/` are the intended entry points for users:
 - `examples/eme_bend.py`: straight-to-bent waveguide mode matching with
   cylindrical bend modes, interface scattering, and a cascaded bend-network
   example built with `scikit-rf`.
-- `examples/eme_taper_cmt.py`: sampled cross-section local-mode CMT for a
+- `examples/eme_taper_cmt.py`: sampled-interface EME with stable S cascading for a
   continuously varying rib-waveguide taper.
 - `examples/eme_taper_pml_modes.py`: targeted complex/PML mode solving,
   tracking, and construction of EME taper sections.

docs/api/fdfd.md

@@ -13,3 +13,17 @@
 ::: meanas.fdfd.scpml
 
 ::: meanas.fdfd.farfield
+
+## SCPML conventions
+
+`stretch_with_scpml` now uses `-1` for the low-index face and `+1` for the
+high-index face, matching CPML and propagation conventions. It returns copied
+widths by default; pass `copy=False` to request mutation explicitly.
+
+Both SCPML builders accept `kappa_function(depth)` for real stretching and
+`max_imag_ratio` to limit the final dimensionless `abs(Im(s))/Re(s)`. Start a
+mode-solver convergence study with `max_imag_ratio=1.0`: overly strong imaginary
+stretching can produce above-light-line PML artifacts. The cap is not a general
+spectral guarantee, and weaker absorption can increase boundary reflections.
+Check domain size, PML thickness, strength, mode localization, and residuals.
+The general scattering profile remains uncapped by default.

docs/api/fdtd.md

@@ -13,3 +13,8 @@
 ::: meanas.fdtd.energy
 
 ::: meanas.fdtd.phasor
+
+CPML derivatives use electric widths for forward differences and magnetic
+widths for backward differences. `cpml_params(..., dx=1.0)` specifies the uniform
+physical cell width normal to a face; arbitrary nonuniform PML grading is not
+implemented by this constructor. Scale the timestep consistently with the grid.

docs/api/waveguides.md

@@ -2,8 +2,90 @@
 
 ::: meanas.fdfd.waveguide_2d
 
-::: meanas.fdfd.waveguide_2d_complex
-
 ::: meanas.fdfd.waveguide_3d
 
 ::: meanas.fdfd.waveguide_cyl
+
+## Overlap measures and normalization
+
+`waveguide_2d.inner_product(..., measure="physical", conj_h=True)` uses real
+Yee widths for physical Poynting power. `measure="stretched"` keeps complex
+widths for coordinate reaction. `lorentz_metric` and `lorentz_overlap` default
+to the stretched measure. Use the same grid and measure when constructing and
+using dual fields. A real overlap grid still defines a consistent constructed
+dual basis, even when the eigenproblem uses a separate PML grid.
+
+`normalized_fields_e/h(..., normalization="lorentz")` normalizes unconjugated
+self-reaction to one. This fixes phase up to a sign and is distinct from unit
+forward power, which remains the default. Zero-power or self-orthogonal modes
+raise `ValueError`; open self-orthogonal modes require paired primal/dual
+normalization through `ModeBasis`.
+
+## Targeted and reproducible mode solving
+
+For modes near an effective index guess, use `solve_modes_near` with
+`target_wavenumbers=[omega * n_guess]`. It targets the electric beta-squared
+operator at `(omega * n_guess)**2`; `solve_modes_near_beta` instead targets the
+first-order signed-beta operator. Both use the LU shift-invert backend and
+support `eigensolver_opts={"v0": starting_vector}`. Factorization cost and
+memory depend on the grid, so this is not universally faster than indexed
+solving. `solve_modes` keeps its global mode-number semantics and now accepts
+`rng` or `v0` for reproducible initialization.
+
+`waveguide_3d.mode_from_2d` accepts a supplied transverse electric eigenvector
+and reuses field reconstruction, normalization, numerical-dispersion
+correction, Yee staggering, rotation, and full-grid embedding without solving
+an eigenproblem.
+
+## Interface and sampled-taper scattering
+
+All scattering matrices use outgoing rows and incoming columns, ordered
+`[left, right]`:
+
+```text
+[b_left ]   [R_left        T_right_to_left] [a_left ]
+[b_right] = [T_left_to_right       R_right] [a_right]
+```
+
+`get_s` and `get_abcd` require explicit `dxes`. Mixed Petrov-Galerkin testing is
+the default: left electric duals test H continuity, and right magnetic duals
+test E continuity. Both incidence directions are solved together with full
+projection matrices. Mixed scattering supports unequal mode counts.
+`formulation="galerkin"` selects equal-count left testing of both equations.
+Neither formulation silently clips transmission. The optional `force_nogain`
+and `force_reciprocal` flags are numerical projections, not general physical
+tests for arbitrary open or evanescent modal coordinates.
+
+`get_taper_s` replaces finite-difference CMT with half propagation in the left
+basis, a midpoint interface, and half propagation in the right basis for each
+sampled interval. Redheffer composition avoids growing passive backward
+transfer branches. Refine section sampling and modal/PML bases to establish
+convergence. The ineffective `rtol`, `atol`, and `max_step` arguments have been
+removed. The existing `eme_taper_cmt.py` example filename is retained.
+
+```python
+prepared = eme.prepare_taper(sections, dxes)
+s_original = eme.evaluate_taper_s(prepared)
+s_longer = eme.evaluate_taper_s(prepared, z_coords=prepared.z_coords * 1.2)
+s_changed_beta = eme.evaluate_taper_s(
+    prepared, wavenumber_matrices=beta_matrices,
+)
+```
+
+`PreparedTaper` owns read-only copies of phase-aligned interface data. Beta
+matrix overrides have shape `(section, mode, mode)` in the original supplied
+section bases; recorded basis transforms are applied automatically. Explicit
+backward branches require backward beta overrides when forward beta is changed.
+Evaluation performs no field overlaps or eigensolves. Length sweeps reuse the
+sampled model; wavelength sweeps with reused fields are explicitly frozen-basis
+approximations. Reprepare when fields, duals, grids, or formulations change.
+
+`cascade_s(left, right, internal_modes=...)` connects spatially ordered S
+matrices, including unequal external ports. Public `abcd_to_s` and `s_to_abcd`
+use `[f_right, b_right] = ABCD @ [f_left, b_left]` and require equal port counts.
+`get_taper_abcd` converts the stable S result; zero transmission or an
+unrepresentable transfer matrix raises `LinAlgError`. Prefer S for long lossy
+devices. Singular interface and cascade systems also raise contextual errors
+rather than silently substituting pseudoinverses.
+
+::: meanas.fdfd.eme

docs/index.md

@@ -28,7 +28,7 @@ Relevant starting examples:
   scattering between two nearby waveguide cross-sections
 - `examples/eme_bend.py` for straight-to-bent mode matching with cylindrical
   bend modes and a cascaded bend-network example
-- `examples/eme_taper_cmt.py` for local-mode CMT through sampled continuously
+- `examples/eme_taper_cmt.py` for sampled-interface EME through continuously
   varying taper cross-sections
 - `examples/fdfd.py` for direct frequency-domain waveguide excitation
 

examples/eme_taper_cmt.py

@@ -1,5 +1,5 @@
 """
-Local-mode CMT example for a continuously varying rib-waveguide taper.
+Sampled-interface EME example for a continuously varying rib-waveguide taper.
 
 This example keeps geometry construction outside `meanas.fdfd.eme`: it samples a
 width taper at several cross-sections, solves and normalizes each local mode with
@@ -112,14 +112,15 @@ def print_summary(
     print('sampled taper effective indices:', ', '.join(f'{value:.5f}' for value in neffs))
     print(f'abrupt endpoint reflection       |S_00|^2          = {abs(abrupt_s[0, 0]) ** 2:.6f}')
     print(f'abrupt endpoint transmission     |S_{n_modes},0|^2 = {abs(abrupt_s[n_modes, 0]) ** 2:.6f}')
-    print(f'taper CMT reflection             |S_00|^2          = {abs(taper_s[0, 0]) ** 2:.6f}')
-    print(f'taper CMT transmission           |S_{n_modes},0|^2 = {abs(taper_s[n_modes, 0]) ** 2:.6f}')
-    print(f'taper CMT total output power                       = {numpy.sum(abs(taper_s[:, 0]) ** 2):.6f}')
+    print(f'taper reflection             |S_00|^2          = {abs(taper_s[0, 0]) ** 2:.6f}')
+    print(f'taper transmission           |S_{n_modes},0|^2 = {abs(taper_s[n_modes, 0]) ** 2:.6f}')
+    print(f'taper total output power                       = {numpy.sum(abs(taper_s[:, 0]) ** 2):.6f}')
 
 
 def main() -> None:
     sections, neffs, _omega, dxes_2d = solve_taper_sections()
-    taper_s = eme.get_taper_s(sections, dxes=dxes_2d)
+    prepared = eme.prepare_taper(sections, dxes=dxes_2d)
+    taper_s = eme.evaluate_taper_s(prepared)
     abrupt_s = eme.get_s(
         sections[0].modes,
         sections[0].wavenumbers,

examples/eme_taper_pml_modes.py

@@ -5,7 +5,7 @@ from __future__ import annotations
 import numpy
 
 from meanas.fdmath import vec
-from meanas.fdfd import eme, scpml, waveguide_2d_complex
+from meanas.fdfd import eme, scpml, waveguide_2d
 
 
 WAVELENGTH = 1.55
@@ -34,6 +34,7 @@ def build_grids() -> tuple[list[list[numpy.ndarray]], list[list[numpy.ndarray]]]
                 omega=OMEGA,
                 epsilon_effective=EPS_CLADDING,
                 thickness=PML_CELLS,
+                max_imag_ratio=1.0,
             )
     return physical, stretched
 
@@ -61,29 +62,48 @@ def main() -> None:
     """Solve, track, and summarize the sampled complex modal bases."""
     physical_dxes, solve_dxes = build_grids()
     z_coords = numpy.array([0.0, 5.0, 10.0])
-    bases = [
-        waveguide_2d_complex.solve_modes_near(
+    epsilon_samples = [build_epsilon(edge_fill) for edge_fill in (0.30, 0.35, 0.40)]
+    solved = [
+        waveguide_2d.solve_modes_near(
             OMEGA,
             solve_dxes,
-            build_epsilon(edge_fill),
+            epsilon,
             target_wavenumbers=[1.45 * OMEGA],
             modes_per_target=3,
             candidate_margin=3,
             overlap_dxes=physical_dxes,
         )
-        for edge_fill in (0.30, 0.35, 0.40)
+        for epsilon in epsilon_samples
+    ]
+    bases = [basis for basis, _diagnostics in solved]
+    tracked = waveguide_2d.track_mode_path(bases, overlap_dxes=physical_dxes)
+    tracked_diagnostics = [
+        waveguide_2d.compute_mode_diagnostics(
+            basis,
+            OMEGA,
+            solve_dxes,
+            epsilon,
+            target_wavenumbers=[1.45 * OMEGA],
+            overlap_dxes=physical_dxes,
+        )
+        for basis, epsilon in zip(tracked, epsilon_samples, strict=True)
     ]
-    tracked = waveguide_2d_complex.track_mode_path(bases)
     sections = eme.mode_sections_from_bases(z_coords, tracked)
+    prepared = eme.prepare_taper(sections, physical_dxes)
+    # Fields and interface overlaps are reused: these are frozen-basis length
+    # sweeps. Wavelength-dependent beta overrides would also freeze the fields.
+    length_sweep = [eme.evaluate_taper_s(prepared, z_coords=z_coords * scale)
+                    for scale in (0.8, 1.0, 1.2)]
 
-    for z_coord, basis in zip(z_coords, tracked, strict=True):
+    for z_coord, basis, diagnostics in zip(z_coords, tracked, tracked_diagnostics, strict=True):
         print(
             f'z={z_coord:5.1f}: '
             f'beta={basis.wavenumbers[0]:.6g}, '
-            f'PML participation={basis.diagnostics.pml_participation[0]:.3f}, '
-            f'right residual={basis.diagnostics.right_residuals[0]:.2e}',
+            f'PML participation={diagnostics.pml_participation[0]:.3f}, '
+            f'right residual={diagnostics.right_residuals[0]:.2e}',
         )
     print(f'constructed {len(sections)} tracked EME sections with {len(sections[0].modes)} modes each')
+    print('frozen-basis length sweep |T_00|:', [float(abs(ss[len(tracked[0]), 0])) for ss in length_sweep])
 
 
 if __name__ == '__main__':

examples/fdtd.py

@@ -145,7 +145,7 @@ def main() -> None:
 
     # PMLs in every direction
     pml_params = [
-        [cpml_params(axis=dd, polarity=pp, dt=dt, thickness=pml_thickness, epsilon_eff=n_air ** 2)
+        [cpml_params(axis=dd, polarity=pp, dt=dt, dx=dx, thickness=pml_thickness, epsilon_eff=n_air ** 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, dtype=complex)
@@ -203,7 +203,7 @@ def main() -> None:
     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)
+            dxes_fdfd = 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)

examples/waveguide.py

@@ -273,7 +273,7 @@ def main2():
 
     # PMLs in every direction
     pml_params = [
-        [cpml_params(axis=dd, polarity=pp, dt=dt, thickness=pml_thickness, epsilon_eff=n_cladding ** 2)
+        [cpml_params(axis=dd, polarity=pp, dt=dt, dx=dx, 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, dtype=complex)
@@ -330,7 +330,7 @@ def main2():
     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)
+            dxes_fdfd = 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)

meanas/test/test_examples_smoke.py

@@ -57,4 +57,4 @@ def test_eme_taper_cmt_example_smoke_runs(tmp_path: Path) -> None:
 
     assert result.returncode == 0, result.stdout + result.stderr
     assert 'sampled taper effective indices:' in result.stdout
-    assert 'taper CMT transmission' in result.stdout
+    assert 'taper transmission' in result.stdout