diff --git a/meanas/fdmath/vectorization.py b/meanas/fdmath/vectorization.py
index f2c01d0..44d8b74 100644
--- a/meanas/fdmath/vectorization.py
+++ b/meanas/fdmath/vectorization.py
@@ -18,35 +18,35 @@ from .types import (
 
 @overload
 def vec(f: None) -> None:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def vec(f: fdfield_t) -> vfdfield_t:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def vec(f: cfdfield_t) -> vcfdfield_t:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def vec(f: fdfield2_t) -> vfdfield2_t:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def vec(f: cfdfield2_t) -> vcfdfield2_t:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def vec(f: fdslice_t) -> vfdslice_t:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def vec(f: cfdslice_t) -> vcfdslice_t:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def vec(f: ArrayLike) -> NDArray:
-    pass
+    pass  # pragma: no cover
 
 def vec(
         f: fdfield_t | cfdfield_t | fdfield2_t | cfdfield2_t | fdslice_t | cfdslice_t | ArrayLike | None,
@@ -70,15 +70,15 @@ def vec(
 
 @overload
 def unvec(v: None, shape: Sequence[int], nvdim: int = 3) -> None:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def unvec(v: vfdfield_t, shape: Sequence[int], nvdim: int = 3) -> fdfield_t:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def unvec(v: vcfdfield_t, shape: Sequence[int], nvdim: int = 3) -> cfdfield_t:
-    pass
+    pass  # pragma: no cover
 
 @overload
 def unvec(v: vfdfield2_t, shape: Sequence[int], nvdim: int = 3) -> fdfield2_t:
diff --git a/meanas/test/test_fdmath_vectorization.py b/meanas/test/test_fdmath_vectorization.py
new file mode 100644
index 0000000..33a9812
--- /dev/null
+++ b/meanas/test/test_fdmath_vectorization.py
@@ -0,0 +1,45 @@
+import numpy
+from numpy.testing import assert_allclose, assert_array_equal
+
+from ..fdmath import unvec, vec
+
+
+SHAPE = (2, 3, 2)
+FIELD = numpy.arange(3 * numpy.prod(SHAPE), dtype=float).reshape((3, *SHAPE), order='C')
+COMPLEX_FIELD = (FIELD + 0.5j * FIELD).astype(complex)
+
+
+def test_vec_and_unvec_return_none_for_none_input() -> None:
+    assert vec(None) is None
+    assert unvec(None, SHAPE) is None
+
+
+def test_real_field_round_trip_preserves_shape_and_values() -> None:
+    vector = vec(FIELD)
+    assert vector is not None
+    restored = unvec(vector, SHAPE)
+    assert restored is not None
+    assert restored.shape == (3, *SHAPE)
+    assert_array_equal(restored, FIELD)
+
+
+def test_complex_field_round_trip_preserves_shape_and_values() -> None:
+    vector = vec(COMPLEX_FIELD)
+    assert vector is not None
+    restored = unvec(vector, SHAPE)
+    assert restored is not None
+    assert restored.shape == (3, *SHAPE)
+    assert_allclose(restored, COMPLEX_FIELD)
+
+
+def test_unvec_with_two_components_round_trips_vector() -> None:
+    vector = numpy.arange(2 * numpy.prod(SHAPE), dtype=float)
+    field = unvec(vector, SHAPE, nvdim=2)
+    assert field is not None
+    assert field.shape == (2, *SHAPE)
+    assert_array_equal(vec(field), vector)
+
+
+def test_vec_accepts_arraylike_input() -> None:
+    arraylike = [[[1, 2], [3, 4]], [[5, 6], [7, 8]]]
+    assert_array_equal(vec(arraylike), numpy.ravel(arraylike, order='C'))
diff --git a/meanas/test/test_fdtd_energy.py b/meanas/test/test_fdtd_energy.py
new file mode 100644
index 0000000..2d15c69
--- /dev/null
+++ b/meanas/test/test_fdtd_energy.py
@@ -0,0 +1,98 @@
+import numpy
+from numpy.testing import assert_allclose
+
+from .. import fdtd
+from ..fdtd import energy as fdtd_energy
+
+
+SHAPE = (2, 2, 2)
+DT = 0.25
+UNIT_DXES = tuple(tuple(numpy.ones(length) for length in SHAPE) for _ in range(2))
+DXES = (
+    (
+        numpy.array([1.0, 1.5]),
+        numpy.array([0.75, 1.25]),
+        numpy.array([1.1, 0.9]),
+    ),
+    (
+        numpy.array([0.8, 1.2]),
+        numpy.array([1.4, 0.6]),
+        numpy.array([0.7, 1.3]),
+    ),
+)
+E0 = numpy.arange(3 * numpy.prod(SHAPE), dtype=float).reshape((3, *SHAPE), order='C')
+E1 = E0 + 0.5
+E2 = E0 + 1.0
+E3 = E0 + 1.5
+H0 = (numpy.arange(3 * numpy.prod(SHAPE), dtype=float).reshape((3, *SHAPE), order='C') + 2.0) / 3.0
+H1 = H0 + 0.25
+H2 = H0 + 0.5
+H3 = H0 + 0.75
+J0 = (E0 + 2.0) / 5.0
+EPSILON = 1.0 + E0 / 20.0
+MU = 1.5 + H0 / 10.0
+
+
+def test_poynting_default_spacing_matches_explicit_unit_spacing() -> None:
+    default_spacing = fdtd.poynting(E1, H1)
+    explicit_spacing = fdtd.poynting(E1, H1, dxes=UNIT_DXES)
+    assert_allclose(default_spacing, explicit_spacing)
+
+
+def test_poynting_divergence_matches_precomputed_poynting_vector() -> None:
+    s = fdtd.poynting(E2, H2, dxes=DXES)
+    from_fields = fdtd.poynting_divergence(e=E2, h=H2, dxes=DXES)
+    from_vector = fdtd.poynting_divergence(s=s)
+    assert_allclose(from_fields, from_vector)
+
+
+def test_delta_energy_h2e_matches_direct_dxmul_formula() -> None:
+    expected = fdtd_energy.dxmul(
+        E2 * (E2 - E0) / DT,
+        H1 * (H3 - H1) / DT,
+        EPSILON,
+        MU,
+        DXES,
+    )
+    actual = fdtd.delta_energy_h2e(
+        dt=DT,
+        e0=E0,
+        h1=H1,
+        e2=E2,
+        h3=H3,
+        epsilon=EPSILON,
+        mu=MU,
+        dxes=DXES,
+    )
+    assert_allclose(actual, expected)
+
+
+def test_delta_energy_e2h_matches_direct_dxmul_formula() -> None:
+    expected = fdtd_energy.dxmul(
+        E1 * (E3 - E1) / DT,
+        H2 * (H2 - H0) / DT,
+        EPSILON,
+        MU,
+        DXES,
+    )
+    actual = fdtd_energy.delta_energy_e2h(
+        dt=DT,
+        h0=H0,
+        e1=E1,
+        h2=H2,
+        e3=E3,
+        epsilon=EPSILON,
+        mu=MU,
+        dxes=DXES,
+    )
+    assert_allclose(actual, expected)
+
+
+def test_delta_energy_j_defaults_to_unit_cell_volume() -> None:
+    expected = (J0 * E1).sum(axis=0)
+    assert_allclose(fdtd.delta_energy_j(j0=J0, e1=E1), expected)
+
+
+def test_dxmul_defaults_to_unit_materials_and_spacing() -> None:
+    expected = E1.sum(axis=0) + H1.sum(axis=0)
+    assert_allclose(fdtd_energy.dxmul(E1, H1), expected)
diff --git a/meanas/test/test_regressions.py b/meanas/test/test_regressions.py
new file mode 100644
index 0000000..8231880
--- /dev/null
+++ b/meanas/test/test_regressions.py
@@ -0,0 +1,238 @@
+import numpy
+import pytest       # type: ignore
+from numpy.testing import assert_allclose
+from scipy import sparse
+
+from ..eigensolvers import power_iteration, rayleigh_quotient_iteration, signed_eigensolve
+from ..fdfd import eme, farfield
+from ..fdtd.boundaries import conducting_boundary
+from ..fdtd.misc import gaussian_beam, gaussian_packet, ricker_pulse
+from ..fdtd.pml import cpml_params, updates_with_cpml
+
+
+def _axis_index(axis: int, index: int) -> tuple[slice | int, ...]:
+    coords: list[slice | int] = [slice(None), slice(None), slice(None)]
+    coords[axis] = index
+    return tuple(coords)
+
+
+@pytest.mark.parametrize('one_sided', [False, True])
+def test_gaussian_packet_accepts_array_input(one_sided: bool) -> None:
+    dt = 0.01
+    source, delay = gaussian_packet(1.55, 0.1, dt, one_sided=one_sided)
+    steps = numpy.array([0, int(numpy.ceil(delay / dt)) + 5])
+    envelope, cc, ss = source(steps)
+
+    assert envelope.shape == (2,)
+    assert numpy.isfinite(envelope).all()
+    assert numpy.isfinite(cc).all()
+    assert numpy.isfinite(ss).all()
+    if one_sided:
+        assert envelope[-1] == pytest.approx(1.0)
+
+
+def test_ricker_pulse_returns_finite_values() -> None:
+    source, delay = ricker_pulse(1.55, 0.01)
+    envelope, cc, ss = source(numpy.array([0, 1, 2]))
+
+    assert numpy.isfinite(delay)
+    assert numpy.isfinite(envelope).all()
+    assert numpy.isfinite(cc).all()
+    assert numpy.isfinite(ss).all()
+
+
+def test_gaussian_beam_centered_grid_is_finite_and_normalized() -> None:
+    beam = gaussian_beam(
+        xyz=[numpy.linspace(-1, 1, 3), numpy.linspace(-1, 1, 3), numpy.linspace(-1, 1, 3)],
+        center=[0, 0, 0],
+        waist_radius=1.0,
+        wl=1.55,
+        )
+
+    row = beam[:, :, beam.shape[2] // 2]
+    assert numpy.isfinite(beam).all()
+    assert numpy.linalg.norm(row) == pytest.approx(1.0)
+
+
+@pytest.mark.parametrize('direction', [0, 1, 2])
+@pytest.mark.parametrize('polarity', [-1, 1])
+def test_conducting_boundary_updates_expected_faces(direction: int, polarity: int) -> None:
+    e = numpy.arange(3 * 4 * 4 * 4, dtype=float).reshape(3, 4, 4, 4)
+    h = e.copy()
+    e0 = e.copy()
+    h0 = h.copy()
+
+    update_e, update_h = conducting_boundary(direction, polarity)
+    update_e(e)
+    update_h(h)
+
+    dirs = [0, 1, 2]
+    dirs.remove(direction)
+    u, v = dirs
+
+    if polarity < 0:
+        boundary = _axis_index(direction, 0)
+        shifted1 = _axis_index(direction, 1)
+
+        assert_allclose(e[direction][boundary], 0)
+        assert_allclose(e[u][boundary], e0[u][shifted1])
+        assert_allclose(e[v][boundary], e0[v][shifted1])
+        assert_allclose(h[direction][boundary], h0[direction][shifted1])
+        assert_allclose(h[u][boundary], 0)
+        assert_allclose(h[v][boundary], 0)
+    else:
+        boundary = _axis_index(direction, -1)
+        shifted1 = _axis_index(direction, -2)
+        shifted2 = _axis_index(direction, -3)
+
+        assert_allclose(e[direction][boundary], -e0[direction][shifted2])
+        assert_allclose(e[direction][shifted1], 0)
+        assert_allclose(e[u][boundary], e0[u][shifted1])
+        assert_allclose(e[v][boundary], e0[v][shifted1])
+        assert_allclose(h[direction][boundary], h0[direction][shifted1])
+        assert_allclose(h[u][boundary], -h0[u][shifted2])
+        assert_allclose(h[u][shifted1], 0)
+        assert_allclose(h[v][boundary], -h0[v][shifted2])
+        assert_allclose(h[v][shifted1], 0)
+
+
+@pytest.mark.parametrize(
+    ('direction', 'polarity'),
+    [(-1, 1), (3, 1), (0, 0)],
+    )
+def test_conducting_boundary_rejects_invalid_arguments(direction: int, polarity: int) -> None:
+    with pytest.raises(Exception):
+        conducting_boundary(direction, polarity)
+
+
+def test_get_abcd_returns_sparse_block_matrix(monkeypatch: pytest.MonkeyPatch) -> None:
+    t = numpy.array([[1.0, 0.0], [0.0, 2.0]])
+    r = numpy.array([[0.0, 0.1], [0.2, 0.0]])
+    monkeypatch.setattr(eme, 'get_tr', lambda *args, **kwargs: (t, r))
+
+    abcd = eme.get_abcd(None, None, None, None)
+
+    assert sparse.issparse(abcd)
+    assert abcd.shape == (4, 4)
+    assert numpy.isfinite(abcd.toarray()).all()
+
+
+def test_get_s_force_reciprocal_symmetrizes_output(monkeypatch: pytest.MonkeyPatch) -> None:
+    left = object()
+    right = object()
+
+    def fake_get_tr(_eh_l, wavenumbers_l, _eh_r, _wavenumbers_r, **kwargs):
+        if wavenumbers_l is left:
+            return (
+                numpy.array([[1.0, 2.0], [0.5, 1.0]]),
+                numpy.array([[0.0, 1.0], [2.0, 0.0]]),
+                )
+        return (
+            numpy.array([[1.0, -1.0], [0.0, 1.0]]),
+            numpy.array([[0.0, 0.5], [1.5, 0.0]]),
+            )
+
+    monkeypatch.setattr(eme, 'get_tr', fake_get_tr)
+    ss = eme.get_s(None, left, None, right, force_reciprocal=True)
+
+    assert_allclose(ss, ss.T)
+
+
+def test_farfield_roundtrip_supports_rectangular_arrays() -> None:
+    e_near = [numpy.zeros((4, 8), dtype=complex), numpy.zeros((4, 8), dtype=complex)]
+    h_near = [numpy.zeros((4, 8), dtype=complex), numpy.zeros((4, 8), dtype=complex)]
+    e_near[0][1, 3] = 1.0 + 0.25j
+    h_near[1][2, 5] = -0.5j
+
+    ff = farfield.near_to_farfield(e_near, h_near, dx=0.2, dy=0.3, padded_size=(4, 8))
+    restored = farfield.far_to_nearfield(
+        [field.copy() for field in ff['E']],
+        [field.copy() for field in ff['H']],
+        ff['dkx'],
+        ff['dky'],
+        padded_size=(4, 8),
+        )
+
+    assert isinstance(ff['dkx'], float)
+    assert isinstance(ff['dky'], float)
+    assert ff['E'][0].shape == (4, 8)
+    assert restored['E'][0].shape == (4, 8)
+    assert restored['H'][0].shape == (4, 8)
+    assert restored['dx'] == pytest.approx(0.2)
+    assert restored['dy'] == pytest.approx(0.3)
+    assert numpy.isfinite(restored['E'][0]).all()
+    assert numpy.isfinite(restored['H'][0]).all()
+
+
+@pytest.mark.parametrize(
+    ('axis', 'polarity', 'thickness', 'epsilon_eff'),
+    [(3, 1, 4, 1.0), (0, 0, 4, 1.0), (0, 1, 2, 1.0), (0, 1, 4, 0.0)],
+    )
+def test_cpml_params_reject_invalid_arguments(axis: int, polarity: int, thickness: int, epsilon_eff: float) -> None:
+    with pytest.raises(Exception):
+        cpml_params(axis=axis, polarity=polarity, dt=0.1, thickness=thickness, epsilon_eff=epsilon_eff)
+
+
+def test_cpml_params_shapes_and_region() -> None:
+    params = cpml_params(axis=1, polarity=1, dt=0.1, thickness=3)
+    p0e, p1e, p2e = params['param_e']
+    p0h, p1h, p2h = params['param_h']
+
+    assert p0e.shape == (1, 3, 1)
+    assert p1e.shape == (1, 3, 1)
+    assert p2e.shape == (1, 3, 1)
+    assert p0h.shape == (1, 3, 1)
+    assert p1h.shape == (1, 3, 1)
+    assert p2h.shape == (1, 3, 1)
+    assert params['region'][1] == slice(-3, None)
+
+
+def test_updates_with_cpml_keeps_zero_fields_zero() -> None:
+    shape = (3, 4, 4, 4)
+    epsilon = numpy.ones(shape, dtype=float)
+    e = numpy.zeros(shape, dtype=float)
+    h = numpy.zeros(shape, dtype=float)
+    dxes = [[numpy.ones(4), numpy.ones(4), numpy.ones(4)] for _ in range(2)]
+    params = [[None, None] for _ in range(3)]
+    params[0][0] = cpml_params(axis=0, polarity=-1, dt=0.1, thickness=3)
+
+    update_e, update_h = updates_with_cpml(params, dt=0.1, dxes=dxes, epsilon=epsilon)
+    update_e(e, h, epsilon)
+    update_h(e, h)
+
+    assert not e.any()
+    assert not h.any()
+
+
+def test_power_iteration_finds_dominant_mode() -> None:
+    operator = sparse.diags([5.0, 3.0, 1.0, -2.0]).tocsr()
+    eigval, eigvec = power_iteration(operator, guess_vector=numpy.ones(4, dtype=complex), iterations=20)
+
+    assert eigval == pytest.approx(5.0, rel=1e-6)
+    assert abs(eigvec[0]) > abs(eigvec[1])
+
+
+def test_rayleigh_quotient_iteration_refines_known_mode() -> None:
+    operator = sparse.diags([5.0, 3.0, 1.0, -2.0]).tocsr()
+
+    def solver(matrix: sparse.spmatrix, rhs: numpy.ndarray) -> numpy.ndarray:
+        return numpy.linalg.lstsq(matrix.toarray(), rhs, rcond=None)[0]
+
+    eigval, eigvec = rayleigh_quotient_iteration(
+        operator,
+        numpy.array([1.0, 0.1, 0.0, 0.0], dtype=complex),
+        iterations=8,
+        solver=solver,
+        )
+
+    residual = numpy.linalg.norm(operator @ eigvec - eigval * eigvec)
+    assert eigval == pytest.approx(3.0, rel=1e-6)
+    assert residual < 1e-8
+
+
+def test_signed_eigensolve_returns_largest_positive_modes() -> None:
+    operator = sparse.diags([5.0, 3.0, 1.0, -2.0]).tocsr()
+    eigvals, eigvecs = signed_eigensolve(operator, how_many=2)
+
+    assert_allclose(eigvals, [3.0, 5.0], atol=1e-6)
+    assert eigvecs.shape == (4, 2)