diff --git a/README.md b/README.md
index d179c9d..73d48b5 100644
--- a/README.md
+++ b/README.md
@@ -172,8 +172,6 @@ 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 1b5ca6b..5475af8 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -28,8 +28,6 @@ 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
deleted file mode 100644
index c5a8f1c..0000000
--- a/examples/eme_taper_cmt.py
+++ /dev/null
@@ -1,134 +0,0 @@
-"""
-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 a9de4dd..af745e8 100644
--- a/meanas/fdfd/eme.py
+++ b/meanas/fdfd/eme.py
@@ -13,9 +13,6 @@ 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
@@ -23,51 +20,19 @@ 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: mode_seq,
+        ehs: Sequence[Sequence[vcfdfield2]],
         wavenumbers: wavenumber_seq,
         ) -> tuple[tuple[int, ...], tuple[int, ...]]:
     if len(ehs) != len(wavenumbers):
@@ -96,274 +61,12 @@ 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: mode_seq,
+        ehLs: Sequence[Sequence[vcfdfield2]],
         wavenumbers_L: wavenumber_seq,
-        ehRs: mode_seq,
+        ehRs: Sequence[Sequence[vcfdfield2]],
         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.
@@ -374,8 +77,6 @@ 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
@@ -389,8 +90,6 @@ 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)
@@ -399,12 +98,11 @@ def get_tr(
     B11 = numpy.zeros((nL,), dtype=complex)
     for ll in range(nL):
         eL, hL = ehLs[ll]
-        eP, hP = projection_ehLs[ll]
-        B11[ll] = inner_product(eL, hP, dxes=dxes, conj_h=False)
+        B11[ll] = inner_product(eL, hL, dxes=dxes, conj_h=False)
         for rr in range(nR):
             eR, hR = ehRs[rr]
-            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)
+            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)
 
     # tt0 = 2 * numpy.linalg.pinv(A21 + numpy.conj(A12))
     tt0, _resid, _rank, _sing = numpy.linalg.lstsq(A21 + A12, numpy.diag(2 * B11), rcond=None)
@@ -421,12 +119,10 @@ def get_tr(
 
 
 def get_abcd(
-        ehLs: mode_seq,
+        ehLs: Sequence[Sequence[vcfdfield2]],
         wavenumbers_L: wavenumber_seq,
-        ehRs: mode_seq,
+        ehRs: Sequence[Sequence[vcfdfield2]],
         wavenumbers_R: wavenumber_seq,
-        dual_ehLs: mode_seq | None = None,
-        dual_ehRs: mode_seq | None = None,
         **kwargs,
         ) -> sparse.sparray:
     """
@@ -439,8 +135,8 @@ def get_abcd(
 
     convention.
     """
-    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)
+    t12, r12 = get_tr(ehLs, wavenumbers_L, ehRs, wavenumbers_R, **kwargs)
+    t21, r21 = get_tr(ehRs, wavenumbers_R, ehLs, wavenumbers_L, **kwargs)
     t21i = numpy.linalg.pinv(t21)
     A = t12 - r21 @ t21i @ r12
     B = r21 @ t21i
@@ -456,14 +152,12 @@ def get_abcd(
 
 
 def get_s(
-        ehLs: mode_seq,
+        ehLs: Sequence[Sequence[vcfdfield2]],
         wavenumbers_L: wavenumber_seq,
-        ehRs: mode_seq,
+        ehRs: Sequence[Sequence[vcfdfield2]],
         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]:
     """
@@ -478,11 +172,9 @@ 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, dual_ehLs=dual_ehLs, **kwargs)
-    t21, r21 = get_tr(ehRs, wavenumbers_R, ehLs, wavenumbers_L, dual_ehLs=dual_ehRs, **kwargs)
+    t12, r12 = get_tr(ehLs, wavenumbers_L, ehRs, wavenumbers_R, **kwargs)
+    t21, r21 = get_tr(ehRs, wavenumbers_R, ehLs, wavenumbers_L, **kwargs)
 
     ss = numpy.block([[r12, t12],
                       [t21, r21]])
@@ -496,93 +188,3 @@ 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 0d1d4d7..f2cb5c3 100644
--- a/meanas/fdfd/waveguide_cyl.py
+++ b/meanas/fdfd/waveguide_cyl.py
@@ -43,9 +43,39 @@ T_b &= \operatorname{diag}(r_b / r_{\min}).
 $$
 
 With the same forward/backward derivative notation used in `waveguide_2d`, the
-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
+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
 `cylindrical_operator(...)`:
 
 $$
@@ -81,33 +111,6 @@ 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
@@ -140,7 +143,6 @@ 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
@@ -174,13 +176,10 @@ 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])
@@ -192,17 +191,12 @@ 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 @ 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))
+    sq0 = omega2 * diag((Tb @ Tb @ eps_x,
+                         Ta @ Ta @ eps_y))
+    lin0 = sparse.vstack((-Tb @ Dby, Ta @ Dbx)) @ Tb @ 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
@@ -215,7 +209,6 @@ 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]]:
     """
@@ -230,7 +223,6 @@ 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.
@@ -241,9 +233,8 @@ 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, mu=mu_real)
+    A_r = cylindrical_operator(numpy.real(omega), dxes_real, numpy.real(epsilon), rmin=rmin)
     eigvals, eigvecs = signed_eigensolve(A_r, max(mode_numbers) + mode_margin)
     keep_inds = -(numpy.array(mode_numbers) + 1)
     e_xys = eigvecs[:, keep_inds].T
@@ -253,7 +244,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, mu=mu)
+    A = cylindrical_operator(omega, dxes, epsilon, rmin=rmin)
     for nn in range(len(mode_numbers)):
         eigvals[nn], e_xys[nn, :] = rayleigh_quotient_iteration(A, e_xys[nn, :])
 
@@ -321,20 +312,12 @@ def linear_wavenumbers(
 
     shape2d = (len(dxes[0][0]), len(dxes[0][1]))
     epsilon2d = unvec(epsilon, shape2d)[:2]
-    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]),
-        ))
+    grid_radii = rmin + numpy.cumsum(dxes[0][0])
     for ii in range(angular_wavenumbers.size):
         efield = unvec(e_xys[ii], shape2d, 2)
         energy = numpy.real((efield * efield.conj()) * epsilon2d)
-        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()
+        energy_vs_x = energy.sum(axis=(0, 2))
+        mode_radii[ii] = (grid_radii * energy_vs_x).sum() / energy_vs_x.sum()
 
     logger.info(f'{mode_radii=}')
     lin_wavenumbers = angular_wavenumbers / mode_radii
@@ -367,11 +350,12 @@ 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, dxes=dxes, rmin=rmin, epsilon=epsilon)
+    return e2hop @ exy2e(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon)
 
 
 def exy2e(
         angular_wavenumber: complex,
+        omega: float,
         dxes: dx_lists2_t,
         rmin: float,
         epsilon: vfdslice,
@@ -387,6 +371,7 @@ 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
@@ -394,22 +379,30 @@ 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
-    exy2ez = (
-        Ta @ epsilon_z_inv
-        @ sparse.hstack((Dbx @ Tb @ epsilon_x,
-                         Dby @ Ta @ epsilon_y))
-        / (1j * wavenumber)
-        )
+    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))
 
     op = sparse.vstack((sparse.eye_array(2 * n_pts),
                         exy2ez))
@@ -455,9 +448,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
@@ -482,14 +475,7 @@ 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.concatenate((
-            numpy.zeros(1, dtype=dxes[1][0].dtype),
-            numpy.cumsum(dxes[1][0][:-1]),
-            ))
-        )                                                     # Radius at Er points
+    rb = rmin + dxes[0][0] / 2.0 + numpy.cumsum(dxes[1][0])   # Radius at Ex points
     ta = ra / rmin
     tb = rb / rmin
 
@@ -541,7 +527,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, dxes=dxes, rmin=rmin, epsilon=epsilon) @ e_xy
+    e = exy2e(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon) @ e_xy
     h = exy2h(angular_wavenumber=angular_wavenumber, omega=omega, dxes=dxes, rmin=rmin, epsilon=epsilon, mu=mu) @ e_xy
     e_norm, h_norm = _normalized_fields(
         e=e, h=h, dxes=dxes, epsilon=epsilon, prop_phase=prop_phase,
@@ -567,16 +553,19 @@ def _normalized_fields(
 
     The normalization procedure is:
 
-    1. Compute the time-averaged forward power with
+    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
        `waveguide_2d.inner_product(..., conj_h=True)`.
-    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
+    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
        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 0d28692..3237c1b 100644
--- a/meanas/test/test_eme_numerics.py
+++ b/meanas/test/test_eme_numerics.py
@@ -77,27 +77,6 @@ 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)
@@ -187,20 +166,6 @@ 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
@@ -262,159 +227,6 @@ 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)
@@ -527,34 +339,6 @@ 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 b3797de..b21f90b 100644
--- a/meanas/test/test_examples_smoke.py
+++ b/meanas/test/test_examples_smoke.py
@@ -45,11 +45,3 @@ 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 dc51de0..d3ec7cd 100644
--- a/meanas/test/test_waveguide_mode_helpers.py
+++ b/meanas/test/test_waveguide_mode_helpers.py
@@ -35,7 +35,6 @@ 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 = [
@@ -44,17 +43,10 @@ def build_waveguide_cyl_fixture(
         ]
     else:
         dxes = [[numpy.ones(5), numpy.ones(5)] for _ in range(2)]
-    epsilon_3d = numpy.ones((3, 5, 5), dtype=float)
-    if asymmetric:
-        epsilon_3d[:, 2, 1] = 2.0
-    epsilon = vec(epsilon_3d)
+    epsilon = vec(numpy.ones((3, 5, 5), dtype=float))
     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
@@ -181,10 +173,8 @@ def test_waveguide_3d_compute_overlap_e_rejects_zero_support_window() -> None:
             )
 
 
-@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
+def test_waveguide_cyl_solved_modes_are_ordered_and_low_residual() -> None:
+    dxes, epsilon, rmin = build_waveguide_cyl_fixture()
 
     e_xys, angular_wavenumbers = waveguide_cyl.solve_modes(
         [0, 1],
@@ -192,9 +182,8 @@ def test_waveguide_cyl_solved_modes_are_ordered_and_low_residual(use_mu: bool) -
         dxes=dxes,
         epsilon=epsilon,
         rmin=rmin,
-        mu=mu,
         )
-    operator = waveguide_cyl.cylindrical_operator(OMEGA, dxes, epsilon, rmin=rmin, mu=mu)
+    operator = waveguide_cyl.cylindrical_operator(OMEGA, dxes, epsilon, rmin=rmin)
 
     assert numpy.all(numpy.diff(numpy.real(angular_wavenumbers)) <= 0)
 
@@ -224,6 +213,7 @@ 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:
@@ -231,127 +221,26 @@ 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.concatenate((numpy.zeros(1), numpy.cumsum(dxes[1][0][:-1])))
-        ) / rmin
+    tb = (rmin + dxes[0][0] / 2 + numpy.cumsum(dxes[1][0])) / 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 = build_waveguide_cyl_mu_profile()
+    mu = vec(2 * numpy.ones((3, 5, 5), dtype=float))
     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,
@@ -376,14 +265,13 @@ 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 = build_waveguide_cyl_mu_profile() if use_mu else None
+    mu = vec(2 * numpy.ones((3, 5, 5), dtype=float)) 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()