[fdfd.eme] enable using leaky modes

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)):
+    test_modes = previous if previous_dual is None else previous_dual
-        overlap = _lorentz_overlap(previous_mode, current_mode, dxes)
+    for index, (previous_mode, current_mode, test_mode) in enumerate(zip(previous, current, test_modes, strict=True)):
-        self_overlap = _lorentz_overlap(previous_mode, previous_mode, dxes)
+        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?
+            A12[ll, rr] = inner_product(eP, hR, dxes=dxes, conj_h=False)    # TODO optimize loop?
-            A21[ll, rr] = inner_product(eR, hL, dxes=dxes, conj_h=False)
+            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)
+    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, **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)
+    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, **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)
 

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