[eme] add analytic tests

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