diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
new file mode 100644
index 0000000..851c3a7
--- /dev/null
+++ b/.github/workflows/tests.yml
@@ -0,0 +1,26 @@
+name: Tests
+
+on:
+  push:
+    branches: [master]
+  pull_request:
+
+jobs:
+  test:
+    runs-on: ubuntu-latest
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Set up Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: "3.11"
+
+      - name: Install dependencies
+        run: |
+          pip install lalsuite pycbc pytest black
+          pip install -e .
+
+      - name: Run tests
+        run: pytest tests/ -v
diff --git a/esigmapy/blend.py b/esigmapy/blend.py
index d2be5ce..7f42c85 100644
--- a/esigmapy/blend.py
+++ b/esigmapy/blend.py
@@ -1,6 +1,6 @@
 # Copyright (C) 2023 Kartikey Sharma, Prayush Kumar
 #
-"""Master function to hybridise any complex timeseries using the 'frequency' as a user input, 
+"""Master function to hybridise any complex timeseries using the 'frequency' as a user input,
 specifically to be used for gravitational waveform hybridisation, fine-tuned for a single mode.
 """
 
@@ -100,17 +100,13 @@ def align_in_phase(
     align_merger_to_inspiral=True,
 ):
     if len(inspiral) == 0:
-        raise IOError(
-            f"""You passed an inspiral waveform of zero length, to align
-                      with merger-ringdown!"""
-        )
+        raise IOError(f"""You passed an inspiral waveform of zero length, to align
+                      with merger-ringdown!""")
     if (t2_index_insp + 1 - t1_index_insp) < 1:
-        raise IOError(
-            f"""You have passed a very narrow window for the inspiral
+        raise IOError(f"""You have passed a very narrow window for the inspiral
                       waveform's hybridization. As per your input, the inspiral
                       waveform from index {t1_index_insp} to {t2_index_insp + 1}
-                      should be used to attach merger-ringdown"""
-        )
+                      should be used to attach merger-ringdown""")
 
     # Function alignes the two waveforms using the phase, optimised over the attachment region
     # m from l,m mode
@@ -237,10 +233,8 @@ def blend_modes(
         Set this to True to enable logging output.
     """
     if frq_width <= 0:
-        raise IOError(
-            f"""You are trying to blend over a frequency window of
-            negative length (= {frq_width}Hz). Fix this."""
-        )
+        raise IOError(f"""You are trying to blend over a frequency window of
+            negative length (= {frq_width}Hz). Fix this.""")
     if blend_using_avg_orbital_frequency:
         if len(inspiral_orbital_frequency) != len(inspiral_modes[mode_to_align_by]):
             raise IOError(
@@ -327,14 +321,12 @@ def blend_modes(
             inspiral_orbital_frequency, (frq_attach + frq_width / 2) / em
         )
         if verbose > 1:
-            print(
-                f"""Hybridizing using orbital frequency. Frequency
+            print(f"""Hybridizing using orbital frequency. Frequency
                   {frq_attach + frq_width / 2}Hz found at {t2_index_insp}.
                   The same frequency would have been found at index
                   {find_last_value_location_in_series(frq_insp[(el, em)],
                   frq_attach + frq_width / 2)} of mode frequency evolution.
-                  """
-            )
+                  """)
     else:
         t2_index_insp = find_last_value_location_in_series(
             frq_insp[(el, em)], frq_attach + frq_width / 2
@@ -344,17 +336,13 @@ def blend_modes(
     t1_index_insp = t2_index_insp - (t2_index_mr - t1_index_mr)
 
     if verbose > 4:
-        print(
-            f"""In the merger-ringdown waveform, the hybridization frequency window
+        print(f"""In the merger-ringdown waveform, the hybridization frequency window
               [{frq_attach - frq_width / 2}, {frq_attach + frq_width / 2}]
               was found at indices [{t1_index_mr}, {t2_index_mr}]
-              """
-        )
-        print(
-            f"""In the inspiral waveform, the same window is to be found at
+              """)
+        print(f"""In the inspiral waveform, the same window is to be found at
               indices [{t1_index_insp}, {t2_index_insp}.]
-              """
-        )
+              """)
     """ 
         Theoretically, we NEED a timeshift to align the waveforms in frequency. 
         Instead of shifting one of the two waveforms for alignment, we are defining
diff --git a/esigmapy/generator.py b/esigmapy/generator.py
index cde63ee..c8332db 100644
--- a/esigmapy/generator.py
+++ b/esigmapy/generator.py
@@ -269,7 +269,7 @@ def get_inspiral_esigma_modes(
 
     if f_ref < f_lower:
         ref_idx = np.searchsorted(
-            (retval[1].data.data ** 1.5) / ((mass1 + mass2) * lal.MTSUN_SI * np.pi),
+            (retval[1].data.data**1.5) / ((mass1 + mass2) * lal.MTSUN_SI * np.pi),
             f_lower,
         )
         new_len = len(retval[0].data.data) - ref_idx
@@ -535,10 +535,8 @@ def _get_transition_frequency_window(
     # frequency, towards the end of waveform
 
     if verbose > 4:
-        print(
-            f"""Transition frequency found at index {transition_idx}
-            in the orbital frequency array of length {len(orbital_freq)}"""
-        )
+        print(f"""Transition frequency found at index {transition_idx}
+            in the orbital frequency array of length {len(orbital_freq)}""")
 
     if blend_using_avg_orbital_frequency:
         # We integrate `orbital_freq` to obtain `orbital_phase`
@@ -640,11 +638,9 @@ def _get_transition_frequency_window(
         )
 
         if verbose > 4:
-            print(
-                f"""The transition is to happen in a window from
+            print(f"""The transition is to happen in a window from
                 frequencies: [{orbital_freq[window_start_idx]}, {orbital_freq[transition_idx]}]Hz,
-                between indices: [{window_start_idx}, {transition_idx}]"""
-            )
+                between indices: [{window_start_idx}, {transition_idx}]""")
     else:
         raise RuntimeError(
             """We should never reach here. Did you set a non-bool value for the
@@ -772,22 +768,16 @@ def get_imr_esigma_modes(
     """
     check_available_mr_approximants(merger_ringdown_approximant)
     if (mean_anomaly is None) and (coa_phase is None):
-        raise IOError(
-            f"""Please specify one of the phase angles, either of
-                      `mean_anomaly` or `coa_phase`."""
-        )
+        raise IOError(f"""Please specify one of the phase angles, either of
+                      `mean_anomaly` or `coa_phase`.""")
     if blend_aligning_merger_to_inspiral and (mean_anomaly is None):
-        raise IOError(
-            f"""If you want to attach ESIGMA inspiral to merger, by
+        raise IOError(f"""If you want to attach ESIGMA inspiral to merger, by
                       phase shifting merger to inspiral, please specify the
-                      phase angle `mean_anomaly`"""
-        )
+                      phase angle `mean_anomaly`""")
     if (not blend_aligning_merger_to_inspiral) and (coa_phase is None):
-        raise IOError(
-            f"""If you want to attach ESIGMA inspiral to merger, by
+        raise IOError(f"""If you want to attach ESIGMA inspiral to merger, by
                       phase shifting inspiral to merger, please specify the
-                      phase angle `coa_phase`"""
-        )
+                      phase angle `coa_phase`""")
     if mean_anomaly is None:
         mean_anomaly = 0
     if coa_phase is None:
@@ -807,13 +797,11 @@ def get_imr_esigma_modes(
             set(available_inspiral_orbital_params)
         )
         if return_orbital_params_user != set(return_orbital_params):
-            print(
-                f"""Warning: You requested the following list of orbital
+            print(f"""Warning: You requested the following list of orbital
 parameters to be returned: {return_orbital_params}, but we reduce it to
 {return_orbital_params_user} as we only have the evolution of the following
 parameters available with us: {available_inspiral_orbital_params}.
-                  """
-            )
+                  """)
     elif not return_orbital_params:
         return_orbital_params = []
         return_orbital_params_user = False
@@ -871,26 +859,21 @@ def get_imr_esigma_modes(
         raise RuntimeError(
             f"""The inspiral modes do not contain the primary
 desired {mode_to_align_by} multipole. It currently holds only the following:
-{modes_inspiral_numpy.keys()}"""
-        )
+{modes_inspiral_numpy.keys()}""")
 
     orbital_eccentricity = retval[-2]["e"]
     # Throw error if eccentricity at the end of inspiral is definitely unsafe
     if orbital_eccentricity[-1] > ECCENTRICITY_LEVEL_ISCO_ERROR:
-        raise IOError(
-            f"""ERROR: You entered a very large initial eccentricity
+        raise IOError(f"""ERROR: You entered a very large initial eccentricity
 {eccentricity}. The orbital eccentricity at the end of inspiral was
 {orbital_eccentricity[-1]}. The merger-ringdown attachment with a
-quasicircular will be questionable."""
-        )
+quasicircular will be questionable.""")
     # Warn user if eccentricity at the end of inspiral is potentially unsafe
     if orbital_eccentricity[-1] > ECCENTRICITY_LEVEL_ISCO_WARNING and verbose:
-        print(
-            f"""WARNING: You entered a very large initial eccentricity
+        print(f"""WARNING: You entered a very large initial eccentricity
 {eccentricity}. The orbital eccentricity at the end of inspiral was
 {orbital_eccentricity[-1]}. The merger-ringdown attachment with a quasicircular
-model might be affected."""
-        )
+model might be affected.""")
 
     if (f_window_mr_transition is None) or failsafe or (verbose > 1):
         if blend_using_avg_orbital_frequency:
@@ -914,14 +897,12 @@ def get_imr_esigma_modes(
             modes_inspiral_numpy[mode_to_align_by]
         )
         mode_frequency = esigmapy.blend.compute_frequency(mode_phase, delta_t)
-        print(
-            f"""DEBUG: Orbital freq at end of inspiral is {orbital_frequency[-1]}Hz,
+        print(f"""DEBUG: Orbital freq at end of inspiral is {orbital_frequency[-1]}Hz,
 mode-{mode_to_align_by} freq at the end of inspiral is {mode_frequency[-1]}Hz, max and min
 mode-{mode_to_align_by} frequencies are {np.max(mode_frequency)}Hz and
 {np.min(mode_frequency)}Hz, and the transition frequency (of {mode_to_align_by}-mode)
 requested is {f_mr_transition}Hz, which should be less than the maximum freq of
-{mode_to_align_by}-mode: {mode_frequency.max()}Hz."""
-        )
+{mode_to_align_by}-mode: {mode_frequency.max()}Hz.""")
         return (
             modes_inspiral_numpy,
             mode_phase,
@@ -940,12 +921,10 @@ def get_imr_esigma_modes(
         mode_frequency = esigmapy.blend.compute_frequency(mode_phase, delta_t)
         if mode_frequency.max() < f_mr_transition:
             if verbose:
-                print(
-                    f"""FAILSAFE: Maximum orbital freq during inspiral is
+                print(f"""FAILSAFE: Maximum orbital freq during inspiral is
 {orbital_frequency.max()}Hz, and max frequency of {mode_to_align_by}-mode is
 {mode_frequency.max()}Hz, so we are resetting transition frequency from
-{f_mr_transition}Hz to {mode_frequency.max()}Hz."""
-                )
+{f_mr_transition}Hz to {mode_frequency.max()}Hz.""")
             f_mr_transition = mode_frequency.max()
 
     # If the user does not provide the width of hybridization window (
@@ -1026,12 +1005,10 @@ def get_imr_esigma_modes(
             verbose=verbose,
         )
     except Exception as exc:
-        print(
-            f"""Inspiral + MergerRingdown attachment failed. Its very likely
+        print(f"""Inspiral + MergerRingdown attachment failed. Its very likely
 that you entered a very large initial eccentricity {eccentricity}. The orbital
 eccentricity at the end of inspiral was {orbital_eccentricity[-1]}
-              """
-        )
+              """)
         raise exc
     modes_imr_numpy = retval[0]
 
diff --git a/esigmapy/legacy.py b/esigmapy/legacy.py
index c69e167..d07363b 100644
--- a/esigmapy/legacy.py
+++ b/esigmapy/legacy.py
@@ -1,7 +1,6 @@
 # Copyright (C) 2023 Prayush Kumar
 #
-"""Legacy code to calibrate transition/attachment of inspiral and merger-ringdown
-"""
+"""Legacy code to calibrate transition/attachment of inspiral and merger-ringdown"""
 
 
 class FitMOmegaIMRAttachmentNonSpinning:
diff --git a/esigmapy/surrogate/generator.py b/esigmapy/surrogate/generator.py
index aee8d63..7754c29 100644
--- a/esigmapy/surrogate/generator.py
+++ b/esigmapy/surrogate/generator.py
@@ -422,13 +422,11 @@ def get_imr_esigmasur_mode(
             set(available_inspiral_orbital_params)
         )
         if return_orbital_params_user != set(return_orbital_params):
-            print(
-                f"""Warning: You requested the following list of orbital
+            print(f"""Warning: You requested the following list of orbital
 parameters to be returned: {return_orbital_params}, but we reduce it to
 {return_orbital_params_user} as we only have the evolution of the following 
 parameters available with us: {available_inspiral_orbital_params}.
-                  """
-            )
+                  """)
     elif not return_orbital_params:
         return_orbital_params = []
         return_orbital_params_user = False
@@ -475,29 +473,23 @@ def get_imr_esigmasur_mode(
     modes_inspiral_numpy = retval[-1]
 
     if mode_to_align_by not in modes_inspiral_numpy:
-        raise RuntimeError(
-            f"""The inspiral modes do not contain the primary 
+        raise RuntimeError(f"""The inspiral modes do not contain the primary 
 desired {mode_to_align_by} multipole. It currently holds only the following:
-{modes_inspiral_numpy.keys()}"""
-        )
+{modes_inspiral_numpy.keys()}""")
 
     orbital_eccentricity = retval[-2]["e"]
     # Throw error if eccentricity at the end of inspiral is definitely unsafe
     if orbital_eccentricity[-1] > ECCENTRICITY_LEVEL_ISCO_ERROR:
-        raise IOError(
-            f"""ERROR: You entered a very large reference eccentricity
+        raise IOError(f"""ERROR: You entered a very large reference eccentricity
 {reference_eccentricity}. The orbital eccentricity at the end of inspiral was
 {orbital_eccentricity[-1]}. The merger-ringdown attachment with a
-quasicircular will be questionable."""
-        )
+quasicircular will be questionable.""")
     # Warn user if eccentricity at the end of inspiral is potentially unsafe
     if orbital_eccentricity[-1] > ECCENTRICITY_LEVEL_ISCO_WARNING and verbose:
-        print(
-            f"""WARNING: You entered a very large reference eccentricity
+        print(f"""WARNING: You entered a very large reference eccentricity
 {reference_eccentricity}. The orbital eccentricity at the end of inspiral was
 {orbital_eccentricity[-1]}. The merger-ringdown attachment with a quasicircular
-model might be affected."""
-        )
+model might be affected.""")
 
     if (f_window_mr_transition is None) or failsafe or (verbose > 1):
         if blend_using_avg_orbital_frequency:
@@ -521,14 +513,12 @@ def get_imr_esigmasur_mode(
             modes_inspiral_numpy[mode_to_align_by]
         )
         mode_frequency = esigmapy.blend.compute_frequency(mode_phase, delta_t)
-        print(
-            f"""DEBUG: Orbital freq at end of inspiral is {orbital_frequency[-1]}Hz,
+        print(f"""DEBUG: Orbital freq at end of inspiral is {orbital_frequency[-1]}Hz,
 mode-{mode_to_align_by} freq at the end of inspiral is {mode_frequency[-1]}Hz, max and min
 mode-{mode_to_align_by} frequencies are {np.max(mode_frequency)}Hz and
 {np.min(mode_frequency)}Hz, and the transition frequency (of {mode_to_align_by}-mode)
 requested is {f_mr_transition}Hz, which should be less than the maximum freq of
-{mode_to_align_by}-mode: {mode_frequency.max()}Hz."""
-        )
+{mode_to_align_by}-mode: {mode_frequency.max()}Hz.""")
         return (
             modes_inspiral_numpy,
             mode_phase,
@@ -547,12 +537,10 @@ def get_imr_esigmasur_mode(
         mode_frequency = esigmapy.blend.compute_frequency(mode_phase, delta_t)
         if mode_frequency.max() < f_mr_transition:
             if verbose:
-                print(
-                    f"""FAILSAFE: Maximum orbital freq during inspiral is
+                print(f"""FAILSAFE: Maximum orbital freq during inspiral is
 {orbital_frequency.max()}Hz, and max frequency of {mode_to_align_by}-mode is
 {mode_frequency.max()}Hz, so we are resetting transition frequency from
-{f_mr_transition}Hz to {mode_frequency.max()}Hz."""
-                )
+{f_mr_transition}Hz to {mode_frequency.max()}Hz.""")
             f_mr_transition = mode_frequency.max()
 
     # If the user does not provide the width of hybridization window (
@@ -633,12 +621,10 @@ def get_imr_esigmasur_mode(
             verbose=verbose,
         )
     except Exception as exc:
-        print(
-            f"""Inspiral + MergerRingdown attachment failed. Its very likely
+        print(f"""Inspiral + MergerRingdown attachment failed. Its very likely
 that you entered a very large reference eccentricity {reference_eccentricity}. The orbital
 eccentricity at the end of inspiral was {orbital_eccentricity[-1]}
-              """
-        )
+              """)
         raise exc
     modes_imr_numpy = retval[0]
 
diff --git a/esigmapy/surrogate/surrogate.py b/esigmapy/surrogate/surrogate.py
index 4f99a2c..8c94311 100644
--- a/esigmapy/surrogate/surrogate.py
+++ b/esigmapy/surrogate/surrogate.py
@@ -170,7 +170,6 @@ def __call__(
         return_amp_phase_only=False,
         return_orbital_variables=False,
     ):
-
         if delta_t is None and times is None:
             raise ValueError("Either delta_t or times must be provided.")
 
@@ -290,9 +289,12 @@ def check_param_range(self, q, e_ref, l_ref, override=False):
                 )
 
     def load_sur_metadata(self):
-        self.sur_total_mass, self.t_ref, self.t_grid_sur, self.l_grid_sur = (
-            self.get_metadata(["M", "t_ref", "t_grid_sur", "l_grid_sur"])
-        )
+        (
+            self.sur_total_mass,
+            self.t_ref,
+            self.t_grid_sur,
+            self.l_grid_sur,
+        ) = self.get_metadata(["M", "t_ref", "t_grid_sur", "l_grid_sur"])
 
     def load_ei_indices(self):
         filename_ei_indices = os.path.join(self.sur_dir, f"ei_indices.npz")
@@ -364,7 +366,6 @@ def __call__(
         remove_start_phase=True,
         return_orbital_variables=False,
     ):
-
         if delta_t is None and times is None:
             raise ValueError("Either delta_t or times must be provided.")
 
diff --git a/tests/__init__.py b/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 0000000..0bcc4d5
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,6 @@
+import logging
+
+import esigmapy.legacy as _legacy_module
+
+# legacy.py uses logging.info without importing logging; patch it so tests work
+_legacy_module.logging = logging
diff --git a/tests/test_blend.py b/tests/test_blend.py
new file mode 100644
index 0000000..b983ff1
--- /dev/null
+++ b/tests/test_blend.py
@@ -0,0 +1,466 @@
+import numpy as np
+import pytest
+
+from esigmapy.blend import (
+    align_in_phase,
+    blend_modes,
+    blend_series,
+    compute_amplitude,
+    compute_frequency,
+    compute_phase,
+    find_first_value_location_in_series,
+    find_last_value_location_in_series,
+    mismatch_discrete,
+)
+
+# ---------------------------------------------------------------------------
+# find_first_value_location_in_series
+# ---------------------------------------------------------------------------
+
+
+class TestFindFirstValue:
+    def test_below_min_raises(self):
+        with pytest.raises(Exception, match="[Ll]ower"):
+            find_first_value_location_in_series(np.array([1.0, 2.0, 3.0]), 0.5)
+
+    def test_above_max_raises(self):
+        with pytest.raises(Exception, match="[Hh]igher"):
+            find_first_value_location_in_series(np.array([1.0, 2.0, 3.0]), 3.5)
+
+    def test_returns_first_bracket_index(self):
+        # 1.5 lies between arr[1]=1 and arr[2]=2; nearest is equidistant → either valid
+        arr = np.array([0.0, 1.0, 2.0, 3.0])
+        idx = find_first_value_location_in_series(arr, 1.5)
+        assert idx in (1, 2)
+
+    def test_nearer_left_element_chosen(self):
+        # 1.2 is closer to arr[1]=1.0 than arr[2]=2.0
+        arr = np.array([0.0, 1.0, 2.0, 3.0])
+        idx = find_first_value_location_in_series(arr, 1.2)
+        assert idx == 1
+
+    def test_nearer_right_element_chosen(self):
+        # 1.8 is closer to arr[2]=2.0 than arr[1]=1.0
+        arr = np.array([0.0, 1.0, 2.0, 3.0])
+        idx = find_first_value_location_in_series(arr, 1.8)
+        assert idx == 2
+
+    def test_exact_match_at_start(self):
+        arr = np.array([1.0, 2.0, 3.0, 4.0])
+        idx = find_first_value_location_in_series(arr, 1.0)
+        assert idx == 0
+
+    def test_first_occurrence_returned_for_non_monotone(self):
+        # arr crosses 0.7 between idx 0→1 (first) and again between 2→3
+        arr = np.array([0.0, 1.0, 0.5, 1.0, 2.0])
+        idx = find_first_value_location_in_series(arr, 0.7)
+        # First bracket is [0→1]; nearest of 0.0 and 1.0 to 0.7 is 1.0 → idx 1
+        assert idx <= 1
+
+
+class TestFindLastValue:
+    def test_below_min_raises(self):
+        with pytest.raises(Exception):
+            find_last_value_location_in_series(np.array([1.0, 2.0, 3.0]), 0.5)
+
+    def test_above_max_raises(self):
+        with pytest.raises(Exception):
+            find_last_value_location_in_series(np.array([1.0, 2.0, 3.0]), 3.5)
+
+    def test_monotone_returns_index_at_least_as_large_as_first(self):
+        arr = np.linspace(0.0, 10.0, 100)
+        target = 5.0
+        idx_first = find_first_value_location_in_series(arr, target)
+        idx_last = find_last_value_location_in_series(arr, target)
+        assert idx_last >= idx_first
+
+    def test_returns_later_occurrence_for_non_monotone(self):
+        # Two crossings of 1.5 in arr: between idx 1-2 and between idx 3-4
+        arr = np.array([0.0, 1.0, 2.0, 1.0, 2.0, 3.0])
+        idx = find_last_value_location_in_series(arr, 1.5)
+        # The last crossing is in the 3-4 region; result should be >= 3
+        assert idx >= 3
+
+    def test_result_within_valid_range(self):
+        arr = np.linspace(1.0, 5.0, 50)
+        idx = find_last_value_location_in_series(arr, 3.0)
+        assert 0 <= idx < len(arr)
+
+
+# ---------------------------------------------------------------------------
+# compute_amplitude
+# ---------------------------------------------------------------------------
+
+
+class TestComputeAmplitude:
+    def test_known_magnitude(self):
+        # |3 + 4j| == 5
+        wave = np.full(10, 3.0 + 4.0j)
+        amp = compute_amplitude(wave)
+        assert np.allclose(amp, 5.0)
+
+    def test_real_valued_input(self):
+        wave = np.array([-3.0, -1.0, 0.0, 1.0, 3.0])
+        assert np.allclose(compute_amplitude(wave), np.abs(wave))
+
+    def test_zero_amplitude(self):
+        wave = np.zeros(5, dtype=complex)
+        assert np.allclose(compute_amplitude(wave), 0.0)
+
+    def test_output_nonnegative(self):
+        rng = np.random.default_rng(0)
+        wave = rng.standard_normal(50) + 1j * rng.standard_normal(50)
+        assert np.all(compute_amplitude(wave) >= 0)
+
+    def test_unit_circle(self):
+        # Points on the unit circle all have amplitude 1
+        phases = np.linspace(0, 2 * np.pi, 100, endpoint=False)
+        wave = np.exp(1j * phases)
+        assert np.allclose(compute_amplitude(wave), 1.0)
+
+
+# ---------------------------------------------------------------------------
+# compute_phase
+# ---------------------------------------------------------------------------
+
+
+class TestComputePhase:
+    def test_linear_phase_recovered(self):
+        # waveform = exp(-i * phi) → compute_phase returns phi (GW convention)
+        N, dt, f = 512, 1.0 / 4096, 80.0
+        t = np.arange(N) * dt
+        phi = 2 * np.pi * f * t
+        wave = np.exp(-1j * phi)
+        recovered = compute_phase(wave)
+        # Allow for unwrapping edge at boundaries; check interior
+        assert np.allclose(recovered[1:-1], phi[1:-1], atol=1e-10)
+
+    def test_constant_phase_is_flat(self):
+        # All samples at the same phase → output is a constant array
+        wave = np.exp(-1j * 0.7) * np.ones(50)
+        phase = compute_phase(wave)
+        assert np.allclose(phase, phase[0])
+
+    def test_opposite_sign_convention(self):
+        # exp(+i*phi) → compute_phase returns -phi (negative sign in formula)
+        N, dt, f = 256, 1.0 / 4096, 50.0
+        t = np.arange(N) * dt
+        phi = 2 * np.pi * f * t
+        wave = np.exp(1j * phi)
+        recovered = compute_phase(wave)
+        assert np.allclose(recovered[1:-1], -phi[1:-1], atol=1e-10)
+
+    def test_output_unwrapped(self):
+        # Phase must not have jumps larger than π between consecutive samples
+        N, dt, f = 1024, 1.0 / 4096, 200.0
+        t = np.arange(N) * dt
+        wave = np.exp(-1j * 2 * np.pi * f * t)
+        phase = compute_phase(wave)
+        assert np.all(np.abs(np.diff(phase)) < np.pi)
+
+
+# ---------------------------------------------------------------------------
+# compute_frequency
+# ---------------------------------------------------------------------------
+
+
+class TestComputeFrequency:
+    def test_constant_frequency_recovered(self):
+        N, dt, f0 = 512, 1.0 / 4096, 100.0
+        t = np.arange(N) * dt
+        phase = 2 * np.pi * f0 * t  # linear phase
+        freq = compute_frequency(phase, dt)
+        # np.gradient is exact for linear functions at all points
+        assert np.allclose(freq, f0, rtol=1e-10)
+
+    def test_negative_frequency_for_decreasing_phase(self):
+        N, dt, f0 = 256, 1.0 / 4096, 50.0
+        t = np.arange(N) * dt
+        phase = -2 * np.pi * f0 * t  # decreasing phase
+        freq = compute_frequency(phase, dt)
+        assert np.allclose(freq, -f0, rtol=1e-10)
+
+    def test_output_length_matches_input(self):
+        phase = np.linspace(0, 10 * np.pi, 200)
+        freq = compute_frequency(phase, 1.0 / 4096)
+        assert len(freq) == len(phase)
+
+    def test_roundtrip_phase_to_frequency(self):
+        # compute_frequency(compute_phase(exp(-i*2*pi*f*t)), dt) ≈ f
+        N, dt, f0 = 1024, 1.0 / 4096, 120.0
+        t = np.arange(N) * dt
+        wave = np.exp(-1j * 2 * np.pi * f0 * t)
+        freq = compute_frequency(compute_phase(wave), dt)
+        # Ignore first and last points where gradient uses one-sided differences
+        assert np.allclose(freq[1:-1], f0, rtol=1e-10)
+
+
+# ---------------------------------------------------------------------------
+# mismatch_discrete
+# ---------------------------------------------------------------------------
+
+
+class TestMismatchDiscrete:
+    @staticmethod
+    def _all_idx(n):
+        return np.arange(n)
+
+    def test_identical_waveforms_zero_mismatch(self):
+        w = np.exp(-1j * np.linspace(0, 4 * np.pi, 50))
+        idx = self._all_idx(len(w))
+        assert mismatch_discrete(w, w, idx, idx) == pytest.approx(0.0, abs=1e-14)
+
+    def test_quadrature_phase_gives_mismatch_one(self):
+        # w2 = i * w1  →  |w1 - w2|^2 = 2|w1|^2  →  mm = 0.5 * 2 = 1.0
+        w1 = np.ones(20, dtype=complex)
+        w2 = 1j * w1
+        idx = self._all_idx(len(w1))
+        assert mismatch_discrete(w1, w2, idx, idx) == pytest.approx(1.0, rel=1e-12)
+
+    def test_antiphase_gives_mismatch_two(self):
+        # w2 = -w1  →  |w1 - w2|^2 = 4|w1|^2  →  mm = 0.5 * 4 = 2.0
+        w1 = np.ones(20, dtype=complex)
+        w2 = -w1
+        idx = self._all_idx(len(w1))
+        assert mismatch_discrete(w1, w2, idx, idx) == pytest.approx(2.0, rel=1e-12)
+
+    def test_nonnegative(self):
+        rng = np.random.default_rng(1)
+        w1 = rng.standard_normal(30) + 1j * rng.standard_normal(30)
+        w2 = rng.standard_normal(30) + 1j * rng.standard_normal(30)
+        idx = self._all_idx(len(w1))
+        assert mismatch_discrete(w1, w2, idx, idx) >= 0
+
+    def test_subsample_indices(self):
+        # Using every other sample should still give 0 for identical waveforms
+        w = np.exp(-1j * np.linspace(0, 2 * np.pi, 40))
+        idx = np.arange(0, 40, 2)
+        assert mismatch_discrete(w, w, idx, idx) == pytest.approx(0.0, abs=1e-14)
+
+
+# ---------------------------------------------------------------------------
+# blend_series
+# ---------------------------------------------------------------------------
+
+
+class TestBlendSeries:
+    def test_left_edge_equals_x1(self):
+        x1 = np.ones(100)
+        x2 = np.zeros(100)
+        result = blend_series(x1, x2, 20, 40, 20, 40)
+        assert result[0] == pytest.approx(1.0)  # tau=0 at left edge
+
+    def test_right_edge_approaches_x2(self):
+        x1 = np.zeros(100)
+        x2 = np.ones(100)
+        result = blend_series(x1, x2, 20, 40, 20, 40)
+        # tau at last point = sin^2(pi/2 * (n-1)/n), very close to 1
+        assert result[-1] > 0.97
+
+    def test_output_length(self):
+        x1 = np.ones(100)
+        x2 = np.zeros(100)
+        t1, t2 = 10, 30
+        result = blend_series(x1, x2, t1, t2, t1, t2)
+        assert len(result) == t2 - t1
+
+    def test_monotone_blend_between_constants(self):
+        # Blending from 0→1 over a window must be monotone non-decreasing
+        x1 = np.zeros(100)
+        x2 = np.ones(100)
+        result = blend_series(x1, x2, 20, 60, 20, 60)
+        assert np.all(np.diff(result) >= -1e-14)
+
+    def test_unequal_window_lengths_raise(self):
+        x1 = np.ones(100)
+        x2 = np.ones(100)
+        with pytest.raises(AssertionError):
+            blend_series(x1, x2, 10, 30, 10, 25)  # insp window=20, mr window=15
+
+    def test_different_offsets_same_length(self):
+        # Windows can be at different positions but must have the same length
+        x1 = np.arange(100, dtype=float)
+        x2 = np.arange(100, dtype=float) * 2
+        result = blend_series(x1, x2, 10, 30, 50, 70)  # both windows length 20
+        assert len(result) == 20
+        assert result[0] == pytest.approx(x1[10])
+
+
+# ---------------------------------------------------------------------------
+# align_in_phase
+# ---------------------------------------------------------------------------
+
+
+class TestAlignInPhase:
+    @staticmethod
+    def _make_chirp(n=200, f0=50.0, dt=1.0 / 4096):
+        t = np.arange(n) * dt
+        return np.exp(-1j * 2 * np.pi * f0 * t)
+
+    @staticmethod
+    def _sample_idx(t1, t2, no_sp=8):
+        return np.linspace(t1, t2, no_sp, dtype=int) - t1
+
+    def test_empty_inspiral_raises(self):
+        with pytest.raises(IOError, match="[Zz]ero length"):
+            align_in_phase(
+                np.array([]),
+                np.ones(10, dtype=complex),
+                np.arange(5),
+                np.arange(5),
+                0,
+                4,
+                0,
+                4,
+            )
+
+    def test_narrow_window_raises(self):
+        with pytest.raises(IOError):
+            align_in_phase(
+                np.ones(10, dtype=complex),
+                np.ones(10, dtype=complex),
+                np.arange(1),
+                np.arange(1),
+                5,
+                4,
+                5,
+                4,  # t2 < t1: zero-width window
+            )
+
+    def test_identical_waveforms_zero_shift(self):
+        wave = self._make_chirp()
+        t1, t2 = 50, 100
+        idx = self._sample_idx(t1, t2)
+        _, _, shift = align_in_phase(wave, wave, idx, idx, t1, t2, t1, t2, m_mode=2)
+        # Optimal shift for identical waveforms is 0 (mod 2π/m)
+        assert np.cos(2 * shift) == pytest.approx(1.0, abs=1e-3)
+
+    def test_reduces_mismatch(self):
+        wave = self._make_chirp()
+        delta = np.pi / 4  # deliberate phase offset
+        mr = wave * np.exp(1j * delta)
+        t1, t2 = 50, 100
+        idx = self._sample_idx(t1, t2)
+
+        mm_before = mismatch_discrete(wave[t1 : t2 + 1], mr[t1 : t2 + 1], idx, idx)
+        insp_a, mr_a, _ = align_in_phase(
+            wave,
+            mr,
+            idx,
+            idx,
+            t1,
+            t2,
+            t1,
+            t2,
+            m_mode=2,
+            align_merger_to_inspiral=True,
+        )
+        mm_after = mismatch_discrete(insp_a[t1 : t2 + 1], mr_a[t1 : t2 + 1], idx, idx)
+        assert mm_after < mm_before
+
+    def test_align_merger_to_inspiral_order(self):
+        wave = self._make_chirp()
+        mr = wave * np.exp(1j * 0.3)
+        t1, t2 = 30, 80
+        idx = self._sample_idx(t1, t2)
+        insp_out, mr_out, _ = align_in_phase(
+            wave,
+            mr,
+            idx,
+            idx,
+            t1,
+            t2,
+            t1,
+            t2,
+            align_merger_to_inspiral=True,
+        )
+        # Inspiral is unchanged; merger-ringdown is shifted toward inspiral
+        assert np.allclose(insp_out, wave)
+        assert not np.allclose(mr_out, mr)  # mr was modified
+        assert np.allclose(mr_out, wave, atol=1e-6)  # and now matches wave
+
+    def test_align_inspiral_to_merger_order(self):
+        wave = self._make_chirp()
+        mr = wave * np.exp(1j * 0.3)
+        t1, t2 = 30, 80
+        idx = self._sample_idx(t1, t2)
+        insp_out, mr_out, _ = align_in_phase(
+            wave,
+            mr,
+            idx,
+            idx,
+            t1,
+            t2,
+            t1,
+            t2,
+            align_merger_to_inspiral=False,
+        )
+        # Merger-ringdown is unchanged; inspiral is shifted
+        assert np.allclose(mr_out, mr)
+        assert not np.allclose(insp_out, wave)
+
+
+# ---------------------------------------------------------------------------
+# blend_modes — input validation guards (pure Python, no LAL required)
+# ---------------------------------------------------------------------------
+
+
+class TestBlendModesValidation:
+    @staticmethod
+    def _modes(n=500, freq_hz=50.0, dt=1.0 / 4096):
+        t = np.arange(n) * dt
+        mode = np.exp(-2j * np.pi * freq_hz * t)
+        return {(2, 2): mode, (2, -2): np.conj(mode)}
+
+    def test_negative_frq_width_raises(self):
+        modes = self._modes()
+        with pytest.raises(IOError, match="negative"):
+            blend_modes(modes, modes, np.ones(500), 50.0, frq_width=-5.0)
+
+    def test_zero_frq_width_raises(self):
+        modes = self._modes()
+        with pytest.raises(IOError):
+            blend_modes(modes, modes, np.ones(500), 50.0, frq_width=0.0)
+
+    def test_mismatched_orbital_freq_length_raises(self):
+        modes = self._modes(n=500)
+        with pytest.raises(IOError):
+            blend_modes(
+                modes,
+                modes,
+                inspiral_orbital_frequency=np.ones(100),  # wrong length
+                frq_attach=50.0,
+                frq_width=5.0,
+                blend_using_avg_orbital_frequency=True,
+            )
+
+    def test_mode_missing_from_inspiral_raises(self):
+        insp = {(2, 2): np.ones(500, dtype=complex)}
+        mr = {(2, 2): np.ones(500, dtype=complex), (3, 3): np.ones(500, dtype=complex)}
+        with pytest.raises(IOError):
+            blend_modes(
+                insp,
+                mr,
+                np.ones(500),
+                50.0,
+                frq_width=5.0,
+                modes_to_blend=[(2, 2), (3, 3)],
+                include_conjugate_modes=False,
+            )
+
+    def test_mode_missing_from_mr_raises(self):
+        insp = {
+            (2, 2): np.ones(500, dtype=complex),
+            (3, 3): np.ones(500, dtype=complex),
+        }
+        mr = {(2, 2): np.ones(500, dtype=complex)}
+        with pytest.raises(IOError):
+            blend_modes(
+                insp,
+                mr,
+                np.ones(500),
+                50.0,
+                frq_width=5.0,
+                modes_to_blend=[(2, 2), (3, 3)],
+                include_conjugate_modes=False,
+            )
diff --git a/tests/test_formatting.py b/tests/test_formatting.py
new file mode 100644
index 0000000..932130a
--- /dev/null
+++ b/tests/test_formatting.py
@@ -0,0 +1,30 @@
+"""Enforce black formatting on the entire source tree.
+
+Running `black --check` here means any unformatted commit will be caught by
+`pytest` (and by CI) before it is merged.  To fix a failure, run:
+
+    black esigmapy/ tests/
+"""
+
+import subprocess
+import sys
+from pathlib import Path
+
+_ROOT = Path(__file__).parent.parent
+
+
+def test_black_formatting():
+    result = subprocess.run(
+        [
+            sys.executable,
+            "-m",
+            "black",
+            "--check",
+            str(_ROOT / "esigmapy"),
+            str(_ROOT / "tests"),
+        ],
+        capture_output=True,
+    )
+    assert result.returncode == 0, (
+        "black found unformatted files:\n" + result.stderr.decode()
+    )
diff --git a/tests/test_generator.py b/tests/test_generator.py
new file mode 100644
index 0000000..6582d36
--- /dev/null
+++ b/tests/test_generator.py
@@ -0,0 +1,363 @@
+import numpy as np
+import pytest
+from scipy import integrate
+from unittest.mock import MagicMock, patch
+
+import lal
+import pycbc.types as pt
+
+from esigmapy.generator import (
+    _get_window_start,
+    _get_transition_frequency_window,
+    get_imr_esigma_modes,
+    get_inspiral_esigma_modes,
+    eccentricity_at_reference_frequency,
+)
+
+# ---------------------------------------------------------------------------
+# Shared helpers for Tier 2 (LAL-mocked) tests
+# ---------------------------------------------------------------------------
+
+
+def _lal_ts(data):
+    """Minimal mock of a LAL REAL8TimeSeries."""
+    ts = MagicMock()
+    ts.data.data = np.asarray(data, dtype=float)
+    return ts
+
+
+def _lal_mode(n, freq_hz=50.0, dt=1.0 / 4096):
+    """Minimal mock of a LAL COMPLEX16TimeSeries for a GW mode."""
+    m = MagicMock()
+    t = np.arange(n) * dt
+    m.data.data = np.exp(-2j * np.pi * freq_hz * t)
+    return m
+
+
+def _fake_dynamics(n=64, dt=1.0 / 4096):
+    """Return a tuple of 8 mock LAL time series mimicking SimInspiralESIGMADynamics output."""
+    t_arr = np.arange(n, dtype=float) * dt
+    return (
+        _lal_ts(t_arr.copy()),  # t  (copied so in-place *= doesn't affect original)
+        _lal_ts(np.ones(n) * 0.1),  # x
+        _lal_ts(np.zeros(n)),  # e
+        _lal_ts(np.zeros(n)),  # l
+        _lal_ts(np.linspace(0, 5, n)),  # phi
+        _lal_ts(np.ones(n) * 100.0),  # phidot
+        _lal_ts(np.ones(n) * 10.0),  # r
+        _lal_ts(np.zeros(n)),  # rdot
+    )
+
+
+# ---------------------------------------------------------------------------
+# _get_window_start
+# ---------------------------------------------------------------------------
+
+
+class TestGetWindowStart:
+    """Tests for the private helper that integrates a frequency series and
+    returns the first index where the cumulative phase exceeds a threshold."""
+
+    _DT = 1.0 / 4096
+    _F0 = 100.0  # Hz — constant frequency used in several tests
+
+    def test_forward_meets_threshold(self):
+        freq = np.ones(500) * self._F0
+        idx = _get_window_start(freq, self._DT, 1.0, direction="forward")
+        assert idx is not None
+        assert abs(integrate.trapezoid(freq[: idx + 1], dx=self._DT)) >= 1.0
+
+    def test_forward_is_first_crossing(self):
+        # The index immediately before should NOT yet meet the threshold
+        freq = np.ones(500) * self._F0
+        idx = _get_window_start(freq, self._DT, 1.0, direction="forward")
+        assert abs(integrate.trapezoid(freq[:idx], dx=self._DT)) < 1.0
+
+    def test_forward_returns_none_when_unreachable(self):
+        # Total integral of [1 Hz × 5 samples × dt] is far below 1000 rad
+        freq = np.ones(5) * 1.0
+        assert _get_window_start(freq, 0.001, 1000.0, direction="forward") is None
+
+    def test_backward_meets_threshold(self):
+        freq = np.ones(500) * self._F0
+        idx = _get_window_start(freq, self._DT, 1.0, direction="backward")
+        assert idx is not None
+        assert abs(integrate.trapezoid(freq[idx:], dx=self._DT)) >= 1.0
+
+    def test_backward_is_last_crossing(self):
+        # The index one step to the right should NOT meet the threshold alone
+        freq = np.ones(500) * self._F0
+        idx = _get_window_start(freq, self._DT, 1.0, direction="backward")
+        assert abs(integrate.trapezoid(freq[idx + 1 :], dx=self._DT)) < 1.0
+
+    def test_backward_returns_none_when_unreachable(self):
+        freq = np.ones(5) * 1.0
+        assert _get_window_start(freq, 0.001, 1000.0, direction="backward") is None
+
+    def test_forward_index_in_valid_range(self):
+        freq = np.linspace(10.0, 200.0, 500)
+        idx = _get_window_start(freq, self._DT, 1.0, direction="forward")
+        assert idx is not None
+        assert 0 < idx < len(freq)
+
+
+# ---------------------------------------------------------------------------
+# _get_transition_frequency_window
+# ---------------------------------------------------------------------------
+
+
+class TestGetTransitionFrequencyWindow:
+    """Tests for the private helper that converts num_hyb_orbits into a
+    hybridization frequency window width."""
+
+    @staticmethod
+    def _setup(n=2000, dt=1.0 / 4096, f_start=10.0, f_end=200.0):
+        freq = np.linspace(f_start, f_end, n)
+        phase = np.cumsum(freq) * dt  # monotonically increasing
+        return freq, phase, dt
+
+    def test_end_mode_returns_positive_width(self):
+        freq, phase, dt = self._setup()
+        result = _get_transition_frequency_window(
+            phase,
+            freq,
+            dt,
+            f_mr_transition=freq[1000],
+            num_hyb_orbits=0.1,
+            keep_f_mr_transition_at_center=False,
+            blend_using_avg_orbital_frequency=False,
+            failsafe=True,
+        )
+        assert result > 0
+
+    def test_center_mode_returns_positive_width(self):
+        freq, phase, dt = self._setup()
+        result = _get_transition_frequency_window(
+            phase,
+            freq,
+            dt,
+            f_mr_transition=freq[1000],
+            num_hyb_orbits=0.1,
+            keep_f_mr_transition_at_center=True,
+            blend_using_avg_orbital_frequency=False,
+            failsafe=True,
+        )
+        assert result > 0
+
+    def test_avg_orbital_frequency_mode_returns_positive_width(self):
+        freq, phase, dt = self._setup()
+        result = _get_transition_frequency_window(
+            phase,
+            freq,
+            dt,
+            f_mr_transition=freq[1000],
+            num_hyb_orbits=0.1,
+            keep_f_mr_transition_at_center=False,
+            blend_using_avg_orbital_frequency=True,
+            failsafe=True,
+        )
+        assert result > 0
+
+    def test_more_orbits_wider_or_equal_window(self):
+        freq, phase, dt = self._setup()
+        f_tr = freq[1000]
+        w_narrow = _get_transition_frequency_window(
+            phase,
+            freq,
+            dt,
+            f_tr,
+            num_hyb_orbits=0.1,
+            keep_f_mr_transition_at_center=False,
+            blend_using_avg_orbital_frequency=False,
+            failsafe=True,
+        )
+        w_wide = _get_transition_frequency_window(
+            phase,
+            freq,
+            dt,
+            f_tr,
+            num_hyb_orbits=0.5,
+            keep_f_mr_transition_at_center=False,
+            blend_using_avg_orbital_frequency=False,
+            failsafe=True,
+        )
+        assert w_wide >= w_narrow
+
+
+# ---------------------------------------------------------------------------
+# get_imr_esigma_modes — input validation guards (fire before any LAL call)
+# ---------------------------------------------------------------------------
+
+
+class TestGetImrEsigmaModesValidation:
+    _BASE = dict(
+        mass1=20.0,
+        mass2=20.0,
+        f_lower=20.0,
+        delta_t=1.0 / 2048,
+        merger_ringdown_approximant="NRSur7dq4",
+    )
+
+    def test_invalid_approximant_raises_before_lal(self):
+        kwargs = {**self._BASE, "merger_ringdown_approximant": "IMRPhenomD"}
+        with pytest.raises(IOError):
+            get_imr_esigma_modes(**kwargs, mean_anomaly=0.0)
+
+    def test_both_phase_angles_none_raises(self):
+        with pytest.raises(IOError):
+            get_imr_esigma_modes(**self._BASE, mean_anomaly=None, coa_phase=None)
+
+    def test_align_merger_without_mean_anomaly_raises(self):
+        # blend_aligning_merger_to_inspiral=True (default) requires mean_anomaly
+        with pytest.raises(IOError):
+            get_imr_esigma_modes(
+                **self._BASE,
+                blend_aligning_merger_to_inspiral=True,
+                mean_anomaly=None,
+                coa_phase=0.0,
+            )
+
+    def test_align_inspiral_without_coa_phase_raises(self):
+        # blend_aligning_merger_to_inspiral=False requires coa_phase
+        with pytest.raises(IOError):
+            get_imr_esigma_modes(
+                **self._BASE,
+                blend_aligning_merger_to_inspiral=False,
+                mean_anomaly=0.0,
+                coa_phase=None,
+            )
+
+
+# ---------------------------------------------------------------------------
+# get_inspiral_esigma_modes — Tier 2 (LAL mocked)
+# ---------------------------------------------------------------------------
+
+
+class TestGetInspiralEsigmaModesWithMock:
+    """Tests for the Python orchestration logic in get_inspiral_esigma_modes.
+    lalsimulation is replaced by a MagicMock so no C library is exercised."""
+
+    @pytest.fixture
+    def ctx(self):
+        n, dt = 64, 1.0 / 4096
+        dyn = _fake_dynamics(n, dt)
+        mode = _lal_mode(n, 50.0, dt)
+        return n, dt, dyn, mode
+
+    def _call(self, mock_ls, dyn, mode, dt, **kwargs):
+        mock_ls.SimInspiralESIGMADynamics.return_value = dyn
+        mock_ls.SimInspiralESIGMAModeFromDynamics.return_value = mode
+        defaults = dict(
+            mass1=10.0,
+            mass2=10.0,
+            f_lower=20.0,
+            delta_t=dt,
+            modes_to_use=[(2, 2)],  # always explicit — avoids mutating the default
+            include_conjugate_modes=False,
+        )
+        return get_inspiral_esigma_modes(**{**defaults, **kwargs})
+
+    def test_returned_modes_contain_requested_key(self, ctx):
+        n, dt, dyn, mode = ctx
+        with patch("esigmapy.generator.ls") as mock_ls:
+            modes = self._call(mock_ls, dyn, mode, dt)
+        assert (2, 2) in modes
+
+    def test_conjugate_modes_added_when_flag_true(self, ctx):
+        n, dt, dyn, mode = ctx
+        with patch("esigmapy.generator.ls") as mock_ls:
+            modes = self._call(mock_ls, dyn, mode, dt, include_conjugate_modes=True)
+        assert (2, -2) in modes
+
+    def test_conjugate_modes_absent_when_flag_false(self, ctx):
+        n, dt, dyn, mode = ctx
+        with patch("esigmapy.generator.ls") as mock_ls:
+            modes = self._call(mock_ls, dyn, mode, dt, include_conjugate_modes=False)
+        assert (2, -2) not in modes
+
+    def test_returns_pycbc_timeseries_by_default(self, ctx):
+        n, dt, dyn, mode = ctx
+        with patch("esigmapy.generator.ls") as mock_ls:
+            modes = self._call(mock_ls, dyn, mode, dt)
+        assert isinstance(modes[(2, 2)], pt.TimeSeries)
+
+    def test_returns_numpy_array_when_not_pycbc(self, ctx):
+        n, dt, dyn, mode = ctx
+        with patch("esigmapy.generator.ls") as mock_ls:
+            t_arr, modes = self._call(
+                mock_ls, dyn, mode, dt, return_pycbc_timeseries=False
+            )
+        assert isinstance(modes[(2, 2)], np.ndarray)
+
+    def test_distance_converted_from_mpc_to_si(self, ctx):
+        n, dt, dyn, mode = ctx
+        distance_mpc = 200.0
+        with patch("esigmapy.generator.ls") as mock_ls:
+            self._call(mock_ls, dyn, mode, dt, distance=distance_mpc)
+            passed_distance = mock_ls.SimInspiralESIGMAModeFromDynamics.call_args.args[
+                -1
+            ]
+        assert passed_distance == pytest.approx(
+            distance_mpc * 1e6 * lal.PC_SI, rel=1e-12
+        )
+
+    def test_mode_generation_called_once_per_mode(self, ctx):
+        n, dt, dyn, mode = ctx
+        with patch("esigmapy.generator.ls") as mock_ls:
+            self._call(
+                mock_ls,
+                dyn,
+                mode,
+                dt,
+                modes_to_use=[(2, 2), (3, 3)],
+                include_conjugate_modes=False,
+            )
+        assert mock_ls.SimInspiralESIGMAModeFromDynamics.call_count == 2
+
+
+# ---------------------------------------------------------------------------
+# eccentricity_at_reference_frequency — Tier 2 (LAL mocked)
+# ---------------------------------------------------------------------------
+
+
+class TestEccentricityAtReferenceFrequencyWithMock:
+    """The function finds the x index closest to x_reference and returns e there.
+    Tests verify that index-finding logic and the return value."""
+
+    @staticmethod
+    def _setup(n=100, target_idx=60, M=20.0, dt=1.0 / 4096):
+        x_arr = np.linspace(0.05, 0.5, n)
+        e_arr = np.linspace(0.20, 0.01, n)
+        # Choose f_reference so x_reference == x_arr[target_idx] exactly
+        f_ref = x_arr[target_idx] ** 1.5 / (np.pi * M * lal.MTSUN_SI)
+        dyn = (
+            _lal_ts(np.arange(n, dtype=float) * dt),
+            _lal_ts(x_arr),
+            _lal_ts(e_arr),
+            _lal_ts(np.zeros(n)),
+            _lal_ts(np.zeros(n)),
+            _lal_ts(np.ones(n) * 100.0),
+            _lal_ts(np.ones(n) * 10.0),
+            _lal_ts(np.zeros(n)),
+        )
+        return dyn, x_arr, e_arr, f_ref
+
+    def test_returns_eccentricity_at_x_reference(self):
+        M, target_idx = 20.0, 60
+        dyn, x_arr, e_arr, f_ref = self._setup(target_idx=target_idx, M=M)
+        with patch("esigmapy.generator.ls") as mock_ls:
+            mock_ls.SimInspiralESIGMADynamics.return_value = dyn
+            result = eccentricity_at_reference_frequency(
+                M / 2, M / 2, 0, 0, 0.1, 0, 20.0, 4096, f_ref
+            )
+        assert result == pytest.approx(e_arr[target_idx], rel=1e-10)
+
+    def test_returns_float(self):
+        M = 20.0
+        dyn, _, _, f_ref = self._setup(M=M)
+        with patch("esigmapy.generator.ls") as mock_ls:
+            mock_ls.SimInspiralESIGMADynamics.return_value = dyn
+            result = eccentricity_at_reference_frequency(
+                M / 2, M / 2, 0, 0, 0.1, 0, 20.0, 4096, f_ref
+            )
+        assert np.isscalar(result) or result.ndim == 0
diff --git a/tests/test_legacy.py b/tests/test_legacy.py
new file mode 100644
index 0000000..7cf887a
--- /dev/null
+++ b/tests/test_legacy.py
@@ -0,0 +1,267 @@
+import pytest
+
+from esigmapy.legacy import FitMOmegaIMRAttachmentNonSpinning as Fit
+
+_BASE = 1.0 / 6**1.5
+
+
+@pytest.fixture(autouse=True)
+def reset_called_once():
+    """Reset the class-level flag before each test so logging fires on first call."""
+    Fit.called_once = False
+    yield
+    Fit.called_once = False
+
+
+# ---------------------------------------------------------------------------
+# fit_quadratic_poly
+# ---------------------------------------------------------------------------
+
+
+class TestFitQuadraticPoly:
+    def test_zero_coeffs_returns_base(self):
+        assert Fit.fit_quadratic_poly(0.25, [0.0, 0.0]) == pytest.approx(_BASE)
+
+    def test_known_value(self):
+        a1, a2, eta = 1.0, 2.0, 0.25
+        expected = _BASE * (1.0 + a1 * eta + a2 * eta**2)
+        assert Fit.fit_quadratic_poly(eta, [a1, a2]) == pytest.approx(
+            expected, rel=1e-14
+        )
+
+    def test_wrong_coeff_count_raises(self):
+        with pytest.raises(AssertionError):
+            Fit.fit_quadratic_poly(0.25, [1.0])
+
+    def test_positive(self):
+        assert Fit.fit_quadratic_poly(0.25, [0.0, 0.0]) > 0
+
+
+# ---------------------------------------------------------------------------
+# fit_cubic_poly
+# ---------------------------------------------------------------------------
+
+
+class TestFitCubicPoly:
+    def test_zero_coeffs_returns_base(self):
+        assert Fit.fit_cubic_poly(0.25, [0.0, 0.0, 0.0]) == pytest.approx(_BASE)
+
+    def test_known_value(self):
+        a1, a2, a3, eta = 2.0, -1.0, 0.5, 0.2
+        expected = _BASE * (1.0 + a1 * eta + a2 * eta**2 + a3 * eta**3)
+        assert Fit.fit_cubic_poly(eta, [a1, a2, a3]) == pytest.approx(
+            expected, rel=1e-14
+        )
+
+    def test_wrong_coeff_count_raises(self):
+        with pytest.raises(AssertionError):
+            Fit.fit_cubic_poly(0.25, [1.0, 2.0])
+
+
+# ---------------------------------------------------------------------------
+# fit_ratio_poly_44
+# ---------------------------------------------------------------------------
+
+
+class TestFitRatioPoly44:
+    _Z6 = [0.0] * 6
+
+    def test_zero_coeffs_returns_base(self):
+        assert Fit.fit_ratio_poly_44(0.25, self._Z6) == pytest.approx(_BASE)
+
+    def test_equal_num_denom_returns_base(self):
+        # num == denom → ratio == 1 → result == _BASE regardless of coeffs
+        cs = [1.0, 2.0, 3.0, 1.0, 2.0, 3.0]
+        assert Fit.fit_ratio_poly_44(0.25, cs) == pytest.approx(_BASE, rel=1e-14)
+
+    def test_known_value(self):
+        a1, a2, a3, b1, b2, b3, eta = 1.0, 0.0, 0.0, 0.5, 0.0, 0.0, 0.2
+        expected = _BASE * (1.0 + a1 * eta) / (1.0 + b1 * eta)
+        assert Fit.fit_ratio_poly_44(eta, [a1, a2, a3, b1, b2, b3]) == pytest.approx(
+            expected, rel=1e-14
+        )
+
+    def test_wrong_coeff_count_raises(self):
+        with pytest.raises(AssertionError):
+            Fit.fit_ratio_poly_44(0.25, [1.0] * 5)
+
+
+# ---------------------------------------------------------------------------
+# fit_ratio_sqrt_poly_44
+# ---------------------------------------------------------------------------
+
+
+class TestFitRatioSqrtPoly44:
+    _Z6 = [0.0] * 6
+
+    def test_zero_coeffs_returns_base(self):
+        assert Fit.fit_ratio_sqrt_poly_44(0.25, self._Z6) == pytest.approx(_BASE)
+
+    def test_uses_sqrt_eta_not_eta(self):
+        # With only a1 nonzero: sqrt variant scales by eta^0.5, linear by eta
+        coeffs = [1.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+        eta = 0.16  # sqrt(0.16) = 0.4, so results differ
+        assert Fit.fit_ratio_sqrt_poly_44(eta, coeffs) != pytest.approx(
+            Fit.fit_ratio_poly_44(eta, coeffs)
+        )
+
+    def test_known_value(self):
+        a1, eta = 1.0, 0.16
+        s_eta = eta**0.5
+        expected = _BASE * (1.0 + a1 * s_eta)
+        assert Fit.fit_ratio_sqrt_poly_44(eta, [a1, 0, 0, 0, 0, 0]) == pytest.approx(
+            expected, rel=1e-14
+        )
+
+    def test_wrong_coeff_count_raises(self):
+        with pytest.raises(AssertionError):
+            Fit.fit_ratio_sqrt_poly_44(0.25, [1.0] * 4)
+
+
+# ---------------------------------------------------------------------------
+# fit_ratio_sqrt_hyb1_poly_44 — documents copy-paste bug
+# ---------------------------------------------------------------------------
+
+
+class TestFitRatioSqrtHyb1Poly44:
+    """fit_ratio_sqrt_hyb1_poly_44 computes s_eta but never uses it;
+    the polynomial uses eta instead, making it identical to fit_ratio_poly_44."""
+
+    _Z6 = [0.0] * 6
+
+    def test_zero_coeffs_returns_base(self):
+        assert Fit.fit_ratio_sqrt_hyb1_poly_44(0.25, self._Z6) == pytest.approx(_BASE)
+
+    def test_identical_to_fit_ratio_poly_44(self):
+        # Bug: s_eta is assigned but eta is used in the body
+        coeffs = [1.0, 0.5, 0.2, 0.3, 0.1, 0.05]
+        eta = 0.2
+        assert Fit.fit_ratio_sqrt_hyb1_poly_44(eta, coeffs) == pytest.approx(
+            Fit.fit_ratio_poly_44(eta, coeffs), rel=1e-14
+        )
+
+    def test_wrong_coeff_count_raises(self):
+        with pytest.raises(AssertionError):
+            Fit.fit_ratio_sqrt_hyb1_poly_44(0.25, [1.0] * 7)
+
+
+# ---------------------------------------------------------------------------
+# fit_ratio_poly_43
+# ---------------------------------------------------------------------------
+
+
+class TestFitRatioPoly43:
+    _Z5 = [0.0] * 5
+
+    def test_zero_coeffs_returns_base(self):
+        assert Fit.fit_ratio_poly_43(0.25, self._Z5) == pytest.approx(_BASE)
+
+    def test_known_value(self):
+        a1, a2, a3, b1, b2, eta = 2.0, 0.0, 0.0, 1.0, 0.0, 0.1
+        expected = _BASE * (1.0 + a1 * eta) / (1.0 + b1 * eta)
+        assert Fit.fit_ratio_poly_43(eta, [a1, a2, a3, b1, b2]) == pytest.approx(
+            expected, rel=1e-14
+        )
+
+    def test_wrong_coeff_count_raises(self):
+        with pytest.raises(AssertionError):
+            Fit.fit_ratio_poly_43(0.25, [1.0] * 6)
+
+
+# ---------------------------------------------------------------------------
+# fit_ratio_sqrt_poly_43
+# ---------------------------------------------------------------------------
+
+
+class TestFitRatioSqrtPoly43:
+    _Z5 = [0.0] * 5
+
+    def test_zero_coeffs_returns_base(self):
+        assert Fit.fit_ratio_sqrt_poly_43(0.25, self._Z5) == pytest.approx(_BASE)
+
+    def test_uses_sqrt_eta_not_eta(self):
+        coeffs = [1.0, 0.0, 0.0, 0.0, 0.0]
+        eta = 0.16
+        assert Fit.fit_ratio_sqrt_poly_43(eta, coeffs) != pytest.approx(
+            Fit.fit_ratio_poly_43(eta, coeffs)
+        )
+
+    def test_known_value(self):
+        a1, eta = 1.0, 0.16
+        s_eta = eta**0.5
+        expected = _BASE * (1.0 + a1 * s_eta)
+        assert Fit.fit_ratio_sqrt_poly_43(eta, [a1, 0, 0, 0, 0]) == pytest.approx(
+            expected, rel=1e-14
+        )
+
+    def test_wrong_coeff_count_raises(self):
+        with pytest.raises(AssertionError):
+            Fit.fit_ratio_sqrt_poly_43(0.25, [1.0] * 6)
+
+
+# ---------------------------------------------------------------------------
+# fit_ratio_sqrt_hyb1_poly_43
+# ---------------------------------------------------------------------------
+
+
+class TestFitRatioSqrtHyb1Poly43:
+    _Z5 = [0.0] * 5
+
+    def test_zero_coeffs_returns_base(self):
+        assert Fit.fit_ratio_sqrt_hyb1_poly_43(0.25, self._Z5) == pytest.approx(_BASE)
+
+    def test_known_value(self):
+        # numerator uses eta * s_eta = eta^(3/2)
+        a1, a2, a3, b1, b2, eta = 2.0, 0.0, 0.0, 0.0, 0.0, 0.25
+        s_eta = eta**0.5
+        expected = _BASE * (1.0 + a1 * eta * s_eta)
+        assert Fit.fit_ratio_sqrt_hyb1_poly_43(
+            eta, [a1, a2, a3, b1, b2]
+        ) == pytest.approx(expected, rel=1e-14)
+
+    def test_distinct_from_fit_ratio_poly_43(self):
+        # numerator exponents differ (3/2, 5/2, 7/2 vs 1, 2, 3) → different results
+        coeffs = [1.0, 0.0, 0.0, 0.0, 0.0]
+        eta = 0.16
+        assert Fit.fit_ratio_sqrt_hyb1_poly_43(eta, coeffs) != pytest.approx(
+            Fit.fit_ratio_poly_43(eta, coeffs)
+        )
+
+    def test_wrong_coeff_count_raises(self):
+        with pytest.raises(AssertionError):
+            Fit.fit_ratio_sqrt_hyb1_poly_43(0.25, [1.0] * 4)
+
+
+# ---------------------------------------------------------------------------
+# fit_ratio_poly_34
+# ---------------------------------------------------------------------------
+
+
+class TestFitRatioPoly34:
+    _Z5 = [0.0] * 5
+
+    def test_zero_coeffs_returns_base(self):
+        assert Fit.fit_ratio_poly_34(0.25, self._Z5) == pytest.approx(_BASE)
+
+    def test_known_value(self):
+        a1, a2, b1, b2, b3, eta = 1.0, 0.5, 2.0, 0.0, 0.0, 0.2
+        expected = _BASE * (1.0 + a1 * eta + a2 * eta**2) / (1.0 + b1 * eta)
+        assert Fit.fit_ratio_poly_34(eta, [a1, a2, b1, b2, b3]) == pytest.approx(
+            expected, rel=1e-14
+        )
+
+    def test_wrong_coeff_count_raises(self):
+        with pytest.raises(AssertionError):
+            Fit.fit_ratio_poly_34(0.25, [1.0] * 6)
+
+
+# ---------------------------------------------------------------------------
+# Logging patch smoke test
+# ---------------------------------------------------------------------------
+
+
+class TestLoggingPatch:
+    def test_first_call_does_not_raise(self):
+        # called_once is False (reset by fixture) → logging.info fires
+        # conftest.py patches logging into legacy's namespace → no NameError
+        Fit.fit_quadratic_poly(0.25, [0.0, 0.0])
diff --git a/tests/test_mr_generator.py b/tests/test_mr_generator.py
new file mode 100644
index 0000000..fbb88d1
--- /dev/null
+++ b/tests/test_mr_generator.py
@@ -0,0 +1,239 @@
+import numpy as np
+import pytest
+from unittest.mock import MagicMock, patch
+
+from esigmapy.mr_generator import (
+    check_available_mr_approximants,
+    _generate_lalsim_modes,
+    get_mr_modes,
+    LALSIM_APPROXIMANTS,
+    PYSEOBNR_APPROXIMANTS,
+    SUPPORTED_MR_APPROXIMANTS,
+)
+
+
+class TestCheckAvailableMrApproximants:
+    def test_lalsim_approximants_valid(self):
+        for approx in LALSIM_APPROXIMANTS:
+            check_available_mr_approximants(approx)  # must not raise
+
+    def test_pyseobnr_approximants_valid(self):
+        for approx in PYSEOBNR_APPROXIMANTS:
+            check_available_mr_approximants(approx)  # must not raise
+
+    def test_unsupported_approximant_raises(self):
+        with pytest.raises(IOError, match="cannot generate"):
+            check_available_mr_approximants("IMRPhenomD")
+
+    def test_empty_string_raises(self):
+        with pytest.raises(IOError):
+            check_available_mr_approximants("")
+
+    def test_case_sensitive(self):
+        with pytest.raises(IOError):
+            check_available_mr_approximants("nrsur7dq4")
+
+
+class TestApproximantConstants:
+    def test_supported_is_union_of_lalsim_and_pyseobnr(self):
+        assert set(SUPPORTED_MR_APPROXIMANTS) == set(LALSIM_APPROXIMANTS) | set(
+            PYSEOBNR_APPROXIMANTS
+        )
+
+    def test_lalsim_and_pyseobnr_disjoint(self):
+        assert not set(LALSIM_APPROXIMANTS) & set(PYSEOBNR_APPROXIMANTS)
+
+    def test_lalsim_contains_nrsur7dq4(self):
+        assert "NRSur7dq4" in LALSIM_APPROXIMANTS
+
+    def test_lalsim_contains_seobnrv4phm(self):
+        assert "SEOBNRv4PHM" in LALSIM_APPROXIMANTS
+
+
+# ---------------------------------------------------------------------------
+# _generate_lalsim_modes — Tier 2 (LAL mocked)
+# ---------------------------------------------------------------------------
+
+
+def _mode_node(l, m, data, nxt=None):
+    """Build one node of the SphHarmTimeSeries linked list that LAL returns."""
+    node = MagicMock()
+    node.l = l
+    node.m = m
+    node.mode.data.data = np.asarray(data, dtype=complex)
+    node.next = nxt
+    return node
+
+
+class TestGenerateLalsimModesWithMock:
+    """Tests for _generate_lalsim_modes: mode filtering and data extraction."""
+
+    _N = 64
+
+    def test_bad_approximant_raises_ioerror(self):
+        # Real lalsimulation does not have this attribute → AttributeError → IOError
+        with pytest.raises(IOError, match="not available"):
+            _generate_lalsim_modes(
+                10,
+                10,
+                20,
+                1.0 / 4096,
+                0,
+                0,
+                0,
+                1,
+                20,
+                modes_to_use=[(2, 2)],
+                approximant="NotARealApproximant99999",
+            )
+
+    def test_returns_only_requested_modes(self):
+        # Linked list has (2,2) and (3,3); only (2,2) is requested
+        node_33 = _mode_node(3, 3, np.ones(self._N, dtype=complex))
+        node_22 = _mode_node(2, 2, np.ones(self._N, dtype=complex), nxt=node_33)
+        with patch("esigmapy.mr_generator.ls") as mock_ls:
+            mock_ls.NRSur7dq4 = 900
+            mock_ls.SimInspiralChooseTDModes.return_value = node_22
+            result = _generate_lalsim_modes(
+                10,
+                10,
+                20,
+                1.0 / 4096,
+                0,
+                0,
+                0,
+                1,
+                20,
+                modes_to_use=[(2, 2)],
+                approximant="NRSur7dq4",
+            )
+        assert set(result.keys()) == {(2, 2)}
+
+    def test_returns_all_requested_modes(self):
+        # Linked list has (2,2) and (3,3); both requested
+        node_33 = _mode_node(3, 3, np.ones(self._N, dtype=complex))
+        node_22 = _mode_node(2, 2, np.ones(self._N, dtype=complex), nxt=node_33)
+        with patch("esigmapy.mr_generator.ls") as mock_ls:
+            mock_ls.NRSur7dq4 = 900
+            mock_ls.SimInspiralChooseTDModes.return_value = node_22
+            result = _generate_lalsim_modes(
+                10,
+                10,
+                20,
+                1.0 / 4096,
+                0,
+                0,
+                0,
+                1,
+                20,
+                modes_to_use=[(2, 2), (3, 3)],
+                approximant="NRSur7dq4",
+            )
+        assert set(result.keys()) == {(2, 2), (3, 3)}
+
+    def test_mode_data_is_numpy_array(self):
+        node = _mode_node(2, 2, np.ones(self._N, dtype=complex))
+        node.next = None
+        with patch("esigmapy.mr_generator.ls") as mock_ls:
+            mock_ls.NRSur7dq4 = 900
+            mock_ls.SimInspiralChooseTDModes.return_value = node
+            result = _generate_lalsim_modes(
+                10,
+                10,
+                20,
+                1.0 / 4096,
+                0,
+                0,
+                0,
+                1,
+                20,
+                modes_to_use=[(2, 2)],
+                approximant="NRSur7dq4",
+            )
+        assert isinstance(result[(2, 2)], np.ndarray)
+
+    def test_unrequested_mode_in_list_is_ignored(self):
+        # Linked list has (4,4) only; request is (2,2) — result should be empty
+        node = _mode_node(4, 4, np.ones(self._N, dtype=complex))
+        node.next = None
+        with patch("esigmapy.mr_generator.ls") as mock_ls:
+            mock_ls.NRSur7dq4 = 900
+            mock_ls.SimInspiralChooseTDModes.return_value = node
+            result = _generate_lalsim_modes(
+                10,
+                10,
+                20,
+                1.0 / 4096,
+                0,
+                0,
+                0,
+                1,
+                20,
+                modes_to_use=[(2, 2)],
+                approximant="NRSur7dq4",
+            )
+        assert (2, 2) not in result
+
+
+# ---------------------------------------------------------------------------
+# get_mr_modes — Tier 2 (internal helpers mocked)
+# ---------------------------------------------------------------------------
+
+
+class TestGetMrModesWithMock:
+    """Tests for the Python routing and defaulting logic in get_mr_modes."""
+
+    _FAKE = {(2, 2): np.ones(64, dtype=complex)}
+
+    def test_f_ref_defaults_to_f_lower(self):
+        with patch("esigmapy.mr_generator._generate_lalsim_modes") as mock_gen:
+            mock_gen.return_value = self._FAKE
+            get_mr_modes(
+                10, 10, 20.0, 1.0 / 4096, modes_to_use=[(2, 2)], approximant="NRSur7dq4"
+            )
+        assert mock_gen.call_args.kwargs["f_ref"] == 20.0
+
+    def test_coa_phase_defaults_to_zero(self):
+        with patch("esigmapy.mr_generator._generate_lalsim_modes") as mock_gen:
+            mock_gen.return_value = self._FAKE
+            get_mr_modes(
+                10,
+                10,
+                20.0,
+                1.0 / 4096,
+                coa_phase=None,
+                modes_to_use=[(2, 2)],
+                approximant="NRSur7dq4",
+            )
+        assert mock_gen.call_args.kwargs["coa_phase"] == 0.0
+
+    def test_lalsim_approximant_routes_to_lalsim(self):
+        with patch(
+            "esigmapy.mr_generator._generate_lalsim_modes"
+        ) as mock_lalsim, patch(
+            "esigmapy.mr_generator._generate_pyseobnr_modes"
+        ) as mock_pyseobnr:
+            mock_lalsim.return_value = self._FAKE
+            get_mr_modes(
+                10, 10, 20.0, 1.0 / 4096, modes_to_use=[(2, 2)], approximant="NRSur7dq4"
+            )
+        mock_lalsim.assert_called_once()
+        mock_pyseobnr.assert_not_called()
+
+    def test_pyseobnr_approximant_routes_to_pyseobnr(self):
+        with patch(
+            "esigmapy.mr_generator._generate_pyseobnr_modes"
+        ) as mock_pyseobnr, patch(
+            "esigmapy.mr_generator._generate_lalsim_modes"
+        ) as mock_lalsim:
+            mock_pyseobnr.return_value = self._FAKE
+            get_mr_modes(
+                10,
+                10,
+                20.0,
+                1.0 / 4096,
+                modes_to_use=[(2, 2)],
+                approximant="SEOBNRv5HM",
+            )
+        mock_pyseobnr.assert_called_once()
+        mock_lalsim.assert_not_called()
diff --git a/tests/test_utils.py b/tests/test_utils.py
new file mode 100644
index 0000000..b2f6d96
--- /dev/null
+++ b/tests/test_utils.py
@@ -0,0 +1,246 @@
+import numpy as np
+import pytest
+import lal
+import pycbc.types as pt
+
+from esigmapy.utils import (
+    f22_from_x,
+    x_from_f22,
+    f_ISCO_spin,
+    get_peak_freqs,
+    get_polarizations_from_multipoles,
+)
+
+_MTSUN = lal.MTSUN_SI
+_PI = np.pi
+
+
+# ---------------------------------------------------------------------------
+# f22_from_x / x_from_f22
+# ---------------------------------------------------------------------------
+
+
+class TestF22X:
+    @pytest.mark.parametrize(
+        "x,M",
+        [
+            (0.05, 20.0),
+            (0.1, 50.0),
+            (0.2, 100.0),
+            (0.3, 200.0),
+        ],
+    )
+    def test_roundtrip_x_to_f_to_x(self, x, M):
+        assert x_from_f22(f22_from_x(x, M), M) == pytest.approx(x, rel=1e-12)
+
+    @pytest.mark.parametrize(
+        "f,M",
+        [
+            (20.0, 30.0),
+            (50.0, 50.0),
+            (100.0, 100.0),
+        ],
+    )
+    def test_roundtrip_f_to_x_to_f(self, f, M):
+        assert f22_from_x(x_from_f22(f, M), M) == pytest.approx(f, rel=1e-12)
+
+    def test_f22_mass_scaling(self):
+        # Doubling total mass at fixed x halves the GW frequency
+        x = 0.1
+        assert f22_from_x(x, 40.0) == pytest.approx(f22_from_x(x, 20.0) / 2, rel=1e-12)
+
+    def test_x_mass_scaling(self):
+        # At fixed f, doubling M scales x by 2^(2/3)
+        f = 100.0
+        ratio = x_from_f22(f, 40.0) / x_from_f22(f, 20.0)
+        assert ratio == pytest.approx(2.0 ** (2.0 / 3.0), rel=1e-12)
+
+    def test_f22_positive(self):
+        assert f22_from_x(0.1, 50.0) > 0
+
+    def test_x_positive(self):
+        assert x_from_f22(100.0, 50.0) > 0
+
+    def test_f22_known_value(self):
+        x, M = 0.1, 50.0
+        expected = x**1.5 / (M * _MTSUN * _PI)
+        assert f22_from_x(x, M) == pytest.approx(expected, rel=1e-14)
+
+    def test_x_known_value(self):
+        f, M = 100.0, 50.0
+        expected = (M * _MTSUN * _PI * f) ** (2.0 / 3.0)
+        assert x_from_f22(f, M) == pytest.approx(expected, rel=1e-14)
+
+
+# ---------------------------------------------------------------------------
+# f_ISCO_spin
+# ---------------------------------------------------------------------------
+
+
+class TestFISCOSpin:
+    def test_positive_for_zero_spin(self):
+        assert f_ISCO_spin(20.0, 20.0, 0.0, 0.0) > 0
+
+    def test_positive_for_nonzero_spin(self):
+        assert f_ISCO_spin(10.0, 20.0, 0.3, -0.2) > 0
+
+    def test_mass_scaling_equal_mass_zero_spin(self):
+        # For equal-mass, zero-spin: frequency ∝ 1/M (eta fixed, Scap=0, SS=1)
+        f1 = f_ISCO_spin(10.0, 10.0, 0.0, 0.0)
+        f2 = f_ISCO_spin(20.0, 20.0, 0.0, 0.0)
+        assert f1 / f2 == pytest.approx(2.0, rel=1e-6)
+
+    def test_prograde_spin_raises_isco_frequency(self):
+        f_zero = f_ISCO_spin(10.0, 10.0, 0.0, 0.0)
+        f_pro = f_ISCO_spin(10.0, 10.0, 0.5, 0.5)
+        assert f_pro > f_zero
+
+    def test_retrograde_spin_lowers_isco_frequency(self):
+        f_zero = f_ISCO_spin(10.0, 10.0, 0.0, 0.0)
+        f_retro = f_ISCO_spin(10.0, 10.0, -0.5, -0.5)
+        assert f_retro < f_zero
+
+    def test_spin_monotone_in_spin_magnitude(self):
+        # Larger prograde spin → higher ISCO frequency
+        f_lo = f_ISCO_spin(10.0, 10.0, 0.3, 0.3)
+        f_hi = f_ISCO_spin(10.0, 10.0, 0.7, 0.7)
+        assert f_hi > f_lo
+
+
+# ---------------------------------------------------------------------------
+# get_peak_freqs
+# ---------------------------------------------------------------------------
+
+
+class TestGetPeakFreqs:
+    @staticmethod
+    def _ts(data, dt=1.0 / 4096):
+        return pt.TimeSeries(data.astype(float), delta_t=dt)
+
+    def test_no_peaks_for_monotone_increasing(self):
+        ts = self._ts(np.linspace(10.0, 200.0, 500))
+        times, peaks = get_peak_freqs(ts)
+        assert len(peaks) == 0
+        assert len(times) == 0
+
+    def test_no_peaks_for_constant_series(self):
+        ts = self._ts(np.full(100, 42.0))
+        times, peaks = get_peak_freqs(ts)
+        assert len(peaks) == 0
+
+    def test_single_isolated_peak(self):
+        data = np.zeros(101)
+        data[50] = 1.0
+        ts = self._ts(data)
+        times, peaks = get_peak_freqs(ts)
+        assert len(peaks) == 1
+        assert peaks[0] == pytest.approx(1.0)
+
+    def test_three_isolated_peaks_count(self):
+        data = np.zeros(301)
+        data[50] = 1.0
+        data[150] = 2.0
+        data[250] = 3.0
+        ts = self._ts(data)
+        _, peaks = get_peak_freqs(ts)
+        assert len(peaks) == 3
+
+    def test_three_isolated_peaks_values(self):
+        data = np.zeros(301)
+        data[50] = 1.0
+        data[150] = 2.0
+        data[250] = 3.0
+        ts = self._ts(data)
+        _, peaks = get_peak_freqs(ts)
+        assert sorted(peaks) == pytest.approx([1.0, 2.0, 3.0])
+
+    def test_boundary_indices_ignored(self):
+        # Peaks at index 0 and n-1 are never reported
+        data = np.zeros(51)
+        data[0] = 10.0
+        data[-1] = 10.0
+        data[25] = 5.0
+        ts = self._ts(data)
+        _, peaks = get_peak_freqs(ts)
+        assert len(peaks) == 1
+        assert peaks[0] == pytest.approx(5.0)
+
+    def test_peak_times_correspond_to_peak_values(self):
+        dt = 1.0 / 1024
+        data = np.zeros(201)
+        data[100] = 7.0
+        ts = self._ts(data, dt=dt)
+        times, peaks = get_peak_freqs(ts)
+        assert len(times) == 1
+        assert times[0] == pytest.approx(100 * dt, rel=1e-6)
+
+
+# ---------------------------------------------------------------------------
+# get_polarizations_from_multipoles
+# ---------------------------------------------------------------------------
+
+
+class TestGetPolarizationsFromMultipoles:
+    def test_zero_mode_gives_zero_polarizations(self):
+        modes = {(2, 2): np.array([0.0 + 0j, 0.0 + 0j])}
+        hp, hc = get_polarizations_from_multipoles(
+            modes, inclination=0.3, coa_phase=0.0
+        )
+        assert np.allclose(hp, 0.0)
+        assert np.allclose(hc, 0.0)
+
+    def test_linear_amplitude_scaling(self):
+        modes1 = {(2, 2): np.array([1.0 + 0j])}
+        modes2 = {(2, 2): np.array([3.0 + 0j])}
+        hp1, hc1 = get_polarizations_from_multipoles(
+            modes1, inclination=0.5, coa_phase=0.0
+        )
+        hp2, hc2 = get_polarizations_from_multipoles(
+            modes2, inclination=0.5, coa_phase=0.0
+        )
+        assert np.allclose(hp2, 3.0 * hp1)
+        assert np.allclose(hc2, 3.0 * hc1)
+
+    def test_polarization_norm_single_mode(self):
+        # hp^2 + hc^2 == |h22|^2 * |Y_{-2,22}|^2  for a single mode
+        inc, coa = 0.5, 0.3
+        h22 = 2.0 + 1.5j
+        modes = {(2, 2): np.array([h22])}
+        hp, hc = get_polarizations_from_multipoles(
+            modes, inclination=inc, coa_phase=coa
+        )
+        ylm = lal.SpinWeightedSphericalHarmonic(inc, coa, -2, 2, 2)
+        expected = abs(h22) ** 2 * abs(ylm) ** 2
+        actual = float(hp[0]) ** 2 + float(hc[0]) ** 2
+        assert actual == pytest.approx(expected, rel=1e-12)
+
+    def test_two_modes_superpose_linearly(self):
+        inc, coa = 0.4, 0.2
+        h22 = np.array([1.0 + 0j])
+        h21 = np.array([0.5 + 0.3j])
+        hp22, hc22 = get_polarizations_from_multipoles({(2, 2): h22}, inc, coa)
+        hp21, hc21 = get_polarizations_from_multipoles({(2, 1): h21}, inc, coa)
+        hp_both, hc_both = get_polarizations_from_multipoles(
+            {(2, 2): h22, (2, 1): h21}, inc, coa
+        )
+        assert np.allclose(hp_both, hp22 + hp21)
+        assert np.allclose(hc_both, hc22 + hc21)
+
+    def test_output_is_real_valued(self):
+        modes = {(2, 2): np.array([1.0 + 2.0j, 3.0 - 1.0j])}
+        hp, hc = get_polarizations_from_multipoles(
+            modes, inclination=0.3, coa_phase=0.0
+        )
+        assert np.all(np.imag(hp) == 0.0)
+        assert np.all(np.imag(hc) == 0.0)
+
+    def test_coa_phase_shifts_polarizations(self):
+        # Different coalescence phases must give different results (not trivially zero-effect)
+        modes = {(2, 2): np.array([1.0 + 0j])}
+        hp0, hc0 = get_polarizations_from_multipoles(
+            modes, inclination=0.5, coa_phase=0.0
+        )
+        hp1, hc1 = get_polarizations_from_multipoles(
+            modes, inclination=0.5, coa_phase=0.5
+        )
+        assert not np.allclose(hp0, hp1) or not np.allclose(hc0, hc1)