diff --git a/README.md b/README.md
index 73d48b5..d179c9d 100644
--- a/README.md
+++ b/README.md
@@ -172,6 +172,8 @@ The tracked examples under `examples/` are the intended entry points for users:
 - `examples/eme_bend.py`: straight-to-bent waveguide mode matching with
   cylindrical bend modes, interface scattering, and a cascaded bend-network
   example built with `scikit-rf`.
+- `examples/eme_taper_cmt.py`: sampled cross-section local-mode CMT for a
+  continuously varying rib-waveguide taper.
 - `examples/fdfd.py`: direct frequency-domain waveguide excitation and overlap /
   Poynting analysis without a time-domain run.
 
diff --git a/docs/index.md b/docs/index.md
index 5475af8..1b5ca6b 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -28,6 +28,8 @@ Relevant starting examples:
   scattering between two nearby waveguide cross-sections
 - `examples/eme_bend.py` for straight-to-bent mode matching with cylindrical
   bend modes and a cascaded bend-network example
+- `examples/eme_taper_cmt.py` for local-mode CMT through sampled continuously
+  varying taper cross-sections
 - `examples/fdfd.py` for direct frequency-domain waveguide excitation
 
 For solver equivalence, prefer the phasor-based examples first. They compare
diff --git a/examples/eme_taper_cmt.py b/examples/eme_taper_cmt.py
new file mode 100644
index 0000000..c5a8f1c
--- /dev/null
+++ b/examples/eme_taper_cmt.py
@@ -0,0 +1,134 @@
+"""
+Local-mode CMT example for a continuously varying rib-waveguide taper.
+
+This example keeps geometry construction outside `meanas.fdfd.eme`: it samples a
+width taper at several cross-sections, solves and normalizes each local mode with
+`waveguide_2d`, then asks `eme.get_taper_s(...)` for the bidirectional taper
+scattering matrix.
+"""
+
+from __future__ import annotations
+
+import numpy
+from numpy import pi
+
+from meanas.fdmath import vec
+from meanas.fdfd import eme, waveguide_2d
+
+
+WL = 1310.0
+DX = 80.0
+TAPER_LENGTH = 4e3
+WIDTH_LEFT = 280.0
+WIDTH_RIGHT = 520.0
+THF = 160.0
+THP = 80.0
+EPS_SI = 3.51 ** 2
+EPS_OX = 1.453 ** 2
+MODE_NUMBERS = numpy.array([0])
+N_SECTIONS = 7
+
+
+def build_dxes(shape: tuple[int, int], dx: float = DX) -> list[list[numpy.ndarray]]:
+    nx, ny = shape
+    return [
+        [numpy.full(nx, dx), numpy.full(ny, dx)],
+        [numpy.full(nx, dx), numpy.full(ny, dx)],
+    ]
+
+
+def build_cross_section(
+        *,
+        width: float,
+        x: numpy.ndarray,
+        y: numpy.ndarray,
+        eps_si: float = EPS_SI,
+        eps_ox: float = EPS_OX,
+        thf: float = THF,
+        thp: float = THP,
+        ) -> numpy.ndarray:
+    epsilon = numpy.full((3, x.size, y.size), eps_ox, dtype=float)
+    xx = x[:, None]
+    yy = y[None, :]
+    slab = (yy >= 0) & (yy <= thp)
+    rib = (abs(xx) <= width / 2) & (yy >= 0) & (yy <= thf)
+    epsilon[:, slab.repeat(x.size, axis=0)] = eps_si
+    epsilon[:, rib] = eps_si
+    return epsilon
+
+
+def solve_cross_section_modes(
+        epsilon: numpy.ndarray,
+        *,
+        omega: float,
+        dxes_2d: list[list[numpy.ndarray]],
+        mode_numbers: numpy.ndarray = MODE_NUMBERS,
+        ) -> tuple[list[tuple[numpy.ndarray, numpy.ndarray]], numpy.ndarray]:
+    epsilon_vec = vec(epsilon)
+    e_xys, wavenumbers = waveguide_2d.solve_modes(
+        epsilon=epsilon_vec,
+        omega=omega,
+        dxes=dxes_2d,
+        mode_numbers=mode_numbers,
+    )
+    eh_fields = [
+        waveguide_2d.normalized_fields_e(
+            e_xy,
+            wavenumber=wavenumber,
+            dxes=dxes_2d,
+            omega=omega,
+            epsilon=epsilon_vec,
+        )
+        for e_xy, wavenumber in zip(e_xys, wavenumbers, strict=True)
+    ]
+    return eh_fields, wavenumbers
+
+
+def solve_taper_sections() -> tuple[list[eme.ModeSection], list[float], float, list[list[numpy.ndarray]]]:
+    omega = 2 * pi / WL
+    x = numpy.arange(-480, 480 + DX, DX)
+    y = numpy.arange(-240, 400 + DX, DX)
+    dxes_2d = build_dxes((x.size, y.size))
+    z_samples = numpy.linspace(0, TAPER_LENGTH, N_SECTIONS)
+    widths = numpy.linspace(WIDTH_LEFT, WIDTH_RIGHT, N_SECTIONS)
+
+    sections = []
+    neffs = []
+    for z_coord, width in zip(z_samples, widths, strict=True):
+        epsilon = build_cross_section(width=float(width), x=x, y=y)
+        modes, wavenumbers = solve_cross_section_modes(epsilon, omega=omega, dxes_2d=dxes_2d)
+        sections.append(eme.ModeSection(float(z_coord), modes, wavenumbers))
+        neffs.append(float(numpy.real(wavenumbers[0] / omega)))
+
+    return sections, neffs, omega, dxes_2d
+
+
+def print_summary(
+        taper_s: numpy.ndarray,
+        abrupt_s: numpy.ndarray,
+        neffs: list[float],
+        ) -> None:
+    n_modes = len(MODE_NUMBERS)
+    print('sampled taper effective indices:', ', '.join(f'{value:.5f}' for value in neffs))
+    print(f'abrupt endpoint reflection       |S_00|^2          = {abs(abrupt_s[0, 0]) ** 2:.6f}')
+    print(f'abrupt endpoint transmission     |S_{n_modes},0|^2 = {abs(abrupt_s[n_modes, 0]) ** 2:.6f}')
+    print(f'taper CMT reflection             |S_00|^2          = {abs(taper_s[0, 0]) ** 2:.6f}')
+    print(f'taper CMT transmission           |S_{n_modes},0|^2 = {abs(taper_s[n_modes, 0]) ** 2:.6f}')
+    print(f'taper CMT total output power                       = {numpy.sum(abs(taper_s[:, 0]) ** 2):.6f}')
+
+
+def main() -> None:
+    sections, neffs, _omega, dxes_2d = solve_taper_sections()
+    taper_s = eme.get_taper_s(sections, dxes=dxes_2d)
+    abrupt_s = eme.get_s(
+        sections[0].modes,
+        sections[0].wavenumbers,
+        sections[-1].modes,
+        sections[-1].wavenumbers,
+        dxes=dxes_2d,
+    )
+    print_summary(taper_s, abrupt_s, neffs)
+
+
+if __name__ == '__main__':
+    main()
diff --git a/meanas/fdfd/eme.py b/meanas/fdfd/eme.py
index af745e8..a9de4dd 100644
--- a/meanas/fdfd/eme.py
+++ b/meanas/fdfd/eme.py
@@ -13,6 +13,9 @@ The returned matrices follow the usual port ordering:
   directional `T/R` solves.
 - `get_s(...)` returns the full block scattering matrix
   `[[R12, T12], [T21, R21]]`.
+- `get_taper_abcd(...)` and `get_taper_s(...)` build the same transfer /
+  scattering objects for sampled continuously varying sections using local-mode
+  CMT.
 
 This module is intentionally a thin library layer rather than an integrated
 simulation suite. It provides the overlap algebra that downstream users can
@@ -20,19 +23,51 @@ compose into larger workflows.
 """
 
 from collections.abc import Sequence
+import dataclasses
 import numpy
 from numpy.typing import NDArray
+from scipy import linalg
 from scipy import sparse
 
 from ..fdmath import dx_lists2_t, vcfdfield2
 from .waveguide_2d import inner_product
 
 type wavenumber_seq = Sequence[complex] | NDArray[numpy.complexfloating] | NDArray[numpy.floating]
+type mode_seq = Sequence[Sequence[vcfdfield2]]
+
+
+@dataclasses.dataclass(frozen=True)
+class ModeSection:
+    """
+    Local modal basis at one longitudinal sample of a continuously varying guide.
+
+    Args:
+        z: Longitudinal coordinate of this section.
+        modes: Forward modes as `(E, H)` field pairs.
+        wavenumbers: Forward propagation constants for `modes`.
+        backward_modes: Optional explicit backward modes. If omitted, backward
+            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
+    modes: mode_seq
+    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(
         name: str,
-        ehs: Sequence[Sequence[vcfdfield2]],
+        ehs: mode_seq,
         wavenumbers: wavenumber_seq,
         ) -> tuple[tuple[int, ...], tuple[int, ...]]:
     if len(ehs) != len(wavenumbers):
@@ -61,12 +96,274 @@ 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,
+        ) -> NDArray[numpy.complex128]:
+    array = numpy.asarray(wavenumbers, dtype=complex)
+    if array.ndim != 1:
+        raise ValueError(f'{name} must be a one-dimensional sequence')
+    if not numpy.isfinite(array).all():
+        raise ValueError(f'{name} must contain only finite values')
+    return array
+
+
+def _as_mode_arrays(
+        ehs: mode_seq,
+        ) -> list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]]:
+    return [
+        (numpy.asarray(e_field, dtype=complex), numpy.asarray(h_field, dtype=complex))
+        for e_field, h_field in ehs
+        ]
+
+
+def _lorentz_overlap(
+        mode_a: tuple[vcfdfield2, vcfdfield2],
+        mode_b: tuple[vcfdfield2, vcfdfield2],
+        dxes: dx_lists2_t,
+        ) -> complex:
+    e_a, h_a = mode_a
+    e_b, h_b = mode_b
+    return 0.5 * (
+        inner_product(e_a, h_b, dxes=dxes, conj_h=False)
+        + inner_product(e_b, h_a, dxes=dxes, conj_h=False)
+        )
+
+
+def _lorentz_derivative_overlap(
+        mode_a: tuple[vcfdfield2, vcfdfield2],
+        derivative_b: tuple[vcfdfield2, vcfdfield2],
+        dxes: dx_lists2_t,
+        ) -> complex:
+    e_a, h_a = mode_a
+    de_b, dh_b = derivative_b
+    return 0.5 * (
+        inner_product(e_a, dh_b, dxes=dxes, conj_h=False)
+        + inner_product(de_b, h_a, dxes=dxes, conj_h=False)
+        )
+
+
+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 = []
+    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 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))
+    return aligned
+
+
+def _section_branches(
+        section: ModeSection,
+        index: int,
+        expected_count: int | None,
+        expected_shape: tuple[int, ...] | None,
+        ) -> tuple[
+            float,
+            list[tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]],
+            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)
+    if not numpy.isfinite(z_coord):
+        raise ValueError(f'sections[{index}].z must be finite')
+
+    shape, _h_shape = _validate_port_modes(f'sections[{index}].modes', section.modes, section.wavenumbers)
+    wavenumbers = _as_wavenumber_array(f'sections[{index}].wavenumbers', section.wavenumbers)
+    if expected_count is not None and len(wavenumbers) != expected_count:
+        raise ValueError('all taper sections must contain the same number of modes')
+    if expected_shape is not None and shape != expected_shape:
+        raise ValueError('all taper section modal fields must share the same E/H shapes')
+
+    if (section.backward_modes is None) != (section.backward_wavenumbers is None):
+        raise ValueError('backward_modes and backward_wavenumbers must be supplied together')
+
+    forward_modes = _as_mode_arrays(section.modes)
+    if section.backward_modes is None:
+        backward_modes = [(e_field.copy(), -h_field.copy()) for e_field, h_field in forward_modes]
+        backward_wavenumbers = -wavenumbers
+    else:
+        backward_shape, _backward_h_shape = _validate_port_modes(
+            f'sections[{index}].backward_modes',
+            section.backward_modes,
+            section.backward_wavenumbers,
+            )
+        if backward_shape != shape:
+            raise ValueError('forward and backward modal fields must share the same E/H shapes')
+        backward_wavenumbers = _as_wavenumber_array(
+            f'sections[{index}].backward_wavenumbers',
+            section.backward_wavenumbers,
+            )
+        backward_modes = _as_mode_arrays(section.backward_modes)
+
+    if len(backward_wavenumbers) != len(wavenumbers):
+        raise ValueError('forward and backward mode counts must match')
+
+    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(
+        sections: Sequence[ModeSection],
+        dxes: dx_lists2_t,
+        ) -> 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,
+            ]:
+    if len(sections) < 2:
+        raise ValueError('at least two taper sections are required')
+
+    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, dual_modes, dual_backward_modes, shape = _section_branches(
+            section,
+            index,
+            expected_count,
+            expected_shape,
+            )
+        if expected_count is None:
+            expected_count = len(forward_wavenumbers)
+            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_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_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],
+        dz: float,
+        dxes: dx_lists2_t,
+        ) -> NDArray[numpy.complex128]:
+    mode_count = len(left_modes)
+    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_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)):
+            derivative = (
+                (right_col_mode[0] - left_col_mode[0]) / dz,
+                (right_col_mode[1] - left_col_mode[1]) / dz,
+                )
+            derivative_overlap[row, col] = _lorentz_derivative_overlap(left_row_mode, derivative, dxes)
+
+    coupling = numpy.linalg.pinv(gram) @ derivative_overlap
+    propagation = numpy.diag(-1j * 0.5 * (left_wavenumbers + right_wavenumbers))
+    return propagation - coupling
+
+
+def _abcd_to_s(
+        abcd: NDArray[numpy.complex128],
+        n_modes: int,
+        ) -> NDArray[numpy.complex128]:
+    A = abcd[:n_modes, :n_modes]
+    B = abcd[:n_modes, n_modes:]
+    C = abcd[n_modes:, :n_modes]
+    D = abcd[n_modes:, n_modes:]
+    D_inv = numpy.linalg.pinv(D)
+    r12 = -D_inv @ C
+    t21 = D_inv
+    t12 = A - B @ D_inv @ C
+    r21 = B @ D_inv
+    return numpy.block([[r12, t12],
+                        [t21, r21]])
+
+
 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.
@@ -77,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
@@ -90,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)
@@ -98,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)
@@ -119,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:
     """
@@ -135,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
@@ -152,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]:
     """
@@ -172,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]])
@@ -188,3 +496,93 @@ def get_s(
         ss = 0.5 * (ss + ss.T)
 
     return ss
+
+
+def get_taper_abcd(
+        sections: Sequence[ModeSection],
+        dxes: dx_lists2_t,
+        *,
+        rtol: float = 1e-7,
+        atol: float = 1e-9,
+        max_step: float | None = None,
+        ) -> sparse.sparray:
+    """
+    Build a bidirectional transfer matrix for a continuously varying taper.
+
+    The taper is represented by local modal bases sampled at increasing `z`
+    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`.
+        dxes: Two-dimensional Yee-cell edge lengths for the shared port plane.
+        rtol: Relative tolerance reserved for future adaptive CMT integrators.
+            Must be positive.
+        atol: Absolute tolerance reserved for future adaptive CMT integrators.
+            Must be positive.
+        max_step: Optional maximum matrix-exponential step inside each sampled
+            interval. This does not change the piecewise-constant interval
+            generator, but can improve conditioning for long intervals.
+
+    Returns:
+        Sparse block transfer matrix ordered as `[forward, backward]`.
+    """
+    if rtol <= 0:
+        raise ValueError('rtol must be positive')
+    if atol <= 0:
+        raise ValueError('atol must be positive')
+    if max_step is not None and max_step <= 0:
+        raise ValueError('max_step must be positive')
+
+    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],
+            float(dz),
+            dxes,
+            )
+        step_count = 1 if max_step is None else max(1, int(numpy.ceil(dz / max_step)))
+        interval_transfer = linalg.expm(generator * (dz / step_count))
+        for _step in range(step_count):
+            transfer = interval_transfer @ transfer
+
+    return sparse.csr_array(transfer)
+
+
+def get_taper_s(
+        sections: Sequence[ModeSection],
+        dxes: dx_lists2_t,
+        *,
+        force_nogain: bool = False,
+        force_reciprocal: bool = False,
+        **kwargs,
+        ) -> NDArray[numpy.complex128]:
+    """
+    Build the full block scattering matrix for a continuously varying taper.
+
+    The returned matrix uses the same ordering as `get_s(...)`:
+    `[[R12, T12], [T21, R21]]`.
+    """
+    _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)
+
+    if force_nogain:
+        U, sing, Vh = numpy.linalg.svd(ss)
+        ss = U @ numpy.diag(numpy.minimum(sing, 1.0)) @ Vh
+
+    if force_reciprocal:
+        ss = 0.5 * (ss + ss.T)
+
+    return ss
diff --git a/meanas/fdfd/waveguide_cyl.py b/meanas/fdfd/waveguide_cyl.py
index f2cb5c3..0d1d4d7 100644
--- a/meanas/fdfd/waveguide_cyl.py
+++ b/meanas/fdfd/waveguide_cyl.py
@@ -43,39 +43,9 @@ T_b &= \operatorname{diag}(r_b / r_{\min}).
 $$
 
 With the same forward/backward derivative notation used in `waveguide_2d`, the
-coordinate-transformed discrete curl equations used here are
-
-$$
-\begin{aligned}
--\imath \omega \mu_{rr} H_r &= \tilde{\partial}_y E_z + \imath \beta T_a^{-1} E_y, \\
--\imath \omega \mu_{yy} H_y &= -\imath \beta T_b^{-1} E_r
-  - T_b^{-1} \tilde{\partial}_r (T_a E_z), \\
--\imath \omega \mu_{zz} H_z &= \tilde{\partial}_r E_y - \tilde{\partial}_y E_r, \\
-\imath \beta H_y &= -\imath \omega T_b \epsilon_{rr} E_r - T_b \hat{\partial}_y H_z, \\
-\imath \beta H_r &= \imath \omega T_a \epsilon_{yy} E_y
-  - T_b T_a^{-1} \hat{\partial}_r (T_b H_z), \\
-\imath \omega E_z &= T_a \epsilon_{zz}^{-1}
-  \left(\hat{\partial}_r H_y - \hat{\partial}_y H_r\right).
-\end{aligned}
-$$
-
-The first three equations are the cylindrical analogue of the straight-guide
-relations for `H_r`, `H_y`, and `H_z`. The next two are the metric-weighted
-versions of the straight-guide identities for `\imath \beta H_y` and
-`\imath \beta H_r`, and the last equation plays the same role as the
-longitudinal `E_z` reconstruction in `waveguide_2d`.
-
-Following the same elimination steps as in `waveguide_2d`, apply
-`\imath \beta \tilde{\partial}_r` and `\imath \beta \tilde{\partial}_y` to the
-equation for `E_z`, substitute for `\imath \beta H_r` and `\imath \beta H_y`,
-and then eliminate `H_z` with
-
-$$
-H_z = \frac{1}{-\imath \omega \mu_{zz}}
-\left(\tilde{\partial}_r E_y - \tilde{\partial}_y E_r\right).
-$$
-
-This yields the transverse electric eigenproblem implemented by
+implementation treats the transverse electric eigenproblem as the canonical
+cylindrical discretization. It reduces to `waveguide_2d.operator_e(...)` in the
+large-radius limit `T_a, T_b \to I`. The eigenproblem implemented by
 `cylindrical_operator(...)`:
 
 $$
@@ -111,6 +81,33 @@ T_a \epsilon_{zz}^{-1}
 \begin{bmatrix} E_r \\ E_y \end{bmatrix}.
 $$
 
+The full fields reconstructed by `exy2e(...)` and `e2h(...)` use the matching
+large-radius-compatible identities
+
+$$
+E_z =
+\frac{1}{\imath \beta} T_a \epsilon_{zz}^{-1}
+\begin{bmatrix}
+\hat{\partial}_r T_b \epsilon_{rr} &
+\hat{\partial}_y T_a \epsilon_{yy}
+\end{bmatrix}
+\begin{bmatrix} E_r \\ E_y \end{bmatrix},
+$$
+
+and
+
+$$
+\begin{bmatrix} H_r \\ H_y \\ H_z \end{bmatrix}
+=
+\frac{1}{-\imath \omega}\mu^{-1}
+\begin{bmatrix}
+0 & \imath\beta T_a^{-1} & \tilde{\partial}_y \\
+-\imath\beta T_b^{-1} & 0 & -T_b^{-1}\tilde{\partial}_r T_a \\
+-\tilde{\partial}_y & \tilde{\partial}_r & 0
+\end{bmatrix}
+\begin{bmatrix} E_r \\ E_y \\ E_z \end{bmatrix}.
+$$
+
 Since `\beta = m / r_{\min}`, the solver implemented in this file returns the
 angular wavenumber `m`, while the operator itself is most naturally written in
 terms of the linear quantity `\beta`. The helpers below reconstruct the full
@@ -143,6 +140,7 @@ def cylindrical_operator(
         dxes: dx_lists2_t,
         epsilon: vfdslice,
         rmin: float,
+        mu: vfdslice | None = None,
         ) -> sparse.sparray:
     r"""
     Cylindrical coordinate waveguide operator of the form
@@ -176,10 +174,13 @@ def cylindrical_operator(
         dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
         epsilon: Vectorized dielectric constant grid
         rmin: Radius at the left edge of the simulation domain (at minimum 'x')
+        mu: Vectorized magnetic permeability grid (default 1 everywhere)
 
     Returns:
         Sparse matrix representation of the operator
     """
+    if mu is None:
+        mu = numpy.ones_like(epsilon)
 
     Dfx, Dfy = deriv_forward(dxes[0])
     Dbx, Dby = deriv_back(dxes[1])
@@ -191,12 +192,17 @@ def cylindrical_operator(
     eps_y = sparse.diags_array(eps_parts[1])
     eps_z_inv = sparse.diags_array(1 / eps_parts[2])
 
+    mu_parts = numpy.split(mu, 3)
+    mu_y = sparse.diags_array(mu_parts[1])
+    mu_x = sparse.diags_array(mu_parts[0])
+    mu_z_inv = sparse.diags_array(1 / mu_parts[2])
+
     omega2 = omega * omega
     diag = sparse.block_diag
 
-    sq0 = omega2 * diag((Tb @ Tb @ eps_x,
-                         Ta @ Ta @ eps_y))
-    lin0 = sparse.vstack((-Tb @ Dby, Ta @ Dbx)) @ Tb @ sparse.hstack((-Dfy, Dfx))
+    sq0 = omega2 * diag((Tb @ Tb @ mu_y @ eps_x,
+                         Ta @ Ta @ mu_x @ eps_y))
+    lin0 = sparse.vstack((-Tb @ mu_y @ Dby, Ta @ mu_x @ Dbx)) @ Tb @ mu_z_inv @ sparse.hstack((-Dfy, Dfx))
     lin1 = sparse.vstack((Dfx, Dfy)) @ Ta @ eps_z_inv @ sparse.hstack((Dbx @ Tb @ eps_x,
                                                                        Dby @ Ta @ eps_y))
     op = sq0 + lin0 + lin1
@@ -209,6 +215,7 @@ def solve_modes(
         dxes: dx_lists2_t,
         epsilon: vfdslice,
         rmin: float,
+        mu: vfdslice | None = None,
         mode_margin: int = 2,
         ) -> tuple[NDArray[numpy.complex128], NDArray[numpy.complex128]]:
     """
@@ -223,6 +230,7 @@ def solve_modes(
         epsilon: Dielectric constant
         rmin: Radius of curvature for the simulation. This should be the minimum value of
             r within the simulation domain.
+        mu: Magnetic permeability (default 1 everywhere)
 
     Returns:
         e_xys: NDArray of vfdfield_t specifying fields. First dimension is mode number.
@@ -233,8 +241,9 @@ def solve_modes(
     # Solve for the largest-magnitude eigenvalue of the real operator
     #
     dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes]
+    mu_real = None if mu is None else numpy.real(mu)
 
-    A_r = cylindrical_operator(numpy.real(omega), dxes_real, numpy.real(epsilon), rmin=rmin)
+    A_r = cylindrical_operator(numpy.real(omega), dxes_real, numpy.real(epsilon), rmin=rmin, mu=mu_real)
     eigvals, eigvecs = signed_eigensolve(A_r, max(mode_numbers) + mode_margin)
     keep_inds = -(numpy.array(mode_numbers) + 1)
     e_xys = eigvecs[:, keep_inds].T
@@ -244,7 +253,7 @@ def solve_modes(
     # Now solve for the eigenvector of the full operator, using the real operator's
     #  eigenvector as an initial guess for Rayleigh quotient iteration.
     #
-    A = cylindrical_operator(omega, dxes, epsilon, rmin=rmin)
+    A = cylindrical_operator(omega, dxes, epsilon, rmin=rmin, mu=mu)
     for nn in range(len(mode_numbers)):
         eigvals[nn], e_xys[nn, :] = rayleigh_quotient_iteration(A, e_xys[nn, :])
 
@@ -312,12 +321,20 @@ def linear_wavenumbers(
 
     shape2d = (len(dxes[0][0]), len(dxes[0][1]))
     epsilon2d = unvec(epsilon, shape2d)[:2]
-    grid_radii = rmin + numpy.cumsum(dxes[0][0])
+    ra = rmin + numpy.cumsum(dxes[0][0])
+    rb = rmin + dxes[0][0] / 2.0 + numpy.concatenate((
+        numpy.zeros(1, dtype=dxes[1][0].dtype),
+        numpy.cumsum(dxes[1][0][:-1]),
+        ))
     for ii in range(angular_wavenumbers.size):
         efield = unvec(e_xys[ii], shape2d, 2)
         energy = numpy.real((efield * efield.conj()) * epsilon2d)
-        energy_vs_x = energy.sum(axis=(0, 2))
-        mode_radii[ii] = (grid_radii * energy_vs_x).sum() / energy_vs_x.sum()
+        er_energy_vs_r = energy[0].sum(axis=1)
+        ey_energy_vs_r = energy[1].sum(axis=1)
+        energy_vs_r = er_energy_vs_r + ey_energy_vs_r
+        mode_radii[ii] = (
+            (rb * er_energy_vs_r).sum() + (ra * ey_energy_vs_r).sum()
+            ) / energy_vs_r.sum()
 
     logger.info(f'{mode_radii=}')
     lin_wavenumbers = angular_wavenumbers / mode_radii
@@ -350,12 +367,11 @@ def exy2h(
         Sparse matrix representing the operator.
     """
     e2hop = e2h(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, mu=mu)
-    return e2hop @ exy2e(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon)
+    return e2hop @ exy2e(angular_wavenumber=angular_wavenumber, dxes=dxes, rmin=rmin, epsilon=epsilon)
 
 
 def exy2e(
         angular_wavenumber: complex,
-        omega: float,
         dxes: dx_lists2_t,
         rmin: float,
         epsilon: vfdslice,
@@ -371,7 +387,6 @@ def exy2e(
         angular_wavenumber: Wavenumber assuming fields have theta-dependence of
             `exp(-i * angular_wavenumber * theta)`. It should satisfy
             `operator_e() @ e_xy == (angular_wavenumber / rmin) ** 2 * e_xy`
-        omega: The angular frequency of the system
         dxes: Grid parameters `[dx_e, dx_h]` as described in `meanas.fdmath.types` (2D)
         rmin: Radius at the left edge of the simulation domain (at minimum 'x')
         epsilon: Vectorized dielectric constant grid
@@ -379,30 +394,22 @@ def exy2e(
     Returns:
         Sparse matrix representing the operator.
     """
-    Dfx, Dfy = deriv_forward(dxes[0])
     Dbx, Dby = deriv_back(dxes[1])
     wavenumber = angular_wavenumber / rmin
 
     Ta, Tb = dxes2T(dxes=dxes, rmin=rmin)
-    Tai = sparse.diags_array(1 / Ta.diagonal())
-    #Tbi = sparse.diags_array(1 / Tb.diagonal())
 
     epsilon_parts = numpy.split(epsilon, 3)
     epsilon_x, epsilon_y = (sparse.diags_array(epsi) for epsi in epsilon_parts[:2])
     epsilon_z_inv = sparse.diags_array(1 / epsilon_parts[2])
 
     n_pts = dxes[0][0].size * dxes[0][1].size
-    zeros = sparse.coo_array((n_pts, n_pts))
-
-    mu_z = numpy.ones(n_pts)
-    mu_z_inv = sparse.diags_array(1 / mu_z)
-    exy2hz = 1 / (-1j * omega) * mu_z_inv @ sparse.hstack((Dfy, -Dfx))
-    hxy2ez = 1 / (1j * omega) * epsilon_z_inv @ sparse.hstack((Dby, -Dbx))
-
-    exy2hy = Tb / (1j * wavenumber) @ (-1j * omega * sparse.hstack((epsilon_x, zeros)) - Dby @ exy2hz)
-    exy2hx = Tb / (1j * wavenumber) @ ( 1j * omega * sparse.hstack((zeros, epsilon_y)) - Tai @ Dbx @ Tb @ exy2hz)
-
-    exy2ez = hxy2ez @ sparse.vstack((exy2hx, exy2hy))
+    exy2ez = (
+        Ta @ epsilon_z_inv
+        @ sparse.hstack((Dbx @ Tb @ epsilon_x,
+                         Dby @ Ta @ epsilon_y))
+        / (1j * wavenumber)
+        )
 
     op = sparse.vstack((sparse.eye_array(2 * n_pts),
                         exy2ez))
@@ -448,9 +455,9 @@ def e2h(
     Tbi = sparse.diags_array(1 / Tb.diagonal())
 
     jB = 1j * angular_wavenumber / rmin
-    op = sparse.block_array([[    None, -jB * Tai,           -Dfy],
-                             [jB * Tbi,      None, Tbi @ Dfx @ Ta],
-                             [     Dfy,      -Dfx,           None]]) / (-1j * omega)
+    op = sparse.block_array([[     None,  jB * Tai,            Dfy],
+                             [-jB * Tbi,       None, -Tbi @ Dfx @ Ta],
+                             [     -Dfy,        Dfx,           None]]) / (-1j * omega)
     if mu is not None:
         op = sparse.diags_array(1 / mu) @ op
     return op
@@ -475,7 +482,14 @@ def dxes2T(
         Sparse diagonal matrices `(T_a, T_b)`.
     """
     ra = rmin + numpy.cumsum(dxes[0][0])                      # Radius at Ey points
-    rb = rmin + dxes[0][0] / 2.0 + numpy.cumsum(dxes[1][0])   # Radius at Ex points
+    rb = (
+        rmin
+        + dxes[0][0] / 2.0
+        + numpy.concatenate((
+            numpy.zeros(1, dtype=dxes[1][0].dtype),
+            numpy.cumsum(dxes[1][0][:-1]),
+            ))
+        )                                                     # Radius at Er points
     ta = ra / rmin
     tb = rb / rmin
 
@@ -527,7 +541,7 @@ def normalized_fields_e(
         fields, then the overall complex phase and sign are fixed so the result is
         reproducible for symmetric modes.
     """
-    e = exy2e(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon) @ e_xy
+    e = exy2e(angular_wavenumber=angular_wavenumber, 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,
@@ -553,19 +567,16 @@ def _normalized_fields(
 
     The normalization procedure is:
 
-    1. Flip the magnetic field sign so the reconstructed `(e, h)` pair follows
-       the same forward-power convention as `waveguide_2d`.
-    2. Compute the time-averaged forward power with
+    1. Compute the time-averaged forward power with
        `waveguide_2d.inner_product(..., conj_h=True)`.
-    3. Scale by `1 / sqrt(S_z)` so the resulting mode has unit forward power.
-    4. Remove the arbitrary complex phase and apply a quadrant-sum sign heuristic
+    2. Scale by `1 / sqrt(S_z)` so the resulting mode has unit forward power.
+    3. Remove the arbitrary complex phase and apply a quadrant-sum sign heuristic
        so symmetric modes choose a stable representative.
 
     `prop_phase` has the same meaning as in `waveguide_2d`: it compensates for
     the half-cell longitudinal staggering between the E and H fields when the
     propagation direction is itself discretized.
     """
-    h *= -1
     shape = [s.size for s in dxes[0]]
 
     # Find time-averaged Sz and normalize to it
diff --git a/meanas/test/test_eme_numerics.py b/meanas/test/test_eme_numerics.py
index 3237c1b..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
@@ -227,6 +262,159 @@ def _build_uniform_mode(index: float) -> tuple[tuple[numpy.ndarray, numpy.ndarra
     return (vec(e_field), vec(h_field)), complex(index * OMEGA)
 
 
+def test_get_taper_abcd_constant_section_is_phase_only() -> None:
+    mode, beta = _build_uniform_mode(1.5)
+    length = 11.0
+
+    abcd = eme.get_taper_abcd(
+        [
+            eme.ModeSection(0.0, [mode], [beta]),
+            eme.ModeSection(length, [mode], [beta]),
+        ],
+        dxes=REAL_DXES,
+    ).toarray()
+
+    assert_close(abcd, _propagation_abcd(beta, length), atol=1e-12, rtol=1e-12)
+
+
+def test_get_taper_s_constant_section_has_no_reflection() -> None:
+    mode, beta = _build_uniform_mode(1.5)
+    length = 11.0
+    phase = numpy.exp(-1j * beta * length)
+
+    ss = eme.get_taper_s(
+        [
+            eme.ModeSection(0.0, [mode], [beta]),
+            eme.ModeSection(length, [mode], [beta]),
+        ],
+        dxes=REAL_DXES,
+    )
+
+    assert_close(ss, numpy.array([[0.0, phase], [phase, 0.0]], dtype=complex), atol=1e-12, rtol=1e-12)
+
+
+def test_get_taper_abcd_is_invariant_to_adjacent_modal_phase() -> None:
+    mode, beta = _build_uniform_mode(1.5)
+    shifted_mode = (mode[0] * numpy.exp(0.73j), mode[1] * numpy.exp(0.73j))
+    length = 11.0
+    base_sections = [
+        eme.ModeSection(0.0, [mode], [beta]),
+        eme.ModeSection(length, [mode], [beta]),
+    ]
+    shifted_sections = [
+        eme.ModeSection(0.0, [mode], [beta]),
+        eme.ModeSection(length, [shifted_mode], [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_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)
+
+    with pytest.raises(ValueError, match='strictly increasing'):
+        eme.get_taper_abcd(
+            [
+                eme.ModeSection(0.0, [mode], [beta]),
+                eme.ModeSection(0.0, [mode], [beta]),
+            ],
+            dxes=REAL_DXES,
+        )
+
+
+def test_get_taper_rejects_mode_count_changes() -> None:
+    mode, beta = _build_uniform_mode(1.5)
+
+    with pytest.raises(ValueError, match='same number of modes'):
+        eme.get_taper_abcd(
+            [
+                eme.ModeSection(0.0, [mode], [beta]),
+                eme.ModeSection(1.0, [mode, mode], [beta, beta]),
+            ],
+            dxes=REAL_DXES,
+        )
+
+
+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)
@@ -339,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
diff --git a/meanas/test/test_examples_smoke.py b/meanas/test/test_examples_smoke.py
index b21f90b..b3797de 100644
--- a/meanas/test/test_examples_smoke.py
+++ b/meanas/test/test_examples_smoke.py
@@ -45,3 +45,11 @@ def test_eme_bend_example_smoke_runs(tmp_path: Path) -> None:
     assert result.returncode == 0, result.stdout + result.stderr
     assert 'straight effective indices:' in result.stdout
     assert 'cascaded bend through power' in result.stdout
+
+
+def test_eme_taper_cmt_example_smoke_runs(tmp_path: Path) -> None:
+    result = _run_example('eme_taper_cmt.py', tmp_path)
+
+    assert result.returncode == 0, result.stdout + result.stderr
+    assert 'sampled taper effective indices:' in result.stdout
+    assert 'taper CMT transmission' in result.stdout
diff --git a/meanas/test/test_waveguide_mode_helpers.py b/meanas/test/test_waveguide_mode_helpers.py
index d3ec7cd..dc51de0 100644
--- a/meanas/test/test_waveguide_mode_helpers.py
+++ b/meanas/test/test_waveguide_mode_helpers.py
@@ -35,6 +35,7 @@ def build_waveguide_3d_mode(
 def build_waveguide_cyl_fixture(
         *,
         nonuniform: bool = False,
+        asymmetric: bool = False,
         ) -> tuple[list[list[numpy.ndarray]], numpy.ndarray, float]:
     if nonuniform:
         dxes = [
@@ -43,10 +44,17 @@ def build_waveguide_cyl_fixture(
         ]
     else:
         dxes = [[numpy.ones(5), numpy.ones(5)] for _ in range(2)]
-    epsilon = vec(numpy.ones((3, 5, 5), dtype=float))
+    epsilon_3d = numpy.ones((3, 5, 5), dtype=float)
+    if asymmetric:
+        epsilon_3d[:, 2, 1] = 2.0
+    epsilon = vec(epsilon_3d)
     return dxes, epsilon, 10.0
 
 
+def build_waveguide_cyl_mu_profile() -> numpy.ndarray:
+    return numpy.linspace(1.5, 2.2, 3 * 5 * 5)
+
+
 def test_waveguide_3d_solve_mode_and_expand_e_are_phase_consistent() -> None:
     epsilon, dxes, slices, result = build_waveguide_3d_mode(slice_start=0, polarity=1)
     axis = 0
@@ -173,8 +181,10 @@ def test_waveguide_3d_compute_overlap_e_rejects_zero_support_window() -> None:
             )
 
 
-def test_waveguide_cyl_solved_modes_are_ordered_and_low_residual() -> None:
-    dxes, epsilon, rmin = build_waveguide_cyl_fixture()
+@pytest.mark.parametrize('use_mu', [False, True])
+def test_waveguide_cyl_solved_modes_are_ordered_and_low_residual(use_mu: bool) -> None:
+    dxes, epsilon, rmin = build_waveguide_cyl_fixture(asymmetric=use_mu)
+    mu = build_waveguide_cyl_mu_profile() if use_mu else None
 
     e_xys, angular_wavenumbers = waveguide_cyl.solve_modes(
         [0, 1],
@@ -182,8 +192,9 @@ def test_waveguide_cyl_solved_modes_are_ordered_and_low_residual() -> None:
         dxes=dxes,
         epsilon=epsilon,
         rmin=rmin,
+        mu=mu,
         )
-    operator = waveguide_cyl.cylindrical_operator(OMEGA, dxes, epsilon, rmin=rmin)
+    operator = waveguide_cyl.cylindrical_operator(OMEGA, dxes, epsilon, rmin=rmin, mu=mu)
 
     assert numpy.all(numpy.diff(numpy.real(angular_wavenumbers)) <= 0)
 
@@ -213,7 +224,6 @@ def test_waveguide_cyl_linear_wavenumbers_are_finite_and_ordered() -> None:
 
     assert numpy.isfinite(linear_wavenumbers).all()
     assert numpy.all(numpy.real(linear_wavenumbers) > 0)
-    assert numpy.all(numpy.diff(numpy.real(linear_wavenumbers)) <= 0)
 
 
 def test_waveguide_cyl_dxes2t_matches_expected_radius_scaling() -> None:
@@ -221,26 +231,127 @@ def test_waveguide_cyl_dxes2t_matches_expected_radius_scaling() -> None:
     Ta, Tb = waveguide_cyl.dxes2T(dxes, rmin)
 
     ta = (rmin + numpy.cumsum(dxes[0][0])) / rmin
-    tb = (rmin + dxes[0][0] / 2 + numpy.cumsum(dxes[1][0])) / rmin
+    tb = (
+        rmin
+        + dxes[0][0] / 2
+        + numpy.concatenate((numpy.zeros(1), numpy.cumsum(dxes[1][0][:-1])))
+        ) / rmin
 
     numpy.testing.assert_allclose(Ta.diagonal(), numpy.repeat(ta, dxes[0][1].size))
     numpy.testing.assert_allclose(Tb.diagonal(), numpy.repeat(tb, dxes[1][1].size))
 
 
+@pytest.mark.parametrize('use_mu', [False, True])
+def test_waveguide_cyl_operator_matches_2d_limit(use_mu: bool) -> None:
+    dxes, epsilon, _rmin = build_waveguide_cyl_fixture(asymmetric=True)
+    mu = build_waveguide_cyl_mu_profile() if use_mu else None
+    rmin = 1e15
+
+    cyl_operator = waveguide_cyl.cylindrical_operator(OMEGA, dxes, epsilon, rmin=rmin, mu=mu)
+    straight_operator = waveguide_2d.operator_e(OMEGA, dxes, epsilon, mu=mu)
+
+    numpy.testing.assert_allclose(
+        cyl_operator.toarray(),
+        straight_operator.toarray(),
+        rtol=1e-9,
+        atol=1e-10,
+        )
+
+
+@pytest.mark.parametrize('use_mu', [False, True])
+def test_waveguide_cyl_reconstruction_matches_2d_limit(use_mu: bool) -> None:
+    dxes, epsilon, _rmin = build_waveguide_cyl_fixture(asymmetric=True)
+    mu = build_waveguide_cyl_mu_profile() if use_mu else None
+    rmin = 1e15
+    e_xy, wavenumber = waveguide_2d.solve_mode(
+        0,
+        omega=OMEGA,
+        dxes=dxes,
+        epsilon=epsilon,
+        mu=mu,
+        )
+    angular_wavenumber = wavenumber * rmin
+
+    cyl_e = waveguide_cyl.exy2e(
+        angular_wavenumber=angular_wavenumber,
+        dxes=dxes,
+        rmin=rmin,
+        epsilon=epsilon,
+        ) @ e_xy
+    straight_e = waveguide_2d.exy2e(
+        wavenumber=wavenumber,
+        dxes=dxes,
+        epsilon=epsilon,
+        ) @ e_xy
+    cyl_h = waveguide_cyl.e2h(
+        angular_wavenumber=angular_wavenumber,
+        omega=OMEGA,
+        dxes=dxes,
+        rmin=rmin,
+        mu=mu,
+        ) @ cyl_e
+    straight_h = waveguide_2d.e2h(
+        wavenumber=wavenumber,
+        omega=OMEGA,
+        dxes=dxes,
+        mu=mu,
+        ) @ straight_e
+
+    numpy.testing.assert_allclose(cyl_e, straight_e, rtol=1e-8, atol=1e-8)
+    numpy.testing.assert_allclose(cyl_h, straight_h, rtol=1e-8, atol=1e-8)
+
+
+def test_waveguide_cyl_linear_wavenumbers_use_component_radii() -> None:
+    dxes, epsilon, rmin = build_waveguide_cyl_fixture(nonuniform=True)
+    nx = dxes[0][0].size
+    ny = dxes[0][1].size
+    angular_wavenumber = numpy.array([2.0])
+
+    ra = rmin + numpy.cumsum(dxes[0][0])
+    rb = (
+        rmin
+        + dxes[0][0] / 2
+        + numpy.concatenate((numpy.zeros(1), numpy.cumsum(dxes[1][0][:-1])))
+        )
+
+    er_only = numpy.zeros((1, 2 * nx * ny), dtype=complex)
+    er_only[0] = vec(numpy.array([numpy.ones((nx, ny)), numpy.zeros((nx, ny))]))
+    ey_only = numpy.zeros_like(er_only)
+    ey_only[0] = vec(numpy.array([numpy.zeros((nx, ny)), numpy.ones((nx, ny))]))
+
+    er_linear = waveguide_cyl.linear_wavenumbers(
+        er_only,
+        angular_wavenumber,
+        epsilon=epsilon,
+        dxes=dxes,
+        rmin=rmin,
+        )
+    ey_linear = waveguide_cyl.linear_wavenumbers(
+        ey_only,
+        angular_wavenumber,
+        epsilon=epsilon,
+        dxes=dxes,
+        rmin=rmin,
+        )
+
+    numpy.testing.assert_allclose(er_linear[0], angular_wavenumber[0] / rb.mean())
+    numpy.testing.assert_allclose(ey_linear[0], angular_wavenumber[0] / ra.mean())
+
+
 def test_waveguide_cyl_exy2e_and_exy2h_return_finite_full_fields() -> None:
     dxes, epsilon, rmin = build_waveguide_cyl_fixture()
-    mu = vec(2 * numpy.ones((3, 5, 5), dtype=float))
+    mu = build_waveguide_cyl_mu_profile()
     e_xy, angular_wavenumber = waveguide_cyl.solve_mode(
         0,
         omega=OMEGA,
         dxes=dxes,
         epsilon=epsilon,
         rmin=rmin,
+        mu=mu,
         )
 
     e_field = waveguide_cyl.exy2e(
         angular_wavenumber=angular_wavenumber,
-        omega=OMEGA,
         dxes=dxes,
         rmin=rmin,
         epsilon=epsilon,
@@ -265,13 +376,14 @@ def test_waveguide_cyl_exy2e_and_exy2h_return_finite_full_fields() -> None:
 @pytest.mark.parametrize('use_mu', [False, True])
 def test_waveguide_cyl_normalized_fields_are_unit_norm_and_silent(use_mu: bool) -> None:
     dxes, epsilon, rmin = build_waveguide_cyl_fixture()
-    mu = vec(2 * numpy.ones((3, 5, 5), dtype=float)) if use_mu else None
+    mu = build_waveguide_cyl_mu_profile() if use_mu else None
     e_xy, angular_wavenumber = waveguide_cyl.solve_mode(
         0,
         omega=OMEGA,
         dxes=dxes,
         epsilon=epsilon,
         rmin=rmin,
+        mu=mu,
         )
 
     output = io.StringIO()