From 39291a8314d65a6abeef68a8939931390f2af7c0 Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Wed, 22 Apr 2026 21:10:18 -0700 Subject: [PATCH] [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