From 35fc67faa3cbc9892abfbbe2a4e9f53bcfa53f77 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Wed, 22 Apr 2026 21:08:12 -0700 Subject: [PATCH 1/3] [compute_overlap_e] remove omega arg (unused) --- examples/waveguide.py | 2 +- meanas/fdfd/waveguide_3d.py | 1 - meanas/test/test_waveguide_fdtd_fdfd.py | 4 ---- meanas/test/test_waveguide_mode_helpers.py | 4 ---- 4 files changed, 1 insertion(+), 10 deletions(-) diff --git a/examples/waveguide.py b/examples/waveguide.py index 7becd59..f243924 100644 --- a/examples/waveguide.py +++ b/examples/waveguide.py @@ -100,7 +100,7 @@ def get_waveguide_mode( # 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) + e_overlap = waveguide_3d.compute_overlap_e(E=wg_results['E'], wavenumber=wg_results['wavenumber'], **wg_args) return J, e_overlap diff --git a/meanas/fdfd/waveguide_3d.py b/meanas/fdfd/waveguide_3d.py index 01db9b1..c77c8d4 100644 --- a/meanas/fdfd/waveguide_3d.py +++ b/meanas/fdfd/waveguide_3d.py @@ -196,7 +196,6 @@ def compute_overlap_e( axis: int, polarity: int, slices: Sequence[slice], - _omega: float, ) -> cfdfield_t: r""" Build an overlap field for projecting another 3D electric field onto a mode. diff --git a/meanas/test/test_waveguide_fdtd_fdfd.py b/meanas/test/test_waveguide_fdtd_fdfd.py index ad9a4e6..0488197 100644 --- a/meanas/test/test_waveguide_fdtd_fdfd.py +++ b/meanas/test/test_waveguide_fdtd_fdfd.py @@ -380,7 +380,6 @@ def _run_straight_waveguide_case(variant: str) -> WaveguideCalibrationResult: axis=0, polarity=1, slices=MONITOR_SLICES, - omega=OMEGA, # type: ignore[call-arg] ) update_e, update_h = fdtd.updates_with_cpml(cpml_params=_build_cpml_params(), dt=DT, dxes=base_dxes, epsilon=epsilon) @@ -488,7 +487,6 @@ def _run_width_step_scattering_case() -> WaveguideScatteringResult: axis=0, polarity=-1, slices=SCATTERING_REFLECT_SLICES, - omega=OMEGA, # type: ignore[call-arg] ) transmitted_mode = waveguide_3d.solve_mode( 0, @@ -506,7 +504,6 @@ def _run_width_step_scattering_case() -> WaveguideScatteringResult: axis=0, polarity=1, slices=SCATTERING_TRANSMIT_SLICES, - omega=OMEGA, # type: ignore[call-arg] ) update_e, update_h = fdtd.updates_with_cpml(cpml_params=_build_cpml_params(), dt=DT, dxes=base_dxes, epsilon=epsilon) @@ -621,7 +618,6 @@ def _run_pulsed_straight_waveguide_case() -> PulsedWaveguideCalibrationResult: axis=0, polarity=1, slices=MONITOR_SLICES, - omega=OMEGA, # type: ignore[call-arg] ) update_e, update_h = fdtd.updates_with_cpml(cpml_params=_build_cpml_params(), dt=DT, dxes=base_dxes, epsilon=epsilon, dtype=complex) diff --git a/meanas/test/test_waveguide_mode_helpers.py b/meanas/test/test_waveguide_mode_helpers.py index 162c634..d3ec7cd 100644 --- a/meanas/test/test_waveguide_mode_helpers.py +++ b/meanas/test/test_waveguide_mode_helpers.py @@ -100,7 +100,6 @@ def test_waveguide_3d_compute_overlap_e_uses_adjacent_window( axis=0, polarity=polarity, slices=slices, - omega=OMEGA, # type: ignore[call-arg] ) nonzero = numpy.argwhere(numpy.abs(overlap) > 0) @@ -130,7 +129,6 @@ def test_waveguide_3d_compute_overlap_e_warns_when_window_is_clipped( axis=0, polarity=polarity, slices=slices, - omega=OMEGA, # type: ignore[call-arg] ) nonzero = numpy.argwhere(numpy.abs(overlap) > 0) @@ -158,7 +156,6 @@ def test_waveguide_3d_compute_overlap_e_rejects_empty_overlap_window( axis=0, polarity=polarity, slices=slices, - omega=OMEGA, # type: ignore[call-arg] ) @@ -173,7 +170,6 @@ def test_waveguide_3d_compute_overlap_e_rejects_zero_support_window() -> None: axis=0, polarity=1, slices=slices, - omega=OMEGA, # type: ignore[call-arg] ) From 061c3f2e907b140b7226ea6aec3231294fbbb7ba Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Wed, 22 Apr 2026 21:09:59 -0700 Subject: [PATCH 2/3] [_normalized_fields] remove unused args --- meanas/fdfd/waveguide_2d.py | 12 ++++-------- meanas/fdfd/waveguide_cyl.py | 8 ++------ 2 files changed, 6 insertions(+), 14 deletions(-) diff --git a/meanas/fdfd/waveguide_2d.py b/meanas/fdfd/waveguide_2d.py index aab7944..e67160e 100644 --- a/meanas/fdfd/waveguide_2d.py +++ b/meanas/fdfd/waveguide_2d.py @@ -373,9 +373,8 @@ def normalized_fields_e( """ e = exy2e(wavenumber=wavenumber, dxes=dxes, epsilon=epsilon) @ e_xy h = exy2h(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) @ e_xy - e_norm, h_norm = _normalized_fields( # type: ignore[call-arg] - e=e, h=h, omega=omega, dxes=dxes, epsilon=epsilon, - mu=mu, prop_phase=prop_phase, + e_norm, h_norm = _normalized_fields( + e=e, h=h, dxes=dxes, epsilon=epsilon, prop_phase=prop_phase, ) return e_norm, h_norm @@ -417,9 +416,8 @@ def normalized_fields_h( """ e = hxy2e(wavenumber=wavenumber, omega=omega, dxes=dxes, epsilon=epsilon, mu=mu) @ h_xy h = hxy2h(wavenumber=wavenumber, dxes=dxes, mu=mu) @ h_xy - e_norm, h_norm = _normalized_fields( # type: ignore[call-arg] - e=e, h=h, omega=omega, dxes=dxes, epsilon=epsilon, - mu=mu, prop_phase=prop_phase, + e_norm, h_norm = _normalized_fields( + e=e, h=h, dxes=dxes, epsilon=epsilon, prop_phase=prop_phase, ) return e_norm, h_norm @@ -427,10 +425,8 @@ def normalized_fields_h( def _normalized_fields( e: vcfdslice, h: vcfdslice, - _omega: complex, dxes: dx_lists2_t, epsilon: vfdslice, - _mu: vfdslice | None = None, prop_phase: float = 0, ) -> tuple[vcfdslice_t, vcfdslice_t]: r""" diff --git a/meanas/fdfd/waveguide_cyl.py b/meanas/fdfd/waveguide_cyl.py index d63accb..f2cb5c3 100644 --- a/meanas/fdfd/waveguide_cyl.py +++ b/meanas/fdfd/waveguide_cyl.py @@ -529,9 +529,8 @@ def normalized_fields_e( """ 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 - e_norm, h_norm = _normalized_fields( # type: ignore[call-arg] - e=e, h=h, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon, - mu=mu, prop_phase=prop_phase, + e_norm, h_norm = _normalized_fields( + e=e, h=h, dxes=dxes, epsilon=epsilon, prop_phase=prop_phase, ) return e_norm, h_norm @@ -539,11 +538,8 @@ def normalized_fields_e( def _normalized_fields( e: vcfdslice, h: vcfdslice, - _omega: complex, dxes: dx_lists2_t, - _rmin: float, # Currently unused, but may want to use cylindrical poynting epsilon: vfdslice, - _mu: vfdslice | None = None, prop_phase: float = 0, ) -> tuple[vcfdslice_t, vcfdslice_t]: r""" From 39291a8314d65a6abeef68a8939931390f2af7c0 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Wed, 22 Apr 2026 21:10:18 -0700 Subject: [PATCH 3/3] [eme] add analytic tests --- meanas/test/test_eme_numerics.py | 248 +++++++++++++++++++++++++++++++ 1 file changed, 248 insertions(+) diff --git a/meanas/test/test_eme_numerics.py b/meanas/test/test_eme_numerics.py index 2949e4c..3237c1b 100644 --- a/meanas/test/test_eme_numerics.py +++ b/meanas/test/test_eme_numerics.py @@ -218,6 +218,81 @@ def _build_bend_mode() -> tuple[tuple[numpy.ndarray, numpy.ndarray], complex]: return (e_field, h_field), linear_wavenumber +def _build_uniform_mode(index: float) -> tuple[tuple[numpy.ndarray, numpy.ndarray], complex]: + area = 25.0 + e_field = numpy.zeros((3, 5, 5), dtype=complex) + h_field = numpy.zeros((3, 5, 5), dtype=complex) + e_field[0] = numpy.sqrt(2.0 / (index * area)) + h_field[1] = numpy.sqrt(2.0 * index / area) + return (vec(e_field), vec(h_field)), complex(index * OMEGA) + + +def _interface_s(n_left: float, n_right: float) -> numpy.ndarray: + left_mode, left_beta = _build_uniform_mode(n_left) + right_mode, right_beta = _build_uniform_mode(n_right) + return eme.get_s([left_mode], [left_beta], [right_mode], [right_beta], dxes=REAL_DXES) + + +def _interface_abcd(n_left: float, n_right: float) -> numpy.ndarray: + left_mode, left_beta = _build_uniform_mode(n_left) + right_mode, right_beta = _build_uniform_mode(n_right) + return eme.get_abcd([left_mode], [left_beta], [right_mode], [right_beta], dxes=REAL_DXES).toarray() + + +def _expected_interface_s(n_left: float, n_right: float) -> numpy.ndarray: + reflection = (n_left - n_right) / (n_left + n_right) + transmission = 2 * numpy.sqrt(n_left * n_right) / (n_left + n_right) + return numpy.array( + [ + [reflection, transmission], + [transmission, -reflection], + ], + dtype=complex, + ) + + +def _propagation_abcd(beta: complex, length: float) -> numpy.ndarray: + phase = numpy.exp(-1j * beta * length) + return numpy.array( + [ + [phase, 0.0], + [0.0, phase ** -1], + ], + dtype=complex, + ) + + +def _abcd_to_s(abcd: numpy.ndarray) -> numpy.ndarray: + aa = abcd[0, 0] + bb = abcd[0, 1] + cc = abcd[1, 0] + dd = abcd[1, 1] + t21 = 1.0 / dd + r21 = bb / dd + r12 = -cc / dd + t12 = aa - bb * cc / dd + return numpy.array( + [ + [r12, t12], + [t21, r21], + ], + dtype=complex, + ) + + +def _expected_bragg_reflector_s(n_low: float, n_high: float, periods: int) -> numpy.ndarray: + ratio = n_high / n_low + reflection = (1 - ratio ** (2 * periods)) / (1 + ratio ** (2 * periods)) + transmission = ((-1) ** periods) * 2 * ratio ** periods / (1 + ratio ** (2 * periods)) + return numpy.array( + [ + [reflection, transmission], + [transmission, -reflection], + ], + dtype=complex, + ) + + def test_get_s_is_near_identity_for_identical_solved_straight_modes() -> None: mode, wavenumber, _epsilon = _build_straight_mode() @@ -241,3 +316,176 @@ def test_get_s_returns_finite_passive_output_for_small_straight_to_bend_fixture( assert ss.shape == (2, 2) assert numpy.isfinite(ss).all() assert numpy.linalg.svd(ss, compute_uv=False).max() <= 1.0 + 1e-10 + + +def test_get_s_matches_analytic_fresnel_interface_for_uniform_one_mode_ports() -> None: + """ + For power-normalized one-mode ports at normal incidence, the interface matrix is + + r12 = (n_left - n_right) / (n_left + n_right) + r21 = -r12 + t12 = t21 = 2 * sqrt(n_left * n_right) / (n_left + n_right) + + so + + S = [[r12, t12], [t21, r21]]. + """ + ss = _interface_s(1.0, 2.0) + expected = _expected_interface_s(1.0, 2.0) + + assert ss.shape == (2, 2) + assert numpy.isfinite(ss).all() + assert_close(ss, expected, atol=1e-6, rtol=1e-6) + assert numpy.linalg.svd(ss, compute_uv=False).max() <= 1.0 + 1e-10 + + +def test_quarter_wave_matching_layer_is_nearly_reflectionless_at_design_frequency() -> None: + """ + A single quarter-wave matching layer with + + n1 = sqrt(n0 * n2), beta1 * L = pi / 2 + + cancels the two interface reflections at the design wavelength, so the + normal-incidence stack should satisfy `r = 0` and `|t| = 1`. + """ + n0 = 1.0 + n1 = numpy.sqrt(2.0) + n2 = 2.0 + interface_01 = _interface_abcd(n0, n1) + interface_12 = _interface_abcd(n1, n2) + _mode_1, beta_1 = _build_uniform_mode(float(n1)) + quarter_wave_length = numpy.pi / (2 * numpy.real(beta_1)) + + stack_abcd = interface_01 @ _propagation_abcd(beta_1, quarter_wave_length) @ interface_12 + ss = _abcd_to_s(stack_abcd) + + assert ss.shape == (2, 2) + assert numpy.isfinite(ss).all() + assert abs(ss[0, 0]) < 1e-12 + assert abs(ss[1, 1]) < 1e-12 + assert abs(abs(ss[0, 1]) - 1.0) < 1e-12 + assert abs(abs(ss[1, 0]) - 1.0) < 1e-12 + assert numpy.linalg.svd(ss, compute_uv=False).max() <= 1.0 + 1e-10 + + +def test_quarter_wave_ar_layer_reduces_reflection_relative_to_abrupt_interface() -> None: + """ + Compare the abrupt interface `n0 -> n2` against the quarter-wave matching-layer + stack `n0 -> sqrt(n0 n2) -> n2` at the same design wavelength. + + For the canonical `n0 = 1`, `n2 = 2` case, the abrupt interface has + + |r_abrupt| = |(n0 - n2) / (n0 + n2)| = 1 / 3, + + while the quarter-wave matching layer should cancel the interface reflections + so that `|r_ar|` is essentially zero and `|t_ar|` is correspondingly larger. + """ + n0 = 1.0 + n2 = 2.0 + abrupt = _interface_s(n0, n2) + + n1 = numpy.sqrt(n0 * n2) + interface_01 = _interface_abcd(n0, n1) + interface_12 = _interface_abcd(n1, n2) + _mode_1, beta_1 = _build_uniform_mode(float(n1)) + quarter_wave_length = numpy.pi / (2 * numpy.real(beta_1)) + ar_stack = _abcd_to_s(interface_01 @ _propagation_abcd(beta_1, quarter_wave_length) @ interface_12) + + abrupt_reflection = abs(abrupt[0, 0]) + abrupt_transmission = abs(abrupt[1, 0]) + ar_reflection = abs(ar_stack[0, 0]) + ar_transmission = abs(ar_stack[1, 0]) + + assert numpy.linalg.svd(abrupt, compute_uv=False).max() <= 1.0 + 1e-10 + assert numpy.linalg.svd(ar_stack, compute_uv=False).max() <= 1.0 + 1e-10 + assert ar_reflection < abrupt_reflection + assert ar_transmission > abrupt_transmission + assert ar_reflection < 1e-12 + assert abs(abrupt_reflection - (1.0 / 3.0)) < 1e-12 + + +def test_half_wave_uniform_slab_restores_unit_transmission_between_matched_media() -> None: + """ + For matched exterior media `n0 = n2`, a half-wave slab with + + beta1 * L = pi + + contributes only a global phase, so the stack returns to `r = 0` and + `|t| = 1` at the design wavelength. + """ + n0 = 1.0 + n1 = 2.0 + interface_01 = _interface_abcd(n0, n1) + interface_10 = _interface_abcd(n1, n0) + _mode_1, beta_1 = _build_uniform_mode(n1) + half_wave_length = numpy.pi / numpy.real(beta_1) + + stack_abcd = interface_01 @ _propagation_abcd(beta_1, half_wave_length) @ interface_10 + ss = _abcd_to_s(stack_abcd) + + assert ss.shape == (2, 2) + assert numpy.isfinite(ss).all() + assert abs(ss[0, 0]) < 1e-12 + assert abs(ss[1, 1]) < 1e-12 + assert abs(abs(ss[0, 1]) - 1.0) < 1e-12 + assert abs(abs(ss[1, 0]) - 1.0) < 1e-12 + assert numpy.linalg.svd(ss, compute_uv=False).max() <= 1.0 + 1e-10 + + +def test_strong_uniform_index_mismatch_behaves_like_near_termination() -> None: + """ + In the large-index-ratio limit, the same Fresnel formulas approach a hard-wall + reflector: + + |r| -> 1, |t| -> 0 as n_right / n_left -> infinity. + """ + ss = _interface_s(1.0, 100.0) + expected = _expected_interface_s(1.0, 100.0) + + assert ss.shape == (2, 2) + assert numpy.isfinite(ss).all() + assert_close(ss, expected, atol=1e-6, rtol=1e-6) + assert abs(ss[0, 0]) > 0.9 + assert abs(ss[1, 0]) < 0.25 + assert numpy.linalg.svd(ss, compute_uv=False).max() <= 1.0 + 1e-10 + + +def test_quarter_wave_bragg_reflector_matches_closed_form_stopband_response() -> None: + """ + For `N` quarter-wave high/low periods at the Bragg wavelength with identical + low-index incident and exit media (`n0 = ns = n_low`), + + M_pair = diag(-(n_low / n_high), -(n_high / n_low)) + M_stack = M_pair ** N + + which yields the closed-form scattering amplitudes + + r = (1 - (n_high / n_low) ** (2N)) / (1 + (n_high / n_low) ** (2N)) + t = 2 * (n_high / n_low) ** N / (1 + (n_high / n_low) ** (2N)). + + The reflector should therefore sit deep in the stopband with `|r|` near 1 and + `|t|` correspondingly small. + """ + n_low = 1.0 + n_high = 2.0 + periods = 5 + interface_lh = _interface_abcd(n_low, n_high) + interface_hl = _interface_abcd(n_high, n_low) + _mode_h, beta_h = _build_uniform_mode(n_high) + _mode_l, beta_l = _build_uniform_mode(n_low) + quarter_wave_high = numpy.pi / (2 * numpy.real(beta_h)) + quarter_wave_low = numpy.pi / (2 * numpy.real(beta_l)) + + stack_abcd = numpy.eye(2, dtype=complex) + for _ in range(periods): + stack_abcd = stack_abcd @ interface_lh @ _propagation_abcd(beta_h, quarter_wave_high) + stack_abcd = stack_abcd @ interface_hl @ _propagation_abcd(beta_l, quarter_wave_low) + ss = _abcd_to_s(stack_abcd) + expected = _expected_bragg_reflector_s(n_low, n_high, periods) + + assert ss.shape == (2, 2) + assert numpy.isfinite(ss).all() + assert_close(ss, expected, atol=1e-12, rtol=1e-12) + assert abs(ss[0, 0]) > 0.99 + assert abs(ss[1, 0]) < 0.1 + assert numpy.linalg.svd(ss, compute_uv=False).max() <= 1.0 + 1e-10