From ff278e6fa174305c86d8e65c390265c2afde398d Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Sun, 19 Apr 2026 16:57:22 -0700 Subject: [PATCH 1/4] [docs] more docs cleanup --- meanas/fdfd/operators.py | 29 ++++++++++++++--------------- meanas/fdfd/waveguide_3d.py | 4 ++-- meanas/fdfd/waveguide_cyl.py | 6 +++--- meanas/fdmath/operators.py | 2 +- meanas/fdmath/vectorization.py | 3 +-- pyproject.toml | 1 + 6 files changed, 22 insertions(+), 23 deletions(-) diff --git a/meanas/fdfd/operators.py b/meanas/fdfd/operators.py index 18aaade..6425a29 100644 --- a/meanas/fdfd/operators.py +++ b/meanas/fdfd/operators.py @@ -64,11 +64,11 @@ def e_full( epsilon: Vectorized dielectric constant mu: Vectorized magnetic permeability (default 1 everywhere). pec: Vectorized mask specifying PEC cells. Any cells where `pec != 0` are interpreted - as containing a perfect electrical conductor (PEC). - The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`) + as containing a perfect electrical conductor (PEC). + The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`) pmc: Vectorized mask specifying PMC cells. Any cells where `pmc != 0` are interpreted - as containing a perfect magnetic conductor (PMC). - The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) + as containing a perfect magnetic conductor (PMC). + The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) Returns: Sparse matrix containing the wave operator. @@ -148,11 +148,11 @@ def h_full( epsilon: Vectorized dielectric constant mu: Vectorized magnetic permeability (default 1 everywhere) pec: Vectorized mask specifying PEC cells. Any cells where `pec != 0` are interpreted - as containing a perfect electrical conductor (PEC). - The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`) + as containing a perfect electrical conductor (PEC). + The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`) pmc: Vectorized mask specifying PMC cells. Any cells where `pmc != 0` are interpreted - as containing a perfect magnetic conductor (PMC). - The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) + as containing a perfect magnetic conductor (PMC). + The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) Returns: Sparse matrix containing the wave operator. @@ -217,11 +217,11 @@ def eh_full( epsilon: Vectorized dielectric constant mu: Vectorized magnetic permeability (default 1 everywhere) pec: Vectorized mask specifying PEC cells. Any cells where `pec != 0` are interpreted - as containing a perfect electrical conductor (PEC). - The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`) + as containing a perfect electrical conductor (PEC). + The PEC is applied per-field-component (i.e. `pec.size == epsilon.size`) pmc: Vectorized mask specifying PMC cells. Any cells where `pmc != 0` are interpreted - as containing a perfect magnetic conductor (PMC). - The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) + as containing a perfect magnetic conductor (PMC). + The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) Returns: Sparse matrix containing the wave operator. @@ -267,8 +267,8 @@ def e2h( dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` mu: Vectorized magnetic permeability (default 1 everywhere) pmc: Vectorized mask specifying PMC cells. Any cells where `pmc != 0` are interpreted - as containing a perfect magnetic conductor (PMC). - The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) + as containing a perfect magnetic conductor (PMC). + The PMC is applied per-field-component (i.e. `pmc.size == epsilon.size`) Returns: Sparse matrix for converting E to H. @@ -483,4 +483,3 @@ def e_boundary_source( # (numpy.roll(mask, +1, axis=2) != mask)) return sparse.diags_array(jmask.astype(int)) @ full - diff --git a/meanas/fdfd/waveguide_3d.py b/meanas/fdfd/waveguide_3d.py index 66cc7cc..e7dfd22 100644 --- a/meanas/fdfd/waveguide_3d.py +++ b/meanas/fdfd/waveguide_3d.py @@ -52,7 +52,7 @@ 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]` must select exactly one item. epsilon: Dielectric constant mu: Magnetic permeability (default 1 everywhere) @@ -62,7 +62,7 @@ def solve_mode( - `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 + propagation axis - `wavenumber_2d`: propagation constant of the reduced 2D eigenproblem Notes: diff --git a/meanas/fdfd/waveguide_cyl.py b/meanas/fdfd/waveguide_cyl.py index caedfaf..201f709 100644 --- a/meanas/fdfd/waveguide_cyl.py +++ b/meanas/fdfd/waveguide_cyl.py @@ -216,13 +216,13 @@ def solve_modes( of the bent waveguide with the specified mode number. Args: - mode_number: Number of the mode, 0-indexed + mode_numbers: Mode numbers to solve, 0-indexed. omega: Angular frequency of the simulation dxes: Grid parameters [dx_e, dx_h] as described in meanas.fdmath.types. - The first coordinate is assumed to be r, the second is y. + The first coordinate is assumed to be r, the second is y. epsilon: Dielectric constant rmin: Radius of curvature for the simulation. This should be the minimum value of - r within the simulation domain. + r within the simulation domain. Returns: e_xys: NDArray of vfdfield_t specifying fields. First dimension is mode number. diff --git a/meanas/fdmath/operators.py b/meanas/fdmath/operators.py index 19ccb80..0c64ae7 100644 --- a/meanas/fdmath/operators.py +++ b/meanas/fdmath/operators.py @@ -158,7 +158,7 @@ def cross( Args: B: List `[Bx, By, Bz]` of sparse matrices corresponding to the x, y, z - portions of the operator on the left side of the cross product. + portions of the operator on the left side of the cross product. Returns: Sparse matrix corresponding to (B x), where x is the cross product. diff --git a/meanas/fdmath/vectorization.py b/meanas/fdmath/vectorization.py index 44d8b74..77b722f 100644 --- a/meanas/fdmath/vectorization.py +++ b/meanas/fdmath/vectorization.py @@ -58,7 +58,7 @@ def vec( Args: f: A vector field, e.g. `[f_x, f_y, f_z]` where each `f_` component is a 1- to - 3-D ndarray (`f_*` should all be the same size). Doesn't fail with `f=None`. + 3-D ndarray (`f_*` should all be the same size). Doesn't fail with `f=None`. Returns: 1D ndarray containing the linearized field (or `None`) @@ -123,4 +123,3 @@ def unvec( if v is None: return None return v.reshape((nvdim, *shape), order='C') # type: ignore - diff --git a/pyproject.toml b/pyproject.toml index 4b81fa6..20a694f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -63,6 +63,7 @@ docs = [ "mkdocs-print-site-plugin>=2.3", "pymdown-extensions>=10.7", "htmlark>=1.0", + "ruff>=0.6", ] examples = [ "matplotlib>=3.10.8", From bb920b8e33466dc9c000dd897714cf0de9b46c98 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Sun, 19 Apr 2026 17:30:04 -0700 Subject: [PATCH 2/4] [tests / fdtd.pml] add pml test --- meanas/test/test_fdtd_pml.py | 202 +++++++++++++++++++++++++++++++++++ 1 file changed, 202 insertions(+) diff --git a/meanas/test/test_fdtd_pml.py b/meanas/test/test_fdtd_pml.py index 6d118c6..9d8aef8 100644 --- a/meanas/test/test_fdtd_pml.py +++ b/meanas/test/test_fdtd_pml.py @@ -1,7 +1,10 @@ import numpy import pytest +from .. import fdtd +from ..fdtd.base import maxwell_e, maxwell_h from ..fdtd.pml import cpml_params, updates_with_cpml +from .utils import assert_close @pytest.mark.parametrize( @@ -42,3 +45,202 @@ def test_updates_with_cpml_keeps_zero_fields_zero() -> None: assert not e.any() assert not h.any() + + +def _unit_dxes(shape: tuple[int, int, int, int]) -> list[list[numpy.ndarray]]: + return [[numpy.ones(n, dtype=float) for n in shape[1:]] for _ in range(2)] + + +def _real_field(shape: tuple[int, int, int, int], start: float) -> numpy.ndarray: + total = numpy.prod(shape, dtype=int) + return numpy.arange(start, start + total, dtype=float).reshape(shape) / total + + +def _complex_field(shape: tuple[int, int, int, int], start: float) -> numpy.ndarray: + real = _real_field(shape, start) + imag = _real_field(shape, start + real.size) + return real + 1j * imag + + +def test_updates_with_cpml_matches_base_updates_when_all_faces_disabled() -> None: + shape = (3, 4, 5, 6) + epsilon = _real_field(shape, 1.0) + 2.0 + mu = _real_field(shape, 4.0) + 1.5 + e = _real_field(shape, 10.0) + h = _real_field(shape, 100.0) + dxes = _unit_dxes(shape) + params = [[None, None] for _ in range(3)] + + update_e_cpml, update_h_cpml = updates_with_cpml(params, dt=0.1, dxes=dxes, epsilon=epsilon) + update_e_base = maxwell_e(dt=0.1, dxes=dxes) + update_h_base = maxwell_h(dt=0.1, dxes=dxes) + + e_cpml = e.copy() + h_cpml = h.copy() + e_base = e.copy() + h_base = h.copy() + + update_e_cpml(e_cpml, h_cpml, epsilon) + update_e_base(e_base, h_base, epsilon) + update_h_cpml(e_cpml, h_cpml, mu) + update_h_base(e_base, h_base, mu) + + assert_close(e_cpml, e_base) + assert_close(h_cpml, h_base) + + +def test_updates_with_cpml_matches_base_updates_with_complex_dtype_when_all_faces_disabled() -> None: + shape = (3, 4, 5, 6) + epsilon = _real_field(shape, 1.0) + 2.0 + mu = _real_field(shape, 4.0) + 1.5 + e = _complex_field(shape, 10.0) + h = _complex_field(shape, 100.0) + dxes = _unit_dxes(shape) + params = [[None, None] for _ in range(3)] + + update_e_cpml, update_h_cpml = updates_with_cpml(params, dt=0.1, dxes=dxes, epsilon=epsilon, dtype=complex) + update_e_base = maxwell_e(dt=0.1, dxes=dxes) + update_h_base = maxwell_h(dt=0.1, dxes=dxes) + + e_cpml = e.copy() + h_cpml = h.copy() + e_base = e.copy() + h_base = h.copy() + + update_e_cpml(e_cpml, h_cpml, epsilon) + update_e_base(e_base, h_base, epsilon) + update_h_cpml(e_cpml, h_cpml, mu) + update_h_base(e_base, h_base, mu) + + assert_close(e_cpml, e_base) + assert_close(h_cpml, h_base) + + +def test_updates_with_cpml_only_changes_the_configured_face_region() -> None: + shape = (3, 6, 6, 6) + epsilon = numpy.ones(shape, dtype=float) + mu = numpy.ones(shape, dtype=float) + e = _real_field(shape, 1.0) + h = _real_field(shape, 100.0) + dxes = _unit_dxes(shape) + thickness = 3 + + params = [[None, None] for _ in range(3)] + params[0][0] = cpml_params(axis=0, polarity=-1, dt=0.1, thickness=thickness) + + update_e_cpml, update_h_cpml = updates_with_cpml(params, dt=0.1, dxes=dxes, epsilon=epsilon) + update_e_base = maxwell_e(dt=0.1, dxes=dxes) + update_h_base = maxwell_h(dt=0.1, dxes=dxes) + + e_cpml = e.copy() + h_cpml = h.copy() + e_base = e.copy() + h_base = h.copy() + + update_e_cpml(e_cpml, h_cpml, epsilon) + update_e_base(e_base, h_base, epsilon) + update_h_cpml(e_cpml, h_cpml, mu) + update_h_base(e_base, h_base, mu) + + e_untouched = slice(thickness, None) + h_untouched = slice(thickness, -1) + assert_close(e_cpml[:, e_untouched, :, :], e_base[:, e_untouched, :, :]) + assert_close(h_cpml[:, h_untouched, :, :], h_base[:, h_untouched, :, :]) + + changed_e = numpy.any(numpy.abs(e_cpml[:, :thickness, :, :] - e_base[:, :thickness, :, :]) > 1e-12) + changed_h = numpy.any(numpy.abs(h_cpml[:, :thickness, :, :] - h_base[:, :thickness, :, :]) > 1e-12) + assert changed_e + assert changed_h + + +def test_cpml_plane_wave_phasor_decays_monotonically_through_outgoing_pml() -> None: + dt = 0.4 + period_steps = 24 + omega = 2 * numpy.pi / (period_steps * dt) + shape = (3, 80, 1, 1) + thickness = 8 + source_x = 16 + warmup_periods = 10 + accumulation_periods = 6 + total_steps = period_steps * (warmup_periods + accumulation_periods) + + epsilon = numpy.ones(shape, dtype=float) + dxes = _unit_dxes(shape) + params = [[None, None] for _ in range(3)] + for polarity_index, polarity in enumerate((-1, 1)): + params[0][polarity_index] = cpml_params(axis=0, polarity=polarity, dt=dt, thickness=thickness) + + update_e, update_h = updates_with_cpml(params, dt=dt, dxes=dxes, epsilon=epsilon) + + e = numpy.zeros(shape, dtype=float) + h = numpy.zeros(shape, dtype=float) + e_accumulator = numpy.zeros((1, *shape), dtype=complex) + + for step in range(total_steps): + update_e(e, h, epsilon) + + source = numpy.cos(omega * (step + 0.5) * dt) + e[1, source_x, 0, 0] -= dt * source + + if step >= period_steps * warmup_periods: + fdtd.accumulate_phasor_e(e_accumulator, omega, dt, e, step + 1) + + update_h(e, h) + + profile = numpy.abs(e_accumulator[0, 1, :, 0, 0]) + right_pml = profile[-thickness:] + interior = profile[-thickness - 6:-thickness] + interior_level = interior.mean() + + assert interior_level > 1.0 + assert right_pml[-1] < interior_level / 100 + assert profile[0] < interior_level / 100 + assert numpy.all(numpy.diff(right_pml) <= interior_level * 1e-3) + + +def test_cpml_point_source_total_energy_reaches_late_time_plateau() -> None: + dt = 0.3 + period_steps = 24 + omega = 2 * numpy.pi / (period_steps * dt) + cycles = 1000 + sample_every_cycles = 50 + sample_stride = period_steps * sample_every_cycles + shape = (3, 9, 9, 9) + thickness = 3 + center = shape[1] // 2 + + epsilon = numpy.ones(shape, dtype=float) + dxes = _unit_dxes(shape) + params = [[None, None] for _ in range(3)] + for axis in range(3): + for polarity_index, polarity in enumerate((-1, 1)): + params[axis][polarity_index] = cpml_params(axis=axis, polarity=polarity, dt=dt, thickness=thickness) + + update_e, update_h = updates_with_cpml(params, dt=dt, dxes=dxes, epsilon=epsilon) + + e = numpy.zeros(shape, dtype=float) + h = numpy.zeros(shape, dtype=float) + sampled_energies: list[float] = [] + + for step in range(period_steps * cycles): + h_before = h.copy() + update_e(e, h, epsilon) + + source = numpy.cos(omega * (step + 0.5) * dt) + e[1, center, center, center] -= dt * source + + update_h(e, h) + + if (step + 1) % sample_stride == 0: + total_energy = fdtd.energy_estep(h0=h_before, e1=e, h2=h, epsilon=epsilon, dxes=dxes).sum().real + sampled_energies.append(total_energy) + + energies = numpy.asarray(sampled_energies) + late_window = energies[-5:] + previous_window = energies[-10:-5] + late_mean = late_window.mean() + + assert energies.size == cycles // sample_every_cycles + assert late_mean > 0.1 + assert (late_window.max() - late_window.min()) / late_mean < 1e-4 + assert abs(late_mean - previous_window.mean()) / late_mean < 1e-4 From 9a0c693848842958c40c8afccb40e377262c128e Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Sun, 19 Apr 2026 20:22:36 -0700 Subject: [PATCH 3/4] [docs] docs dark mode --- docs/stylesheets/extra.css | 39 ++++++++++++++++++++++++++++++++++++++ mkdocs.yml | 13 +++++++++++++ 2 files changed, 52 insertions(+) diff --git a/docs/stylesheets/extra.css b/docs/stylesheets/extra.css index d70a077..090b5e7 100644 --- a/docs/stylesheets/extra.css +++ b/docs/stylesheets/extra.css @@ -11,3 +11,42 @@ .md-typeset h3 code { word-break: break-word; } + +[data-md-color-scheme="slate"] { + --md-default-bg-color: #0f141c; + --md-default-fg-color: #e8eef7; + --md-default-fg-color--light: #b3bfd1; + --md-default-fg-color--lighter: #7f8ba0; + --md-default-fg-color--lightest: #5d6880; + --md-code-bg-color: #111923; + --md-code-fg-color: #e4edf8; + --md-accent-fg-color: #7dd3fc; +} + +[data-md-color-scheme="slate"] .md-header, +[data-md-color-scheme="slate"] .md-tabs { + background: linear-gradient(90deg, #111923 0%, #162235 100%); +} + +[data-md-color-scheme="slate"] .md-typeset pre > code, +[data-md-color-scheme="slate"] .md-typeset code { + border: 1px solid rgba(125, 211, 252, 0.14); +} + +[data-md-color-scheme="slate"] .md-typeset table:not([class]) { + background: rgba(255, 255, 255, 0.015); +} + +[data-md-color-scheme="slate"] .md-typeset table:not([class]) th { + background: rgba(125, 211, 252, 0.08); +} + +[data-md-color-scheme="slate"] .md-typeset .admonition, +[data-md-color-scheme="slate"] .md-typeset details { + background: rgba(255, 255, 255, 0.02); + border-color: rgba(125, 211, 252, 0.2); +} + +[data-md-color-scheme="slate"] .md-typeset .arithmatex { + padding: 0.1rem 0; +} diff --git a/mkdocs.yml b/mkdocs.yml index a0dcd2c..379af9b 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -10,6 +10,19 @@ strict: false theme: name: material font: false + palette: + - scheme: slate + primary: blue grey + accent: cyan + toggle: + icon: material/weather-sunny + name: Switch to light mode + - scheme: default + primary: teal + accent: indigo + toggle: + icon: material/weather-night + name: Switch to dark mode features: - navigation.indexes - navigation.sections From f8ad0250d11db85d607c41b0ec18ca2b85f4a182 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Mon, 20 Apr 2026 10:15:25 -0700 Subject: [PATCH 4/4] [eme / examples] add EME examples --- README.md | 5 + docs/index.md | 4 + examples/eme.py | 217 ++++++++++++++++++++++++++++++ examples/eme_bend.py | 310 +++++++++++++++++++++++++++++++++++++++++++ pyproject.toml | 1 + 5 files changed, 537 insertions(+) create mode 100644 examples/eme.py create mode 100644 examples/eme_bend.py diff --git a/README.md b/README.md index ce9bcc1..b400796 100644 --- a/README.md +++ b/README.md @@ -163,6 +163,11 @@ The tracked examples under `examples/` are the intended entry points for users: guide, with late-time monitor slices, guided-core windows, and mode-weighted errors compared directly against real fields reconstructed from the matching FDFD solution, plus a guided-mode / orthogonal-residual split. +- `examples/eme.py`: straight-interface mode matching / EME, including port + mode solving, interface scattering, and modal field visualization. +- `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/fdfd.py`: direct frequency-domain waveguide excitation and overlap / Poynting analysis without a time-domain run. diff --git a/docs/index.md b/docs/index.md index bc4cdeb..5475af8 100644 --- a/docs/index.md +++ b/docs/index.md @@ -24,6 +24,10 @@ Relevant starting examples: - `examples/waveguide_real.py` for real-valued continuous-wave FDTD compared against real fields reconstructed from an FDFD solution, including guided-core, mode-weighted, and guided-mode / residual comparisons +- `examples/eme.py` for straight-interface mode matching / EME and modal + 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/fdfd.py` for direct frequency-domain waveguide excitation For solver equivalence, prefer the phasor-based examples first. They compare diff --git a/examples/eme.py b/examples/eme.py new file mode 100644 index 0000000..6215cbc --- /dev/null +++ b/examples/eme.py @@ -0,0 +1,217 @@ +""" +Mode-matching / EME example for a straight rib-waveguide interface. + +This example shows the intended user-facing workflow for `meanas.fdfd.eme` on a +simple straight interface: + +1. build two nearby waveguide cross-sections on a Yee grid, +2. solve a small set of guided modes on each side, +3. normalize those modes into E/H port fields, +4. assemble the interface scattering matrix with `meanas.fdfd.eme.get_s(...)`, +5. inspect the dominant modal coupling numerically and visually. +""" + +from __future__ import annotations + +import importlib + +import numpy +from numpy import pi + +import gridlock +from gridlock import Extent + +from meanas.fdfd import eme, waveguide_2d +from meanas.fdmath import unvec + + +WL = 1310.0 +DX = 40.0 +WIDTH = 400.0 +THF = 161.0 +THP = 77.0 +EPS_SI = 3.51 ** 2 +EPS_OX = 1.453 ** 2 +MODE_NUMBERS = numpy.array([0]) + + +def require_optional(name: str, package_name: str | None = None): + package_name = package_name or name + try: + return importlib.import_module(name) + except ImportError as exc: # pragma: no cover - user environment guard + raise SystemExit( + f"This example requires the optional package '{package_name}'. " + "Install example dependencies with `pip install -e './meanas[examples]'`.", + ) from exc + + +def build_geometry( + *, + dx: float = DX, + width: float = WIDTH, + thf: float = THF, + thp: float = THP, + eps_si: float = EPS_SI, + eps_ox: float = EPS_OX, + ) -> tuple[gridlock.Grid, numpy.ndarray, list[list[numpy.ndarray]], float]: + x0 = (width / 2) % dx + omega = 2 * pi / WL + + grid = gridlock.Grid( + [ + numpy.arange(-800, 800 + dx, dx), + numpy.arange(-400, 400 + dx, dx), + numpy.arange(-2 * dx, 2 * dx + dx, dx), + ], + periodic=True, + ) + epsilon = grid.allocate(eps_ox) + + grid.draw_cuboid( + epsilon, + foreground=eps_si, + x=Extent(center=x0, span=width + 1200), + y=Extent(min=0, max=thf), + z=Extent(min=-1e6, max=0), + ) + grid.draw_cuboid( + epsilon, + foreground=eps_ox, + x=Extent(max=-width / 2, span=300), + y=Extent(min=thp, max=1e6), + z=Extent(min=-1e6, max=0), + ) + grid.draw_cuboid( + epsilon, + foreground=eps_ox, + x=Extent(min=width / 2, span=300), + y=Extent(min=thp, max=1e6), + z=Extent(min=-1e6, max=0), + ) + + grid.draw_cuboid( + epsilon, + foreground=eps_si, + x=Extent(max=-(width / 2 + 600), span=240), + y=Extent(min=0, max=thf), + z=Extent(min=0, max=1e6), + ) + grid.draw_cuboid( + epsilon, + foreground=eps_si, + x=Extent(max=width / 2 + 600, span=240), + y=Extent(min=0, max=thf), + z=Extent(min=0, max=1e6), + ) + + dxes = [grid.dxyz, grid.autoshifted_dxyz()] + dxes_2d = [[d[0], d[1]] for d in dxes] + return grid, epsilon, dxes_2d, omega + + +def solve_cross_section_modes( + epsilon_slice: numpy.ndarray, + *, + omega: float, + dxes_2d: list[list[numpy.ndarray]], + mode_numbers: numpy.ndarray = MODE_NUMBERS, + ) -> tuple[list[tuple[numpy.ndarray, numpy.ndarray]], numpy.ndarray]: + e_xys, wavenumbers = waveguide_2d.solve_modes( + epsilon=epsilon_slice.ravel(), + omega=omega, + dxes=dxes_2d, + mode_numbers=mode_numbers, + ) + eh_fields = [ + waveguide_2d.normalized_fields_e( + e_xy, + wavenumber=wavenumber, + dxes=dxes_2d, + omega=omega, + epsilon=epsilon_slice.ravel(), + ) + for e_xy, wavenumber in zip(e_xys, wavenumbers, strict=True) + ] + return eh_fields, wavenumbers + + +def print_summary(ss: numpy.ndarray, wavenumbers_left: numpy.ndarray, wavenumbers_right: numpy.ndarray, omega: float) -> None: + n_left = len(wavenumbers_left) + left_neff = numpy.real(wavenumbers_left / omega) + right_neff = numpy.real(wavenumbers_right / omega) + + print('left effective indices:', ', '.join(f'{value:.5f}' for value in left_neff[:4])) + print('right effective indices:', ', '.join(f'{value:.5f}' for value in right_neff[:4])) + + reflection = abs(ss[0, 0]) ** 2 + transmission = abs(ss[n_left, 0]) ** 2 + total_output = numpy.sum(abs(ss[:, 0]) ** 2) + print(f'fundamental left-incident reflection |S_00|^2 = {reflection:.6f}') + print(f'fundamental left-to-right transmission |S_{n_left},0|^2 = {transmission:.6f}') + print(f'fundamental left-incident total output power = {total_output:.6f}') + + strongest = numpy.argsort(abs(ss[n_left:, 0]) ** 2)[::-1][:3] + print('dominant transmitted right-side modes for left mode 0:') + for mode_index in strongest: + print(f' mode {mode_index}: |S|^2 = {abs(ss[n_left + mode_index, 0]) ** 2:.6f}') + + +def plot_results( + *, + pyplot, + ss: numpy.ndarray, + left_mode: tuple[numpy.ndarray, numpy.ndarray], + right_mode: tuple[numpy.ndarray, numpy.ndarray], + shape_2d: tuple[int, int], + ) -> None: + fig_s, ax_s = pyplot.subplots() + image = ax_s.imshow(abs(ss) ** 2, origin='lower', cmap='magma') + fig_s.colorbar(image, ax=ax_s) + ax_s.set_title('Interface scattering magnitude |S|^2') + ax_s.set_xlabel('Incident mode index') + ax_s.set_ylabel('Outgoing mode index') + + e_left = unvec(left_mode[0], shape_2d) + e_right = unvec(right_mode[0], shape_2d) + left_intensity = numpy.sum(abs(e_left) ** 2, axis=0).T + right_intensity = numpy.sum(abs(e_right) ** 2, axis=0).T + + fig_modes, axes = pyplot.subplots(1, 2, figsize=(10, 4)) + left_plot = axes[0].imshow(left_intensity, origin='lower', cmap='viridis') + fig_modes.colorbar(left_plot, ax=axes[0]) + axes[0].set_title('Left fundamental mode |E|^2') + right_plot = axes[1].imshow(right_intensity, origin='lower', cmap='viridis') + fig_modes.colorbar(right_plot, ax=axes[1]) + axes[1].set_title('Right fundamental mode |E|^2') + if pyplot.get_backend().lower().endswith('agg'): + pyplot.close(fig_s) + pyplot.close(fig_modes) + else: + pyplot.show() + + +def main() -> None: + pyplot = require_optional('matplotlib.pyplot', package_name='matplotlib') + + grid, epsilon, dxes_2d, omega = build_geometry() + left_slice = epsilon[:, :, :, 1] + right_slice = epsilon[:, :, :, -2] + + left_modes, wavenumbers_left = solve_cross_section_modes(left_slice, omega=omega, dxes_2d=dxes_2d) + right_modes, wavenumbers_right = solve_cross_section_modes(right_slice, omega=omega, dxes_2d=dxes_2d) + + ss = eme.get_s(left_modes, wavenumbers_left, right_modes, wavenumbers_right, dxes=dxes_2d) + + print_summary(ss, wavenumbers_left, wavenumbers_right, omega) + plot_results( + pyplot=pyplot, + ss=ss, + left_mode=left_modes[0], + right_mode=right_modes[0], + shape_2d=grid.shape[:2], + ) + + +if __name__ == '__main__': + main() diff --git a/examples/eme_bend.py b/examples/eme_bend.py new file mode 100644 index 0000000..caff4df --- /dev/null +++ b/examples/eme_bend.py @@ -0,0 +1,310 @@ +""" +Mode-matching / EME example for a straight-to-bent waveguide interface. + +This example demonstrates a cylindrical-waveguide EME workflow: + +1. build a rib-waveguide cross-section, +2. solve straight port modes with `waveguide_2d`, +3. solve bend modes with `waveguide_cyl`, +4. assemble the straight-to-bend interface scattering matrix with + `meanas.fdfd.eme.get_s(...)`, +5. optionally cascade a straight section, bend section, and interface pair into + a compact multiport network using `scikit-rf`. +""" + +from __future__ import annotations + +import importlib + +import numpy +from numpy import pi +from scipy import sparse + +import gridlock +from gridlock import Extent + +from meanas.fdfd import eme, waveguide_2d, waveguide_cyl +from meanas.fdmath import unvec + + +WL = 1310.0 +DX = 40.0 +RADIUS = 6e3 +WIDTH = 400.0 +THF = 161.0 +THP = 77.0 +EPS_SI = 3.51 ** 2 +EPS_OX = 1.453 ** 2 +MODE_NUMBERS = numpy.array([0]) +STRAIGHT_SECTION_LENGTH = 12e3 +BEND_ANGLE = pi / 2 + + +def require_optional(name: str, package_name: str | None = None): + package_name = package_name or name + try: + return importlib.import_module(name) + except ImportError as exc: # pragma: no cover - user environment guard + raise SystemExit( + f"This example requires the optional package '{package_name}'. " + "Install example dependencies with `pip install -e './meanas[examples]'`.", + ) from exc + + +def build_geometry( + *, + dx: float = DX, + width: float = WIDTH, + thf: float = THF, + thp: float = THP, + eps_si: float = EPS_SI, + eps_ox: float = EPS_OX, + ) -> tuple[gridlock.Grid, numpy.ndarray, list[list[numpy.ndarray]], float]: + x0 = (width / 2) % dx + omega = 2 * pi / WL + + grid = gridlock.Grid( + [ + numpy.arange(-800, 800 + dx, dx), + numpy.arange(-400, 400 + dx, dx), + numpy.arange(-2 * dx, 2 * dx + dx, dx), + ], + periodic=True, + ) + epsilon = grid.allocate(eps_ox) + + grid.draw_cuboid( + epsilon, + foreground=eps_si, + x=Extent(center=x0, span=width + 1200), + y=Extent(min=0, max=thf), + z=Extent(min=-1e6, max=0), + ) + grid.draw_cuboid( + epsilon, + foreground=eps_ox, + x=Extent(max=-width / 2, span=300), + y=Extent(min=thp, max=1e6), + z=Extent(min=-1e6, center=0), + ) + grid.draw_cuboid( + epsilon, + foreground=eps_ox, + x=Extent(min=width / 2, span=300), + y=Extent(min=thp, max=1e6), + z=Extent(min=-1e6, center=0), + ) + + dxes = [grid.dxyz, grid.autoshifted_dxyz()] + dxes_2d = [[d[0], d[1]] for d in dxes] + return grid, epsilon, dxes_2d, omega + + +def solve_straight_modes( + epsilon_slice: numpy.ndarray, + *, + omega: float, + dxes_2d: list[list[numpy.ndarray]], + mode_numbers: numpy.ndarray = MODE_NUMBERS, + ) -> tuple[list[tuple[numpy.ndarray, numpy.ndarray]], numpy.ndarray]: + e_xys, wavenumbers = waveguide_2d.solve_modes( + epsilon=epsilon_slice.ravel(), + omega=omega, + dxes=dxes_2d, + mode_numbers=mode_numbers, + ) + eh_fields = [ + waveguide_2d.normalized_fields_e( + e_xy, + wavenumber=wavenumber, + dxes=dxes_2d, + omega=omega, + epsilon=epsilon_slice.ravel(), + ) + for e_xy, wavenumber in zip(e_xys, wavenumbers, strict=True) + ] + return eh_fields, wavenumbers + + +def solve_bend_modes( + epsilon_slice: numpy.ndarray, + *, + omega: float, + dxes_2d: list[list[numpy.ndarray]], + rmin: float, + mode_numbers: numpy.ndarray = MODE_NUMBERS, + ) -> tuple[list[tuple[numpy.ndarray, numpy.ndarray]], numpy.ndarray, numpy.ndarray]: + e_xys, angular_wavenumbers = waveguide_cyl.solve_modes( + epsilon=epsilon_slice.ravel(), + omega=omega, + dxes=dxes_2d, + mode_numbers=mode_numbers, + rmin=rmin, + ) + linear_wavenumbers = waveguide_cyl.linear_wavenumbers( + e_xys=e_xys, + angular_wavenumbers=angular_wavenumbers, + rmin=rmin, + epsilon=epsilon_slice.ravel(), + dxes=dxes_2d, + ) + eh_fields = [ + waveguide_cyl.normalized_fields_e( + e_xy, + angular_wavenumber=angular_wavenumber, + dxes=dxes_2d, + omega=omega, + epsilon=epsilon_slice.ravel(), + rmin=rmin, + ) + for e_xy, angular_wavenumber in zip(e_xys, angular_wavenumbers, strict=True) + ] + return eh_fields, linear_wavenumbers, angular_wavenumbers + + +def build_cascaded_network( + skrf, + *, + interface_s: numpy.ndarray, + straight_wavenumbers: numpy.ndarray, + bend_angular_wavenumbers: numpy.ndarray, + n_modes: int, + ) -> object: + net_sb = skrf.Network(f=[1 / WL], s=interface_s[numpy.newaxis, ...]) + net_bs = net_sb.copy() + net_bs.renumber(numpy.arange(2 * n_modes), numpy.roll(numpy.arange(2 * n_modes), n_modes)) + + straight_phase = sparse.diags_array(numpy.exp(-1j * straight_wavenumbers[:n_modes] * STRAIGHT_SECTION_LENGTH)) + bend_phase = sparse.diags_array(numpy.exp(-1j * bend_angular_wavenumbers[:n_modes] * BEND_ANGLE)) + zero = numpy.zeros((n_modes, n_modes), dtype=complex) + + straight_s = numpy.block([[zero, straight_phase.toarray()], [straight_phase.toarray(), zero]]) + bend_s = numpy.block([[zero, bend_phase.toarray()], [bend_phase.toarray(), zero]]) + net_straight = skrf.Network(f=[1 / WL], s=straight_s[numpy.newaxis, ...]) + net_bend = skrf.Network(f=[1 / WL], s=bend_s[numpy.newaxis, ...]) + + return skrf.network.cascade_list([net_straight, net_sb, net_bend, net_bs, net_straight]) + + +def print_summary( + interface_s: numpy.ndarray, + cascaded_s: numpy.ndarray, + straight_wavenumbers: numpy.ndarray, + bend_linear_wavenumbers: numpy.ndarray, + bend_angular_wavenumbers: numpy.ndarray, + omega: float, + n_modes: int, + ) -> None: + straight_neff = numpy.real(straight_wavenumbers / omega) + bend_neff = numpy.real(bend_linear_wavenumbers / omega) + print('straight effective indices:', ', '.join(f'{value:.5f}' for value in straight_neff[:4])) + print('bend effective indices :', ', '.join(f'{value:.5f}' for value in bend_neff[:4])) + print('bend angular wavenumbers :', ', '.join(f'{value:.5e}' for value in bend_angular_wavenumbers[:4])) + + interface_transmission = abs(interface_s[n_modes, 0]) ** 2 + interface_reflection = abs(interface_s[0, 0]) ** 2 + print(f'interface fundamental transmission |S_{n_modes},0|^2 = {interface_transmission:.6f}') + print(f'interface fundamental reflection |S_00|^2 = {interface_reflection:.6f}') + + total_cascaded_output = numpy.sum(abs(cascaded_s[:, 0]) ** 2) + bend_through = abs(cascaded_s[n_modes, 0]) ** 2 + bend_reflection = abs(cascaded_s[0, 0]) ** 2 + print(f'cascaded bend through power |S_{n_modes},0|^2 = {bend_through:.6f}') + print(f'cascaded bend reflection |S_00|^2 = {bend_reflection:.6f}') + print(f'cascaded left-incident total output power = {total_cascaded_output:.6f}') + + +def plot_results( + *, + pyplot, + interface_s: numpy.ndarray, + cascaded_s: numpy.ndarray, + straight_mode: tuple[numpy.ndarray, numpy.ndarray], + bend_mode: tuple[numpy.ndarray, numpy.ndarray], + shape_2d: tuple[int, int], + ) -> None: + fig_s, axes = pyplot.subplots(1, 2, figsize=(12, 4)) + interface_plot = axes[0].imshow(abs(interface_s) ** 2, origin='lower', cmap='magma') + fig_s.colorbar(interface_plot, ax=axes[0]) + axes[0].set_title('Straight-to-bend interface |S|^2') + axes[0].set_xlabel('Incident mode index') + axes[0].set_ylabel('Outgoing mode index') + + cascaded_plot = axes[1].imshow(abs(cascaded_s) ** 2, origin='lower', cmap='magma') + fig_s.colorbar(cascaded_plot, ax=axes[1]) + axes[1].set_title('Cascaded bend network |S|^2') + axes[1].set_xlabel('Incident mode index') + axes[1].set_ylabel('Outgoing mode index') + + straight_e = unvec(straight_mode[0], shape_2d) + bend_e = unvec(bend_mode[0], shape_2d) + straight_intensity = numpy.sum(abs(straight_e) ** 2, axis=0).T + bend_intensity = numpy.sum(abs(bend_e) ** 2, axis=0).T + + fig_modes, axes_modes = pyplot.subplots(1, 2, figsize=(10, 4)) + straight_plot = axes_modes[0].imshow(straight_intensity, origin='lower', cmap='viridis') + fig_modes.colorbar(straight_plot, ax=axes_modes[0]) + axes_modes[0].set_title('Straight fundamental mode |E|^2') + bend_plot = axes_modes[1].imshow(bend_intensity, origin='lower', cmap='viridis') + fig_modes.colorbar(bend_plot, ax=axes_modes[1]) + axes_modes[1].set_title('Bent fundamental mode |E|^2') + if pyplot.get_backend().lower().endswith('agg'): + pyplot.close(fig_s) + pyplot.close(fig_modes) + else: + pyplot.show() + + +def main() -> None: + pyplot = require_optional('matplotlib.pyplot', package_name='matplotlib') + skrf = require_optional('skrf', package_name='scikit-rf') + + grid, epsilon, dxes_2d, omega = build_geometry() + epsilon_slice = epsilon[:, :, :, 2] + rmin = RADIUS + grid.xyz[0].min() + + straight_modes, straight_wavenumbers = solve_straight_modes(epsilon_slice, omega=omega, dxes_2d=dxes_2d) + bend_modes, bend_linear_wavenumbers, bend_angular_wavenumbers = solve_bend_modes( + epsilon_slice, + omega=omega, + dxes_2d=dxes_2d, + rmin=rmin, + ) + + interface_s = eme.get_s( + straight_modes, + straight_wavenumbers, + bend_modes, + bend_linear_wavenumbers, + dxes=dxes_2d, + ) + cascaded = build_cascaded_network( + skrf, + interface_s=interface_s, + straight_wavenumbers=straight_wavenumbers, + bend_angular_wavenumbers=bend_angular_wavenumbers, + n_modes=len(MODE_NUMBERS), + ) + cascaded_s = cascaded.s[0] + + print_summary( + interface_s, + cascaded_s, + straight_wavenumbers, + bend_linear_wavenumbers, + bend_angular_wavenumbers, + omega, + len(MODE_NUMBERS), + ) + plot_results( + pyplot=pyplot, + interface_s=interface_s, + cascaded_s=cascaded_s, + straight_mode=straight_modes[0], + bend_mode=bend_modes[0], + shape_2d=grid.shape[:2], + ) + + +if __name__ == '__main__': + main() diff --git a/pyproject.toml b/pyproject.toml index 20a694f..013631a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -67,6 +67,7 @@ docs = [ ] examples = [ "matplotlib>=3.10.8", + "scikit-rf>=1.0", ] test = ["pytest", "coverage"]