diff --git a/examples/waveguide.py b/examples/waveguide.py index f243924..7becd59 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) + e_overlap = waveguide_3d.compute_overlap_e(E=wg_results['E'], wavenumber=wg_results['wavenumber'], **wg_args, omega=omega) return J, e_overlap diff --git a/meanas/fdfd/waveguide_2d.py b/meanas/fdfd/waveguide_2d.py index e67160e..aab7944 100644 --- a/meanas/fdfd/waveguide_2d.py +++ b/meanas/fdfd/waveguide_2d.py @@ -373,8 +373,9 @@ 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( - e=e, h=h, dxes=dxes, epsilon=epsilon, prop_phase=prop_phase, + 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, ) return e_norm, h_norm @@ -416,8 +417,9 @@ 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( - e=e, h=h, dxes=dxes, epsilon=epsilon, prop_phase=prop_phase, + 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, ) return e_norm, h_norm @@ -425,8 +427,10 @@ 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_3d.py b/meanas/fdfd/waveguide_3d.py index c77c8d4..01db9b1 100644 --- a/meanas/fdfd/waveguide_3d.py +++ b/meanas/fdfd/waveguide_3d.py @@ -196,6 +196,7 @@ 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/fdfd/waveguide_cyl.py b/meanas/fdfd/waveguide_cyl.py index f2cb5c3..d63accb 100644 --- a/meanas/fdfd/waveguide_cyl.py +++ b/meanas/fdfd/waveguide_cyl.py @@ -529,8 +529,9 @@ 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( - e=e, h=h, dxes=dxes, epsilon=epsilon, prop_phase=prop_phase, + 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, ) return e_norm, h_norm @@ -538,8 +539,11 @@ 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""" diff --git a/meanas/test/test_eme_numerics.py b/meanas/test/test_eme_numerics.py index 3237c1b..2949e4c 100644 --- a/meanas/test/test_eme_numerics.py +++ b/meanas/test/test_eme_numerics.py @@ -218,81 +218,6 @@ 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() @@ -316,176 +241,3 @@ 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 diff --git a/meanas/test/test_waveguide_fdtd_fdfd.py b/meanas/test/test_waveguide_fdtd_fdfd.py index 0488197..ad9a4e6 100644 --- a/meanas/test/test_waveguide_fdtd_fdfd.py +++ b/meanas/test/test_waveguide_fdtd_fdfd.py @@ -380,6 +380,7 @@ 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) @@ -487,6 +488,7 @@ 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, @@ -504,6 +506,7 @@ 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) @@ -618,6 +621,7 @@ 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 d3ec7cd..162c634 100644 --- a/meanas/test/test_waveguide_mode_helpers.py +++ b/meanas/test/test_waveguide_mode_helpers.py @@ -100,6 +100,7 @@ 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) @@ -129,6 +130,7 @@ 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) @@ -156,6 +158,7 @@ def test_waveguide_3d_compute_overlap_e_rejects_empty_overlap_window( axis=0, polarity=polarity, slices=slices, + omega=OMEGA, # type: ignore[call-arg] ) @@ -170,6 +173,7 @@ def test_waveguide_3d_compute_overlap_e_rejects_zero_support_window() -> None: axis=0, polarity=1, slices=slices, + omega=OMEGA, # type: ignore[call-arg] )