From 444ae49a7409a31f92a1fc38e52b68a7bc1eabfd Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Mon, 4 May 2026 19:03:30 -0700 Subject: [PATCH] [fdfd.eme] enable using leaky modes --- meanas/fdfd/eme.py | 128 ++++++++++++++++++++++------ meanas/test/test_eme_numerics.py | 140 +++++++++++++++++++++++++++++++ 2 files changed, 244 insertions(+), 24 deletions(-) diff --git a/meanas/fdfd/eme.py b/meanas/fdfd/eme.py index 6844476..a9de4dd 100644 --- a/meanas/fdfd/eme.py +++ b/meanas/fdfd/eme.py @@ -49,6 +49,11 @@ class ModeSection: modes are synthesized as `(E, -H)`. backward_wavenumbers: Optional propagation constants for `backward_modes`. If omitted, they are synthesized as `-wavenumbers`. + dual_modes: Optional forward dual / adjoint projection modes. If + omitted, `modes` are used as their own projection basis. + dual_backward_modes: Optional backward dual / adjoint projection modes. + If omitted, they are synthesized from `dual_modes` when available, + otherwise from `backward_modes`. """ z: float @@ -56,6 +61,8 @@ class ModeSection: wavenumbers: wavenumber_seq backward_modes: mode_seq | None = None backward_wavenumbers: wavenumber_seq | None = None + dual_modes: mode_seq | None = None + dual_backward_modes: mode_seq | None = None def _validate_port_modes( @@ -89,6 +96,21 @@ def _validate_port_modes( return e_shape, h_shape +def _validate_dual_modes( + name: str, + dual_ehs: mode_seq | None, + reference_shape: tuple[int, ...], + wavenumbers: wavenumber_seq, + ) -> mode_seq | None: + if dual_ehs is None: + return None + + dual_e_shape, dual_h_shape = _validate_port_modes(name, dual_ehs, wavenumbers) + if dual_e_shape != reference_shape or dual_h_shape != reference_shape: + raise ValueError(f'{name} modal fields must share the same E/H shapes as the corresponding modes') + return dual_ehs + + def _as_wavenumber_array( name: str, wavenumbers: wavenumber_seq, @@ -140,15 +162,17 @@ def _phase_align_modes( previous: Sequence[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]], current: Sequence[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]], dxes: dx_lists2_t, + previous_dual: Sequence[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]] | None = None, ) -> list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]]: aligned = [] - for index, (previous_mode, current_mode) in enumerate(zip(previous, current, strict=True)): - overlap = _lorentz_overlap(previous_mode, current_mode, dxes) - self_overlap = _lorentz_overlap(previous_mode, previous_mode, dxes) + test_modes = previous if previous_dual is None else previous_dual + for index, (previous_mode, current_mode, test_mode) in enumerate(zip(previous, current, test_modes, strict=True)): + overlap = _lorentz_overlap(test_mode, current_mode, dxes) + self_overlap = _lorentz_overlap(test_mode, previous_mode, dxes) if overlap == 0: raise ValueError(f'cannot phase-track mode {index}: adjacent section overlap is zero') if self_overlap == 0: - raise ValueError(f'cannot phase-track mode {index}: mode self-overlap is zero') + raise ValueError(f'cannot phase-track mode {index}: mode dual-overlap is zero') phase = (overlap / abs(overlap)) / (self_overlap / abs(self_overlap)) e_field, h_field = current_mode aligned.append((e_field / phase, h_field / phase)) @@ -166,6 +190,8 @@ def _section_branches( NDArray[numpy.complex128], list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]], NDArray[numpy.complex128], + list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]], + list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]], tuple[int, ...], ]: z_coord = float(section.z) @@ -203,7 +229,37 @@ def _section_branches( if len(backward_wavenumbers) != len(wavenumbers): raise ValueError('forward and backward mode counts must match') - return z_coord, forward_modes, wavenumbers, backward_modes, backward_wavenumbers, shape + if section.dual_modes is None: + dual_modes = forward_modes + else: + dual_shape, _dual_h_shape = _validate_port_modes( + f'sections[{index}].dual_modes', + section.dual_modes, + section.wavenumbers, + ) + if dual_shape != shape: + raise ValueError('modal fields and dual modal fields must share the same E/H shapes') + dual_modes = _as_mode_arrays(section.dual_modes) + + if section.dual_backward_modes is None: + if section.dual_modes is None and section.backward_modes is not None: + dual_backward_modes = backward_modes + else: + dual_backward_modes = [(e_field.copy(), -h_field.copy()) for e_field, h_field in dual_modes] + else: + dual_backward_shape, _dual_backward_h_shape = _validate_port_modes( + f'sections[{index}].dual_backward_modes', + section.dual_backward_modes, + section.backward_wavenumbers if section.backward_wavenumbers is not None else backward_wavenumbers, + ) + if dual_backward_shape != shape: + raise ValueError('backward modal fields and dual backward modal fields must share the same E/H shapes') + dual_backward_modes = _as_mode_arrays(section.dual_backward_modes) + + if len(dual_modes) != len(forward_modes) or len(dual_backward_modes) != len(backward_modes): + raise ValueError('dual mode counts must match the corresponding modal basis counts') + + return z_coord, forward_modes, wavenumbers, backward_modes, backward_wavenumbers, dual_modes, dual_backward_modes, shape def _validate_taper_sections( @@ -212,6 +268,7 @@ def _validate_taper_sections( ) -> tuple[ NDArray[numpy.float64], list[list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]]], + list[list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]]], list[NDArray[numpy.complex128]], int, ]: @@ -220,12 +277,14 @@ def _validate_taper_sections( z_coords: list[float] = [] branch_modes: list[list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]]] = [] + branch_dual_modes: list[list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]]] = [] branch_wavenumbers: list[NDArray[numpy.complex128]] = [] + explicit_duals: list[bool] = [] expected_count: int | None = None expected_shape: tuple[int, ...] | None = None for index, section in enumerate(sections): - z_coord, forward_modes, forward_wavenumbers, backward_modes, backward_wavenumbers, shape = _section_branches( + z_coord, forward_modes, forward_wavenumbers, backward_modes, backward_wavenumbers, dual_modes, dual_backward_modes, shape = _section_branches( section, index, expected_count, @@ -236,21 +295,26 @@ def _validate_taper_sections( expected_shape = shape z_coords.append(z_coord) branch_modes.append([*forward_modes, *backward_modes]) + branch_dual_modes.append([*dual_modes, *dual_backward_modes]) branch_wavenumbers.append(numpy.concatenate((forward_wavenumbers, backward_wavenumbers))) + explicit_duals.append(section.dual_modes is not None or section.dual_backward_modes is not None) z_array = numpy.asarray(z_coords, dtype=float) if not (numpy.diff(z_array) > 0).all(): raise ValueError('taper section z coordinates must be strictly increasing') for index in range(1, len(branch_modes)): - branch_modes[index] = _phase_align_modes(branch_modes[index - 1], branch_modes[index], dxes) + branch_modes[index] = _phase_align_modes(branch_modes[index - 1], branch_modes[index], dxes, branch_dual_modes[index - 1]) + if not explicit_duals[index]: + branch_dual_modes[index] = branch_modes[index] assert expected_count is not None - return z_array, branch_modes, branch_wavenumbers, expected_count + return z_array, branch_modes, branch_dual_modes, branch_wavenumbers, expected_count def _taper_interval_generator( left_modes: Sequence[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]], + left_dual_modes: Sequence[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]], right_modes: Sequence[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]], left_wavenumbers: NDArray[numpy.complex128], right_wavenumbers: NDArray[numpy.complex128], @@ -261,7 +325,7 @@ def _taper_interval_generator( gram = numpy.zeros((mode_count, mode_count), dtype=complex) derivative_overlap = numpy.zeros((mode_count, mode_count), dtype=complex) - for row, left_row_mode in enumerate(left_modes): + for row, left_row_mode in enumerate(left_dual_modes): for col, left_col_mode in enumerate(left_modes): gram[row, col] = _lorentz_overlap(left_row_mode, left_col_mode, dxes) for col, (left_col_mode, right_col_mode) in enumerate(zip(left_modes, right_modes, strict=True)): @@ -294,11 +358,12 @@ def _abcd_to_s( def get_tr( - ehLs: Sequence[Sequence[vcfdfield2]], + ehLs: mode_seq, wavenumbers_L: wavenumber_seq, - ehRs: Sequence[Sequence[vcfdfield2]], + ehRs: mode_seq, wavenumbers_R: wavenumber_seq, dxes: dx_lists2_t, + dual_ehLs: mode_seq | None = None, ) -> tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]: """ Compute left-incident transmission and reflection matrices. @@ -309,6 +374,8 @@ def get_tr( ehRs: Right-port modes as `(E, H)` field pairs. wavenumbers_R: Propagation constants for `ehRs`. dxes: Two-dimensional Yee-cell edge lengths for the shared port plane. + dual_ehLs: Optional left-port dual / adjoint projection modes. If + omitted, `ehLs` are used as their own projection basis. Returns: `(T12, R12)` where columns index left-incident modes and rows index @@ -322,6 +389,8 @@ def get_tr( right_e_shape, right_h_shape = _validate_port_modes('ehRs', ehRs, wavenumbers_R) if left_e_shape != right_e_shape or left_h_shape != right_h_shape: raise ValueError('left and right modal fields must share the same E/H shapes') + dual_projection_ehLs = _validate_dual_modes('dual_ehLs', dual_ehLs, left_e_shape, wavenumbers_L) + projection_ehLs = ehLs if dual_projection_ehLs is None else dual_projection_ehLs nL = len(wavenumbers_L) nR = len(wavenumbers_R) @@ -330,11 +399,12 @@ def get_tr( B11 = numpy.zeros((nL,), dtype=complex) for ll in range(nL): eL, hL = ehLs[ll] - B11[ll] = inner_product(eL, hL, dxes=dxes, conj_h=False) + eP, hP = projection_ehLs[ll] + B11[ll] = inner_product(eL, hP, dxes=dxes, conj_h=False) for rr in range(nR): eR, hR = ehRs[rr] - A12[ll, rr] = inner_product(eL, hR, dxes=dxes, conj_h=False) # TODO optimize loop? - A21[ll, rr] = inner_product(eR, hL, dxes=dxes, conj_h=False) + A12[ll, rr] = inner_product(eP, hR, dxes=dxes, conj_h=False) # TODO optimize loop? + A21[ll, rr] = inner_product(eR, hP, dxes=dxes, conj_h=False) # tt0 = 2 * numpy.linalg.pinv(A21 + numpy.conj(A12)) tt0, _resid, _rank, _sing = numpy.linalg.lstsq(A21 + A12, numpy.diag(2 * B11), rcond=None) @@ -351,10 +421,12 @@ def get_tr( def get_abcd( - ehLs: Sequence[Sequence[vcfdfield2]], + ehLs: mode_seq, wavenumbers_L: wavenumber_seq, - ehRs: Sequence[Sequence[vcfdfield2]], + ehRs: mode_seq, wavenumbers_R: wavenumber_seq, + dual_ehLs: mode_seq | None = None, + dual_ehRs: mode_seq | None = None, **kwargs, ) -> sparse.sparray: """ @@ -367,8 +439,8 @@ def get_abcd( convention. """ - t12, r12 = get_tr(ehLs, wavenumbers_L, ehRs, wavenumbers_R, **kwargs) - t21, r21 = get_tr(ehRs, wavenumbers_R, ehLs, wavenumbers_L, **kwargs) + t12, r12 = get_tr(ehLs, wavenumbers_L, ehRs, wavenumbers_R, dual_ehLs=dual_ehLs, **kwargs) + t21, r21 = get_tr(ehRs, wavenumbers_R, ehLs, wavenumbers_L, dual_ehLs=dual_ehRs, **kwargs) t21i = numpy.linalg.pinv(t21) A = t12 - r21 @ t21i @ r12 B = r21 @ t21i @@ -384,12 +456,14 @@ def get_abcd( def get_s( - ehLs: Sequence[Sequence[vcfdfield2]], + ehLs: mode_seq, wavenumbers_L: wavenumber_seq, - ehRs: Sequence[Sequence[vcfdfield2]], + ehRs: mode_seq, wavenumbers_R: wavenumber_seq, force_nogain: bool = False, force_reciprocal: bool = False, + dual_ehLs: mode_seq | None = None, + dual_ehRs: mode_seq | None = None, **kwargs, ) -> NDArray[numpy.complex128]: """ @@ -404,9 +478,11 @@ def get_s( scattering matrix to at most one. force_reciprocal: If `True`, symmetrize the assembled matrix as `0.5 * (S + S.T)`. + dual_ehLs: Optional left-port dual / adjoint projection modes. + dual_ehRs: Optional right-port dual / adjoint projection modes. """ - t12, r12 = get_tr(ehLs, wavenumbers_L, ehRs, wavenumbers_R, **kwargs) - t21, r21 = get_tr(ehRs, wavenumbers_R, ehLs, wavenumbers_L, **kwargs) + t12, r12 = get_tr(ehLs, wavenumbers_L, ehRs, wavenumbers_R, dual_ehLs=dual_ehLs, **kwargs) + t21, r21 = get_tr(ehRs, wavenumbers_R, ehLs, wavenumbers_L, dual_ehLs=dual_ehRs, **kwargs) ss = numpy.block([[r12, t12], [t21, r21]]) @@ -437,6 +513,9 @@ def get_taper_abcd( coordinates. Adjacent modal phases are tracked with the same non-conjugated Lorentz/Poynting overlap used by the abrupt-interface helpers, then each interval is propagated with a finite-difference local-mode CMT generator. + If a `ModeSection` supplies dual / adjoint modes, those modes are used for + the CMT projection. This supports leaky or radiative mode sets whose natural + projection basis is biorthogonal rather than self-projected. Args: sections: Local modal samples ordered by increasing `z`. @@ -459,13 +538,14 @@ def get_taper_abcd( if max_step is not None and max_step <= 0: raise ValueError('max_step must be positive') - z_coords, branch_modes, branch_wavenumbers, n_modes = _validate_taper_sections(sections, dxes) + z_coords, branch_modes, branch_dual_modes, branch_wavenumbers, n_modes = _validate_taper_sections(sections, dxes) branch_count = 2 * n_modes transfer = numpy.eye(branch_count, dtype=complex) for index, dz in enumerate(numpy.diff(z_coords)): generator = _taper_interval_generator( branch_modes[index], + branch_dual_modes[index], branch_modes[index + 1], branch_wavenumbers[index], branch_wavenumbers[index + 1], @@ -494,7 +574,7 @@ def get_taper_s( The returned matrix uses the same ordering as `get_s(...)`: `[[R12, T12], [T21, R21]]`. """ - _z_coords, _branch_modes, _branch_wavenumbers, n_modes = _validate_taper_sections(sections, dxes) + _z_coords, _branch_modes, _branch_dual_modes, _branch_wavenumbers, n_modes = _validate_taper_sections(sections, dxes) abcd = get_taper_abcd(sections, dxes, **kwargs).toarray() ss = _abcd_to_s(abcd, n_modes) diff --git a/meanas/test/test_eme_numerics.py b/meanas/test/test_eme_numerics.py index ed4708a..0d28692 100644 --- a/meanas/test/test_eme_numerics.py +++ b/meanas/test/test_eme_numerics.py @@ -77,6 +77,27 @@ def test_get_tr_returns_finite_bounded_transfer_matrices() -> None: assert (singular_values <= 1.0 + 1e-12).all() +def test_get_tr_accepts_scaled_dual_projection_modes() -> None: + left_modes, right_modes = _mode_sets() + dual_left_modes = [ + (mode[0] * (0.5 + 0.25j), mode[1] * (0.5 + 0.25j)) + for mode in left_modes + ] + + plain_t, plain_r = eme.get_tr(left_modes, WAVENUMBERS_L, right_modes, WAVENUMBERS_R, dxes=DXES) + dual_t, dual_r = eme.get_tr( + left_modes, + WAVENUMBERS_L, + right_modes, + WAVENUMBERS_R, + dxes=DXES, + dual_ehLs=dual_left_modes, + ) + + assert_close(dual_t, plain_t) + assert_close(dual_r, plain_r) + + def test_get_abcd_matches_explicit_block_formula() -> None: left_modes, right_modes = _mode_sets() t12, r12 = eme.get_tr(left_modes, WAVENUMBERS_L, right_modes, WAVENUMBERS_R, dxes=DXES) @@ -166,6 +187,20 @@ def test_get_tr_rejects_incompatible_field_shapes() -> None: eme.get_tr(left_modes, [1.0], right_modes, [1.0], dxes=DXES) +def test_get_tr_rejects_dual_mode_length_mismatches() -> None: + left_modes, right_modes = _mode_sets() + + with pytest.raises(ValueError, match='same length'): + eme.get_tr( + left_modes, + WAVENUMBERS_L, + right_modes, + WAVENUMBERS_R, + dxes=DXES, + dual_ehLs=left_modes[:1], + ) + + def _build_real_epsilon() -> numpy.ndarray: epsilon = numpy.ones((3, 5, 5), dtype=float) epsilon[:, 2, 1] = 2.0 @@ -277,6 +312,70 @@ def test_get_taper_abcd_is_invariant_to_adjacent_modal_phase() -> None: assert_close(shifted, base, atol=1e-12, rtol=1e-12) +def test_get_taper_abcd_is_invariant_to_modal_phase_across_multiple_sections() -> None: + mode, beta = _build_uniform_mode(1.5) + mid_length = 5.0 + length = 11.0 + base_sections = [ + eme.ModeSection(0.0, [mode], [beta]), + eme.ModeSection(mid_length, [mode], [beta]), + eme.ModeSection(length, [mode], [beta]), + ] + shifted_sections = [ + eme.ModeSection(0.0, [mode], [beta]), + eme.ModeSection(mid_length, [(mode[0] * numpy.exp(0.41j), mode[1] * numpy.exp(0.41j))], [beta]), + eme.ModeSection(length, [(mode[0] * numpy.exp(-0.28j), mode[1] * numpy.exp(-0.28j))], [beta]), + ] + + base = eme.get_taper_abcd(base_sections, dxes=REAL_DXES).toarray() + shifted = eme.get_taper_abcd(shifted_sections, dxes=REAL_DXES).toarray() + + assert_close(shifted, base, atol=1e-12, rtol=1e-12) + + +def test_get_taper_accepts_complex_leaky_wavenumber() -> None: + mode, beta = _build_uniform_mode(1.5) + leaky_beta = beta - 0.02j + length = 3.0 + + abcd = eme.get_taper_abcd( + [ + eme.ModeSection(0.0, [mode], [leaky_beta]), + eme.ModeSection(length, [mode], [leaky_beta]), + ], + dxes=REAL_DXES, + ).toarray() + + assert_close(abcd, _propagation_abcd(leaky_beta, length), atol=1e-12, rtol=1e-12) + + +def test_get_taper_uses_supplied_dual_modes_for_phase_tracking() -> None: + mode, beta = _build_uniform_mode(1.5) + primal_scale = numpy.exp(0.42j) + dual_scale = 0.31 * numpy.exp(-0.77j) + dual_mode = (mode[0] * dual_scale, mode[1] * dual_scale) + shifted_mode = (mode[0] * primal_scale, mode[1] * primal_scale) + shifted_dual_mode = (dual_mode[0] * 2.3j, dual_mode[1] * 2.3j) + length = 11.0 + + base = eme.get_taper_abcd( + [ + eme.ModeSection(0.0, [mode], [beta], dual_modes=[dual_mode]), + eme.ModeSection(length, [mode], [beta], dual_modes=[dual_mode]), + ], + dxes=REAL_DXES, + ).toarray() + shifted = eme.get_taper_abcd( + [ + eme.ModeSection(0.0, [mode], [beta], dual_modes=[dual_mode]), + eme.ModeSection(length, [shifted_mode], [beta], dual_modes=[shifted_dual_mode]), + ], + dxes=REAL_DXES, + ).toarray() + + assert_close(shifted, base, atol=1e-12, rtol=1e-12) + + def test_get_taper_rejects_nonmonotonic_sections() -> None: mode, beta = _build_uniform_mode(1.5) @@ -303,6 +402,19 @@ def test_get_taper_rejects_mode_count_changes() -> None: ) +def test_get_taper_rejects_dual_mode_count_changes() -> None: + mode, beta = _build_uniform_mode(1.5) + + with pytest.raises(ValueError, match='same length'): + eme.get_taper_abcd( + [ + eme.ModeSection(0.0, [mode], [beta], dual_modes=[mode]), + eme.ModeSection(1.0, [mode], [beta], dual_modes=[mode, mode]), + ], + dxes=REAL_DXES, + ) + + 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) @@ -415,6 +527,34 @@ def test_get_s_matches_analytic_fresnel_interface_for_uniform_one_mode_ports() - assert numpy.linalg.svd(ss, compute_uv=False).max() <= 1.0 + 1e-10 +def test_get_s_with_dual_modes_matches_analytic_fresnel_interface() -> None: + left_mode, left_beta = _build_uniform_mode(1.0) + right_mode, right_beta = _build_uniform_mode(2.0) + left_dual = (left_mode[0] * (0.25 + 0.5j), left_mode[1] * (0.25 + 0.5j)) + right_dual = (right_mode[0] * (-0.75 + 0.125j), right_mode[1] * (-0.75 + 0.125j)) + + ss = eme.get_s( + [left_mode], + [left_beta], + [right_mode], + [right_beta], + dxes=REAL_DXES, + dual_ehLs=[left_dual], + dual_ehRs=[right_dual], + ) + expected = _expected_interface_s(1.0, 2.0) + + assert_close(ss, expected, atol=1e-6, rtol=1e-6) + + +def test_get_s_accepts_complex_leaky_wavenumbers_for_abrupt_interface() -> None: + mode, beta = _build_uniform_mode(1.5) + + ss = eme.get_s([mode], [beta - 0.02j], [mode], [beta - 0.03j], dxes=REAL_DXES) + + assert_close(ss, numpy.array([[0.0, 1.0], [1.0, 0.0]], dtype=complex), atol=1e-12, rtol=1e-12) + + def test_quarter_wave_matching_layer_is_nearly_reflectionless_at_design_frequency() -> None: """ A single quarter-wave matching layer with