diff --git a/gw_eccentricity/eccDefinition.py b/gw_eccentricity/eccDefinition.py
index de8884da..f1bf023f 100644
--- a/gw_eccentricity/eccDefinition.py
+++ b/gw_eccentricity/eccDefinition.py
@@ -8,6 +8,7 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
+from copy import deepcopy
 import warnings
 from .load_data import get_coprecessing_data_dict
 from .utils import peak_time_via_quadratic_fit, check_kwargs_and_set_defaults
@@ -16,6 +17,8 @@
 from .utils import interpolate
 from .utils import get_interpolant
 from .utils import get_default_spline_kwargs
+from .utils import get_rational_fit
+from .utils import get_default_rational_fit_kwargs
 from .utils import debug_message
 from .plot_settings import use_fancy_plotsettings, colorsDict, labelsDict
 from .plot_settings import figWidthsTwoColDict, figHeightsDict
@@ -199,13 +202,30 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2,
         extra_kwargs: dict
             A dictionary of any extra kwargs to be passed. Allowed kwargs
             are:
-            spline_kwargs: dict
+            special_interp_kwargs_for_extrema: dict
+                A dictionary with a single key matching the current
+                `omega_gw_extrema_interpolation_method`. See under
+                `omega_gw_extrema_interpolation_method` for more details on
+                possible values of interpolation methods.
+
+                Example structure:
+                    {
+                        "method": {"param1": value1, "param2": value2},
+                    }
+
+                currently, the available options for "method" are: 
+
+                - "spline": default kwargs are set using
+                  `utils.get_default_spline_kwargs`
+                - "rational_fit": default kwargs are set using
+                  `utils.get_default_rational_fit_kwargs`
+
+            general_interp_kwargs: dict
                 Dictionary of arguments to be passed to the spline
                 interpolation routine
                 (scipy.interpolate.InterpolatedUnivariateSpline) used to
-                compute quantities like omega_gw_pericenters(t) and
-                omega_gw_apocenters(t).
-                Defaults are set using utils.get_default_spline_kwargs
+                interpolate data other than the omega_gw extrema.
+                Defaults are set using `utils.get_default_spline_kwargs`
 
             extrema_finding_kwargs: dict
                 Dictionary of arguments to be passed to the extrema finder,
@@ -305,6 +325,31 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2,
                 eccentricity is the cause, and set the returned eccentricity
                 and mean anomaly to zero.
                 USE THIS WITH CAUTION!
+
+            omega_gw_extrema_interpolation_method : str, default="rational_fit"
+                Specifies the method used to build the interpolations for
+                `omega_gw_pericenters_interp(t)` or
+                `omega_gw_apocenters_interp(t)`. The available options are:
+
+                - `spline`: Uses `scipy.interpolate.InterpolatedUnivariateSpline`.
+                    - Best suited for cleaner data, such as when waveform modes
+                      are generated using models like SEOB or TEOB.
+                    - Faster to construct and evaluate.
+                    - Since it fits through every data point, it may exhibit
+                      oscillatory behavior, particularly near the merger,
+                      especially for noisy NR data.
+        
+                - `rational_fit`: Uses `polyrat.StabilizedSKRationalApproximation`.
+                    - Can handle both clean and noisy data, e.g., waveform
+                        modes from numerical simulations.
+                    - Better monotonic behaviour, particularly near the merger.
+                    - Significantly slower compared to the `spline` method.
+                      This is because finding optimal numerator and denominator
+                      degree needs several iterations
+                    - Can suppress pathologies in the waveform that might be
+                      visible with `spline`.
+
+                Default value: `"rational_fit"`.
         """
         self.precessing = precessing
         self.frame = frame
@@ -337,24 +382,52 @@ def __init__(self, dataDict, num_orbits_to_exclude_before_merger=2,
             = self.get_amp_phase_omega_gw(self.dataDict)
         # Sanity check various kwargs and set default values
         self.extra_kwargs = check_kwargs_and_set_defaults(
-            extra_kwargs, self.get_default_extra_kwargs(),
+            deepcopy(extra_kwargs), self.get_default_extra_kwargs(),
             "extra_kwargs",
             "eccDefinition.get_default_extra_kwargs()")
+        self.debug_level = self.extra_kwargs["debug_level"]
         self.extrema_finding_kwargs = check_kwargs_and_set_defaults(
             self.extra_kwargs['extrema_finding_kwargs'],
             self.get_default_extrema_finding_kwargs(min_width_for_extrema),
             "extrema_finding_kwargs",
             "eccDefinition.get_default_extrema_finding_kwargs()")
-        self.spline_kwargs = check_kwargs_and_set_defaults(
-            self.extra_kwargs["spline_kwargs"],
+        # set omega extrema interpolation method
+        self.omega_gw_extrema_interpolation_method = \
+            self.extra_kwargs["omega_gw_extrema_interpolation_method"]
+        # check extrema interpolation method
+        self.available_omega_gw_extrema_interpolation_methods = \
+            self.get_available_omega_gw_extrema_interpolation_methods()
+        if self.omega_gw_extrema_interpolation_method \
+            not in self.available_omega_gw_extrema_interpolation_methods:
+            raise Exception(
+                "Unknown omega_gw_extrema_interpolation_method "
+                f"`{self.omega_gw_extrema_interpolation_method}`. Should be "
+                f"of {self.available_omega_gw_extrema_interpolation_methods.keys()}")
+        # check extrema interp kwargs
+        self.check_special_interp_kwargs_for_extrema(extra_kwargs)
+        # set extrema interpolation kwargs
+        self.special_interp_kwargs_for_extrema = check_kwargs_and_set_defaults(
+            self.extra_kwargs[
+                "special_interp_kwargs_for_extrema"][
+                    self.omega_gw_extrema_interpolation_method],
+            self.get_default_special_interp_kwargs_for_extrema(),
+            "special_interp_kwargs_for_extrema",
+            "eccDefinition.get_default_special_interp_kwargs_for_extrema()")
+        # set other variables required for the omega_gw extrema interpolation
+        # method
+        self.set_other_variables_for_extrema_interpolation()
+        # set spline interpolation kwargs to be used for interpolating data
+        # other than the extrema.
+        self.general_interp_kwargs = check_kwargs_and_set_defaults(
+            self.extra_kwargs["general_interp_kwargs"],
             get_default_spline_kwargs(),
-            "spline_kwargs",
+            "general_interp_kwargs",
             "utils.get_default_spline_kwargs()")
         self.available_averaging_methods \
             = self.get_available_omega_gw_averaging_methods()
-        self.debug_level = self.extra_kwargs["debug_level"]
         self.debug_plots = self.extra_kwargs["debug_plots"]
-        self.return_zero_if_small_ecc_failure = self.extra_kwargs["return_zero_if_small_ecc_failure"]
+        self.return_zero_if_small_ecc_failure \
+            = self.extra_kwargs["return_zero_if_small_ecc_failure"]
         # check if there are unrecognized keys in the dataDict
         self.recognized_dataDict_keys = self.get_recognized_dataDict_keys()
         for kw in dataDict.keys():
@@ -980,11 +1053,94 @@ def get_default_extrema_finding_kwargs(self, width):
             "rel_height": 0.5,
             "plateau_size": None}
         return default_extrema_finding_kwargs
+    
+    def get_default_special_interp_kwargs_for_extrema(self):
+        """Get default kwargs to be passed to extrema interpolating method.
+        
+        Accurate interpolation of omega_gw extrema is crucial for obtaining
+        monotonic eccentricity evoultion with time. Depending on the
+        interpolation method, this function returns the default kwargs to be
+        passed to the interpolating function. 
+        """
+        allowed_methods = self.extra_kwargs["special_interp_kwargs_for_extrema"].keys()
+        if self.omega_gw_extrema_interpolation_method == "spline":
+            kwargs = get_default_spline_kwargs()
+        elif self.omega_gw_extrema_interpolation_method == "rational_fit":
+            kwargs = get_default_rational_fit_kwargs()
+        else:
+            raise Exception("Unknown omega_gw extrema interpolation method "
+                            f"`{self.omega_gw_extrema_interpolation_method}`. "
+                            f"Allowed methods are {allowed_methods}")
+        return kwargs
+    
+    def check_special_interp_kwargs_for_extrema(self, extra_kwargs):
+        """Check extrema interp kwargs provided in extra_kwargs.
+        
+        The `special_interp_kwargs_for_extrema` should have exactly one
+        key matching the current `omega_gw_extrema_interpolation_method`.
+        """
+        if (extra_kwargs is not None 
+            and "special_interp_kwargs_for_extrema" in extra_kwargs):
+            # check that kwargs for only one extrema interpolation method is
+            # provided.
+            common_message = (
+                "The input `special_interp_kwargs_for_extrema` must contain "
+                "a single key, matching the current "
+                "`omega_gw_extrema_interpolation_method` "
+                f"'{self.omega_gw_extrema_interpolation_method}', "
+                "with a dictionary of kwargs for "
+                f"'{self.omega_gw_extrema_interpolation_method}' "
+                "as its value.")
+            extrema_interp_keys \
+                = list(extra_kwargs["special_interp_kwargs_for_extrema"].keys())
+            if len(extrema_interp_keys) == 0:
+                raise Exception(
+                    "Dictionay provided via 'special_interp_kwargs_for_extrema'"
+                    " in 'extra_kwargs' can not be empty.\n"
+                    f"{common_message}"
+                    )
+            # check that the key in special_interp_kwargs_for_extrema matches
+            # `omega_gw_extrema_interpolation_method`
+            if  len(extrema_interp_keys) == 1 and (
+                extrema_interp_keys[0] 
+                != self.omega_gw_extrema_interpolation_method):
+                raise Exception(
+                    "`omega_gw_extrema_interpolation_method` is "
+                    f"{self.omega_gw_extrema_interpolation_method} but the "
+                    "kwargs in 'special_interp_kwargs_for_extrema' via "
+                    f"'extra_kwargs' is for {extrema_interp_keys[0]}.\n"
+                    f"{common_message}"
+                )
+            if len(extrema_interp_keys) > 1:
+                raise Exception(f"{common_message}")
+        
+    def set_other_variables_for_extrema_interpolation(self):
+        """Set other variables required for extrema interpolation.
+        
+        Depending on the omega_gw_extrema_interpolation_method, additional
+        varibles may be required for obtaining good interpolant and/or for
+        debugging purpose. These variables may vary from method to method.
+        """
+        if self.omega_gw_extrema_interpolation_method == "rational_fit":
+            # set verbose to debug_level. If verbose is True, then it prints
+            # information of each iteration for rational fits for a given
+            # degree.
+            self.special_interp_kwargs_for_extrema["verbose"] = self.debug_level
+            # keep history of rational fit degree and nonmonotonicity of the
+            # corresponding fits
+            self.rational_fit_nonmonotonicity_history = {
+                "pericenters": {}, "apocenters": {}}
+            # Store degrees used to construct the final fits
+            self.rational_fit_degrees = {
+                "pericenters": None, "apocenters": None}
 
     def get_default_extra_kwargs(self):
         """Defaults for additional kwargs."""
         default_extra_kwargs = {
-            "spline_kwargs": {},
+            "general_interp_kwargs": {},
+            "special_interp_kwargs_for_extrema": {
+                "spline": {},
+                "rational_fit": {}},
             "extrema_finding_kwargs": {},   # Gets overridden in methods like
                                             # eccDefinitionUsingAmplitude
             "debug_level": 0,
@@ -994,6 +1150,7 @@ def get_default_extra_kwargs(self):
             "refine_extrema": False,
             "kwargs_for_fits_methods": {},  # Gets overriden in fits methods
             "return_zero_if_small_ecc_failure": False,
+            "omega_gw_extrema_interpolation_method": "rational_fit"
         }
         return default_extra_kwargs
 
@@ -1285,33 +1442,220 @@ def get_good_extrema(self, pericenters, apocenters,
             pericenters, apocenters)
         return pericenters, apocenters
 
-    def get_interp(self, oldX, oldY, allowExtrapolation=False,
-                   interpolator="spline"):
-        """Get interpolant.
-
-        A wrapper of utils.get_interpolant with check_kwargs=False.
-        This is to make sure that the checking of kwargs is not performed
-        everytime the interpolation function is called. Instead, the kwargs
-        are checked once in the init and passed to the interpolation
-        function without repeating checks.
+    def get_rational_fit_for_extrema(self, x, y, data_name=None):
+        """Get rational fit with adaptive numerator and denominator degree.
+
+        This function ensures that the rational fit uses the optimal degree for
+        the numerator and denominator. It first checks for nonmonotonicity in
+        the fit's derivative and lowers the degrees if necessary. Only if the
+        fit for the initial degree was monotonic (no reduction of degree
+        needed), it attempts to increase the degrees to find a higher-degree
+        fit. If increasing the degrees causes nonmonotonicity, it reverts to
+        the previous valid monotonic fit.
+
+        The initial degrees for the rational fit can be specified through
+        `special_interp_kwargs_for_extrema` in `extra_kwargs` with the key
+        "rational_fit". (See `special_interp_kwargs_for_extrema` under
+        `extra_kwargs` for more details) Default values are provided by
+        `get_default_rational_fit_kwargs`, where both degrees are set to
+        `None`. If both degrees remain `None`, appropriate starting values are
+        determined using `self.get_approximate_degree_for_rational_fit`.
         """
-        return get_interpolant(oldX, oldY, allowExtrapolation, interpolator,
-                               spline_kwargs=self.spline_kwargs,
-                               check_kwargs=False)
-
-    def interp(self, newX, oldX, oldY, allowExtrapolation=False,
-               interpolator="spline"):
-        """Get interpolated values.
+        # make sure that data_name is not None. A data_name is needed to
+        # update the optimal values of the numerator and denominator degrees
+        # used to build the final rational fit
+        if data_name is None:
+            raise Exception(
+                "Please provide a data_name for which to build a "
+                "rational fit. `data_name` can not be None. For example, "
+                "it can be 'apocenters', 'pericenters' or 'omega_gw_average'.")
+        # Set initial degrees for the rational fit if not already specified. If
+        # either `num_degree` or `denom_degree` is specified, use that value to
+        # set the other. If both are unspecified, we use the following
+        # procedure to set appropriate values for `num_degree` and
+        # `denom_degree`.
+        
+        # The optimal degrees for the numerator and denominator can vary with
+        # the eccentricity and duration of the waveform. To estimate these
+        # values, we first count the approximate number of orbits, then use a
+        # predefined mapping from orbit count to optimal degrees, provided by
+        # `get_approximate_degree_for_rational_fit`. 
+        
+        # While this provides a good initial estimate, it may not be optimal,
+        # as the ideal degree can change based on the eccentricity which we do
+        # not know apriori. The true optimal values are determined through
+        # further iterations.
+        if self.special_interp_kwargs_for_extrema["num_degree"] is None \
+            and self.special_interp_kwargs_for_extrema["denom_degree"] is None:
+            self.special_interp_kwargs_for_extrema["num_degree"], \
+                self.special_interp_kwargs_for_extrema["denom_degree"] = \
+                self.get_approximate_degree_for_rational_fit()
+        elif self.special_interp_kwargs_for_extrema["num_degree"] is None:
+            self.special_interp_kwargs_for_extrema["num_degree"] \
+                = self.special_interp_kwargs_for_extrema["denom_degree"]
+        elif self.special_interp_kwargs_for_extrema["denom_degree"] is None:
+            self.special_interp_kwargs_for_extrema["denom_degree"] \
+                = self.special_interp_kwargs_for_extrema["num_degree"]
+
+        rat_fit = get_rational_fit(
+            x, y, rational_fit_kwargs=self.special_interp_kwargs_for_extrema,
+            check_kwargs=False)
+        x_test = np.arange(x[0], x[-1], self.t[1] - self.t[0])
+        # save the degrees for checks at each step of iterations for finding
+        # the optimal degrees
+        old_num_degree = self.special_interp_kwargs_for_extrema["num_degree"]
+        old_denom_degree = self.special_interp_kwargs_for_extrema["denom_degree"]
+
+        # Check for nonmonotonicity and lower degrees if needed
+        fit_is_nonmonotonic = self.check_if_first_derivative_is_not_strictly_monotonic(
+            x_test, rat_fit(x_test), 1.0, data_name)
+        while fit_is_nonmonotonic:
+            if self.special_interp_kwargs_for_extrema["num_degree"] > 1:
+                self.special_interp_kwargs_for_extrema["num_degree"] -= 1
+            if self.special_interp_kwargs_for_extrema["denom_degree"] > 1:
+                self.special_interp_kwargs_for_extrema["denom_degree"] -= 1
+            if self.special_interp_kwargs_for_extrema["num_degree"] == 1 \
+                and self.special_interp_kwargs_for_extrema["denom_degree"] == 1:
+                debug_message(
+                    "Both numerator and denominator degrees are equal to 1 "
+                    "and cannot be lowered further.",
+                    debug_level=self.debug_level, important=False)
+                break
 
-        A wrapper of utils.interpolate with check_kwargs=False for
-        reasons explained in the documentation of get_interp function.
+            debug_message(
+                "Lowering degrees to "
+                "num_degree="
+                f"{self.special_interp_kwargs_for_extrema['num_degree']}, "
+                "denom_degree="
+                f"{self.special_interp_kwargs_for_extrema['denom_degree']} "
+                "and retrying.", debug_level=self.debug_level, important=False)
+
+            # build new fit and check monotonicity
+            rat_fit = get_rational_fit(
+                x, y,
+                rational_fit_kwargs=self.special_interp_kwargs_for_extrema,
+                check_kwargs=False)
+            fit_is_nonmonotonic \
+                = self.check_if_first_derivative_is_not_strictly_monotonic(
+                x_test, rat_fit(x_test), 1.0, data_name)
+        # If fit with initial degree is monotonic, try increasing the degree
+        # for a better fit
+        if not fit_is_nonmonotonic:
+            last_monotonic_rat_fit = rat_fit  # Track last monotonic fit
+            last_monotonic_num_degree = old_num_degree
+            last_monotonic_denom_degree = old_denom_degree
+
+            while not fit_is_nonmonotonic:
+                # Increase the degrees for both numerator and denominator
+                new_num_degree \
+                    = self.special_interp_kwargs_for_extrema["num_degree"] + 1
+                new_denom_degree \
+                    = self.special_interp_kwargs_for_extrema["denom_degree"] + 1
+                self.special_interp_kwargs_for_extrema["num_degree"],\
+                    self.special_interp_kwargs_for_extrema["denom_degree"] \
+                    = new_num_degree, new_denom_degree
+                # build new fit and check monotonicity
+                new_rat_fit = get_rational_fit(
+                    x, y,
+                    rational_fit_kwargs=self.special_interp_kwargs_for_extrema,
+                    check_kwargs=False)
+                fit_is_nonmonotonic \
+                    = self.check_if_first_derivative_is_not_strictly_monotonic(
+                        x_test, new_rat_fit(x_test), 1.0, data_name)
+                if fit_is_nonmonotonic:
+                    # Revert to previous fit and degrees if nonmonotonicity is
+                    # detected
+                    debug_message(
+                        "Increasing degrees caused nonmonotonicity. "
+                        "Reverting to last monotonic fit with "
+                        f"num_degree={last_monotonic_num_degree} and "
+                        f"denom_degree={last_monotonic_denom_degree}.",
+                        debug_level=self.debug_level, important=False)
+                    self.special_interp_kwargs_for_extrema["num_degree"] \
+                        = last_monotonic_num_degree
+                    self.special_interp_kwargs_for_extrema["denom_degree"] \
+                        = last_monotonic_denom_degree
+                    rat_fit = last_monotonic_rat_fit
+                    break
+
+                last_monotonic_rat_fit = new_rat_fit
+                last_monotonic_num_degree = new_num_degree
+                last_monotonic_denom_degree = new_denom_degree
+                rat_fit = new_rat_fit
+
+        # update final degrees used to build the fit
+        self.rational_fit_degrees[data_name] = (
+            self.special_interp_kwargs_for_extrema["num_degree"],
+            self.special_interp_kwargs_for_extrema["denom_degree"])
+        return rat_fit
+
+    def get_approximate_degree_for_rational_fit(self):
+        """Get approximate degree based on the number of extrema.
+
+        Assign degree based on number of extrema found. The degree is increased
+        as the number of extrema increases.
+        """
+        # TODO: Optimize this.
+        # assign degree based on the number of extrema if user provided
+        # degree is None.
+        approximate_num_orbits = max(len(self.pericenters_location),
+                                     len(self.apocenters_location))
+        if approximate_num_orbits <= 5:
+            num_degree, denom_degree = 1, 1
+        elif approximate_num_orbits <= 20:
+            num_degree, denom_degree =  2, 2
+        elif approximate_num_orbits <= 50:
+            num_degree, denom_degree = 3, 3
+        elif approximate_num_orbits <= 100:
+            num_degree, denom_degree = 4, 4
+        elif approximate_num_orbits <= 150:
+            num_degree, denom_degree = 5, 5
+        elif approximate_num_orbits <= 200:
+            num_degree, denom_degree = 6, 6
+        else:
+            num_degree, denom_degree \
+                = (5 + int(np.log10(approximate_num_orbits)),
+                5 + int(np.log10(approximate_num_orbits)))
+        return num_degree, denom_degree
+
+    def check_if_first_derivative_is_not_strictly_monotonic(
+        self, x, y, tol=1.0, data_name=None):
+        """Check if the first derivative of data is not strictly monotonic.
         """
-        return interpolate(newX, oldX, oldY, allowExtrapolation, interpolator,
-                           spline_kwargs=self.spline_kwargs,
-                           check_kwargs=False)
+        dy_dx = np.gradient(y, x[1] - x[0])
+        is_not_strictly_monotonic = any(dy_dx[1:]/dy_dx[:-1] < tol)
+        # update history
+        if self.extra_kwargs["omega_gw_extrema_interpolation_method"] \
+            == "rational_fit":
+            if data_name in self.rational_fit_nonmonotonicity_history:
+                self.rational_fit_nonmonotonicity_history[data_name].update(
+                    {(self.special_interp_kwargs_for_extrema['num_degree'],
+                      self.special_interp_kwargs_for_extrema['num_degree']):
+                      is_not_strictly_monotonic})
+        return is_not_strictly_monotonic
+
+    def get_available_omega_gw_extrema_interpolation_methods(self):
+        """Return available omega_gw extrema interpolation methods."""
+        available_methods = {
+            "spline": get_interpolant,
+            "rational_fit": self.get_rational_fit_for_extrema,
+        }
+        return available_methods
 
-    def interp_extrema(self, extrema_type="pericenters"):
-        """Build interpolant through extrema.
+    def get_omega_gw_extrema_interpolant(self, extrema_type="pericenters"):
+        """Build interpolant through omega_gw extrema.
+
+        Calculating the eccentricity using the gravitational wave frequency
+        (`omega_gw`), requires building an interpolant through its extrema
+        (pericenters and apocenters).
+        
+        The method for constructing this interpolant can be specified by the
+        user via `omega_gw_extrema_interpolation_method` in
+        `extra_kwargs`. Currently, two options are supported: `"spline"` and
+        `"rational_fit"`.
+        
+        Based on the chosen method, an interpolation will be constructed for
+        `omega_gw` at the pericenter/apocenter points.
 
         parameters:
         -----------
@@ -1330,8 +1674,16 @@ def interp_extrema(self, extrema_type="pericenters"):
             raise Exception("extrema_type must be either "
                             "'pericenrers' or 'apocenters'.")
         if len(extrema) >= 2:
-            return self.get_interp(self.t[extrema],
-                                   self.omega_gw[extrema])
+            method = self.available_omega_gw_extrema_interpolation_methods[
+                self.omega_gw_extrema_interpolation_method]
+            if self.omega_gw_extrema_interpolation_method == "rational_fit":
+                return method(self.t[extrema], self.omega_gw[extrema],
+                extrema_type)
+            if self.omega_gw_extrema_interpolation_method == "spline":
+                return method(
+                    self.t[extrema], self.omega_gw[extrema], 
+                    spline_kwargs=self.special_interp_kwargs_for_extrema,
+                    check_kwargs=False)
         else:
             raise Exception(
                 f"Sufficient number of {extrema_type} are not found."
@@ -1634,11 +1986,11 @@ def measure_ecc(self, tref_in=None, fref_in=None):
             = self.check_num_extrema(apocenters, "apocenters")
 
         # If the eccentricity is too small for a method to find the extrema,
-        # and `return_zero_if_small_ecc_failure` is true, then we set the eccentricity and
-        # mean anomaly to zero and return them. In this case, the rest of the
-        # code in this function is not executed, and therefore, many variables
-        # that are needed for making diagnostic plots are not computed. Thus,
-        # in such cases, the diagnostic plots may not work.
+        # and `return_zero_if_small_ecc_failure` is true, then we set the
+        # eccentricity and mean anomaly to zero and return them. In this case,
+        # the rest of the code in this function is not executed, and therefore,
+        # many variables that are needed for making diagnostic plots are not
+        # computed. Thus, in such cases, the diagnostic plots may not work.
         if any([insufficient_pericenters_but_long_waveform,
                 insufficient_apocenters_but_long_waveform]) \
                 and self.return_zero_if_small_ecc_failure:
@@ -1672,14 +2024,11 @@ def measure_ecc(self, tref_in=None, fref_in=None):
             = self.check_extrema_separation(self.apocenters_location,
                                             "apocenters")
 
-        # Build the interpolants of omega_gw at the extrema
-        self.omega_gw_pericenters_interp = self.interp_extrema("pericenters")
-        self.omega_gw_apocenters_interp = self.interp_extrema("apocenters")
-
         self.t_pericenters = self.t[self.pericenters_location]
         self.t_apocenters = self.t[self.apocenters_location]
         self.tmax = min(self.t_pericenters[-1], self.t_apocenters[-1])
         self.tmin = max(self.t_pericenters[0], self.t_apocenters[0])
+        
         if self.domain == "frequency":
             # get the tref_in and fref_out from fref_in
             self.tref_in, self.fref_out \
@@ -1728,19 +2077,26 @@ def measure_ecc(self, tref_in=None, fref_in=None):
            or self.tref_out[-1] > self.t_pericenters[-1]:
             raise Exception("Reference time must be within two pericenters.")
 
+        # Build omega_gw extrema interpolants
+        self.omega_gw_pericenters_interp \
+            = self.get_omega_gw_extrema_interpolant("pericenters")
+        self.omega_gw_apocenters_interp \
+            = self.get_omega_gw_extrema_interpolant("apocenters")
+        # check monotonicity of the interpolants
+        self.check_omega_gw_extrema_interpolants()
         # compute eccentricity at self.tref_out
         self.eccentricity = self.compute_eccentricity(self.tref_out)
         # Compute mean anomaly at tref_out
         self.mean_anomaly = self.compute_mean_anomaly(self.tref_out)
 
+        # check if eccentricity is nonmonotonic
+        self.check_monotonicity_and_convexity()
+        
         # check if eccentricity is positive
         if any(self.eccentricity < 0):
             debug_message("Encountered negative eccentricity.",
                           self.debug_level, point_to_verbose_output=True)
 
-        # check if eccentricity is monotonic and convex
-        self.check_monotonicity_and_convexity()
-
         if self.debug_plots:
             # make a plot for diagnostics
             fig, axes = self.make_diagnostic_plots()
@@ -1749,6 +2105,29 @@ def measure_ecc(self, tref_in=None, fref_in=None):
         # return measured eccentricity, mean anomaly and reference time or
         # frequency where these are measured.
         return self.make_return_dict_for_eccentricity_and_mean_anomaly()
+    
+    def check_omega_gw_extrema_interpolants(self):
+        """
+        check monotonicity of the omega_gw extrema interpolants.
+        """
+        # Verify the monotonicity of the first derivative of the omega_gw
+        # interpolant with spline.
+        if self.extra_kwargs["omega_gw_extrema_interpolation_method"] \
+            == "spline":
+        # Check if the first derivative of omega_gw at pericenters or
+        # apocenters is non-monotonic
+            if (self.check_if_first_derivative_is_not_strictly_monotonic(
+                    self.t_for_checks,
+                    self.omega_gw_pericenters_interp(self.t_for_checks)) or
+                self.check_if_first_derivative_is_not_strictly_monotonic(
+                    self.t_for_checks,
+                    self.omega_gw_apocenters_interp(self.t_for_checks))):
+                debug_message(
+                    "Nonmonotonic time derivative detected in the spline "
+                    "interpolant through extrema. Using rational fit by "
+                    "setting 'omega_gw_extrema_interpolation_method' to "
+                    "'rational_fit' may provide better result.",
+                    debug_level=self.debug_level, important=True)
 
     def set_eccentricity_and_mean_anomaly_to_zero(self):
         """Set eccentricity and mean_anomaly to zero."""
@@ -1881,8 +2260,9 @@ def derivative_of_eccentricity(self, t, n=1):
                 self.t_for_checks)
 
         if self.ecc_interp is None:
-            self.ecc_interp = self.get_interp(self.t_for_checks,
-                                              self.ecc_for_checks)
+            self.ecc_interp = get_interpolant(
+                self.t_for_checks, self.ecc_for_checks,
+                spline_kwargs=self.general_interp_kwargs, check_kwargs=False)
         # Get derivative of ecc(t) using spline
         return self.ecc_interp.derivative(n=n)(t)
 
@@ -2205,13 +2585,15 @@ def compute_res_amp_gw_and_res_omega_gw(self):
         # residual quantities can be computed. Above, we check that this
         # extrapolation does not happen before t_merger, which is where
         # eccentricity is normally measured.
-        self.amp_gw_zeroecc_interp = self.interp(
+        self.amp_gw_zeroecc_interp = interpolate(
             self.t, self.t_zeroecc_shifted, self.amp_gw_zeroecc,
-            allowExtrapolation=True)
+            allowExtrapolation=True, spline_kwargs=self.general_interp_kwargs,
+            check_kwargs=False)
         self.res_amp_gw = self.amp_gw - self.amp_gw_zeroecc_interp
-        self.omega_gw_zeroecc_interp = self.interp(
+        self.omega_gw_zeroecc_interp = interpolate(
             self.t, self.t_zeroecc_shifted, self.omega_gw_zeroecc,
-            allowExtrapolation=True)
+            allowExtrapolation=True, spline_kwargs=self.general_interp_kwargs,
+            check_kwargs=False)
         self.res_omega_gw = (self.omega_gw - self.omega_gw_zeroecc_interp)
 
     def get_t_average_for_orbit_averaged_omega_gw(self):
@@ -2342,21 +2724,21 @@ def compute_orbit_averaged_omega_gw_between_extrema(self, t):
             orbit_averaged_omega_gw,
             "omega_gw averaged [apocenter to apocenter] and "
             "[pericenter to pericenter]")
-        return self.interp(
-            t, self.t_for_orbit_averaged_omega_gw, orbit_averaged_omega_gw)
+        return interpolate(
+            t, self.t_for_orbit_averaged_omega_gw, orbit_averaged_omega_gw,
+            spline_kwargs=self.general_interp_kwargs, check_kwargs=False)
 
     def check_monotonicity_of_omega_gw_average(self,
                                               omega_gw_average,
-                                              description="omega_gw average"):
+                                              data_name="omega_gw average"):
         """Check that omega average is monotonically increasing.
 
         Parameters
         ----------
         omega_gw_average : array-like
             1d array of omega_gw averages to check for monotonicity.
-        description : str
-            String to describe what the the which omega_gw average we are
-            looking at.
+        data_name : str
+            String to describe which omega_gw average we are looking at.
         """
         idx_non_monotonic = np.where(
             np.diff(omega_gw_average) <= 0)[0]
@@ -2412,7 +2794,7 @@ def check_monotonicity_of_omega_gw_average(self,
                 fig.tight_layout()
                 figName = (
                     "./gwecc_"
-                    f"{self.method}_{description.replace(' ', '_')}.pdf")
+                    f"{self.method}_{data_name.replace(' ', '_')}.pdf")
                 # fig.savefig(figName)
                 self.save_debug_fig(fig, figName)
                 plt.close(fig)
@@ -2422,7 +2804,7 @@ def check_monotonicity_of_omega_gw_average(self,
                              "for diagnostic plot use `debug_plots=True` in "
                              "extra_kwargs")
             raise Exception(
-                f"{description} are non-monotonic.\n"
+                f"{data_name} are non-monotonic.\n"
                 f"First non-monotonicity occurs at peak number {first_idx},"
                 f" where omega_gw drops from {omega_gw_average[first_idx]} to"
                 f" {omega_gw_average[first_idx+1]}, a decrease by"
@@ -2447,8 +2829,9 @@ def compute_mean_of_extrema_interpolants(self, t):
 
     def compute_omega_gw_zeroecc(self, t):
         """Find omega_gw from zeroecc data."""
-        return self.interp(
-            t, self.t_zeroecc_shifted, self.omega_gw_zeroecc)
+        return interpolate(
+            t, self.t_zeroecc_shifted, self.omega_gw_zeroecc,
+            spline_kwargs=self.general_interp_kwargs, check_kwargs=False)
 
     def get_available_omega_gw_averaging_methods(self):
         """Return available omega_gw averaging methods."""
@@ -2592,16 +2975,19 @@ def compute_tref_in_and_fref_out_from_fref_in(self, fref_in):
             # of omega_gw_average.
             # We get omega_gw_average by evaluating the omega_gw_average(t)
             # on t, from tmin_for_fref to tmax_for_fref
-            self.t_for_omega_gw_average, self.omega_gw_average = self.get_omega_gw_average(method)
+            self.t_for_omega_gw_average, self.omega_gw_average \
+                = self.get_omega_gw_average(method)
 
             # check that omega_gw_average is monotonically increasing
             self.check_monotonicity_of_omega_gw_average(
                 self.omega_gw_average, "Interpolated omega_gw_average")
 
             # Get tref_in using interpolation
-            tref_in = self.interp(fref_out,
+            tref_in = interpolate(fref_out,
                                   self.omega_gw_average/(2 * np.pi),
-                                  self.t_for_omega_gw_average)
+                                  self.t_for_omega_gw_average,
+                                  spline_kwargs=self.general_interp_kwargs,
+                                  check_kwargs=False)
             # check if tref_in is monotonically increasing
             if any(np.diff(tref_in) <= 0):
                 debug_message(f"tref_in from fref_in using method {method} is"
@@ -2640,7 +3026,8 @@ def get_fref_bounds(self, method=None):
             frequency.
         """
         if self.omega_gw_average is None:
-            self.t_for_omega_gw_average, self.omega_gw_average = self.get_omega_gw_average(method)
+            self.t_for_omega_gw_average, self.omega_gw_average \
+                = self.get_omega_gw_average(method)
         return [min(self.omega_gw_average)/2/np.pi,
                 max(self.omega_gw_average)/2/np.pi]
 
diff --git a/gw_eccentricity/gw_eccentricity.py b/gw_eccentricity/gw_eccentricity.py
index caaab318..79ca22c7 100644
--- a/gw_eccentricity/gw_eccentricity.py
+++ b/gw_eccentricity/gw_eccentricity.py
@@ -338,12 +338,30 @@ def measure_eccentricity(tref_in=None,
         Default value is "inertial".
 
     extra_kwargs: A dict of any extra kwargs to be passed. Allowed kwargs are:
-        spline_kwargs:
+        special_interp_kwargs_for_extrema: dict
+            A dictionary with a single key matching the current
+                `omega_gw_extrema_interpolation_method`. See under
+                `omega_gw_extrema_interpolation_method` for more details on
+                possible values of interpolation methods.
+
+                Example structure:
+                    {
+                        "method": {"param1": value1, "param2": value2},
+                    }
+
+                currently, the available options for "method" are: 
+
+                - "spline": default kwargs are set using
+                  `utils.get_default_spline_kwargs`
+                - "rational_fit": default kwargs are set using
+                  `utils.get_default_rational_fit_kwargs`
+
+        general_interp_kwargs: dict
             Dictionary of arguments to be passed to the spline interpolation
             routine (scipy.interpolate.InterpolatedUnivariateSpline) used to
-            compute quantities like omega_gw_pericenters(t) and
-            omega_gw_apocenters(t).
-            Defaults are set using utils.get_default_spline_kwargs
+            interpolate data other than the omega_gw extrema.
+            
+            Defaults are set using `utils.get_default_spline_kwargs`
 
         extrema_finding_kwargs:
             Dictionary of arguments to be passed to the extrema finder,
@@ -443,6 +461,29 @@ def measure_eccentricity(tref_in=None,
             mean anomaly to zero.
             USE THIS WITH CAUTION!
 
+        omega_gw_extrema_interpolation_method : str, default="rational_fit"
+            Specifies the method used to build the interpolations for 
+            `omega_gw_pericenters_interp(t)` or `omega_gw_apocenters_interp(t)`.
+            The available options are:
+
+            - `spline`: Uses `scipy.interpolate.InterpolatedUnivariateSpline`.
+                - Best suited for cleaner data, such as when waveform modes are generated
+                    using models like SEOB or TEOB.
+                - Faster to construct and evaluate.
+                - Since it fits through every data point, it may exhibit oscillatory 
+                    behavior, particularly near the merger.
+
+            - `rational_fit`: Uses `polyrat.StabilizedSKRationalApproximation`.
+                - Can handle both clean and noisy data, e.g., waveform modes
+                    from numerical simulations.
+                - Better monotonic behaviour, particularly near the merger.
+                - Significantly slower compared to the `spline` method. This is because
+                    finding optimal numerator and denominator degree needs several iterations
+                - Can suppress pathologies in the waveform that might be visible with 
+                    `spline`.
+
+            Default value: `"rational_fit"`.
+
     Returns
     -------
     A dictionary containing the following keys
diff --git a/gw_eccentricity/utils.py b/gw_eccentricity/utils.py
index e66b3647..3411b19a 100644
--- a/gw_eccentricity/utils.py
+++ b/gw_eccentricity/utils.py
@@ -1,6 +1,7 @@
 """Useful functions for the project."""
 import numpy as np
 import argparse
+from polyrat import StabilizedSKRationalApproximation
 from scipy.interpolate import InterpolatedUnivariateSpline
 from scipy.interpolate import PchipInterpolator
 import warnings
@@ -317,6 +318,42 @@ def get_interpolant(oldX,
     return interpolant
 
 
+def get_default_rational_fit_kwargs():
+    """Get default kwargs for rational fit."""
+    default_rational_fit_kwargs = {
+        "num_degree": None,
+        "denom_degree": None,
+        "norm": 2,
+        "maxiter": 20,
+        "verbose": False,
+        "xtol": 1e-07,
+    }
+    return default_rational_fit_kwargs
+
+
+def get_rational_fit(x, y, rational_fit_kwargs=None, check_kwargs=True):
+    """Get rational fit for the data set (x, y).
+
+    We use `polyrat.StabilizedSKRationalApproximation` to obtain
+    rational fit to the data.
+    """
+    if check_kwargs:
+        rational_fit_kwargs = check_kwargs_and_set_defaults(
+            rational_fit_kwargs,
+            get_default_rational_fit_kwargs(),
+            "rational_fit_kwargs",
+            "utils.get_default_rational_fit_kwargs"
+        )
+    # We use a rational approximation based on Stabilized Sanathanan-Koerner Iteration
+    # described in arXiv:2009.10803 and implemented in `polyrat.StabilizedSKRationalApproximation`.
+    rat = StabilizedSKRationalApproximation(**rational_fit_kwargs)
+    # The input x-axis data must be 2-dimensional
+    rat.fit(x.reshape(-1, 1), y)
+    # Using a lambda function to change the input 1-d data to
+    # float64 and reshape to 2-d as required by the fit.
+    return lambda t: rat(t.astype('float64').reshape(-1, 1))
+
+
 def debug_message(message, debug_level, important=True,
                   point_to_verbose_output=False):
     """Show message based on debug_level.
diff --git a/pytest.ini b/pytest.ini
new file mode 100644
index 00000000..8a199c10
--- /dev/null
+++ b/pytest.ini
@@ -0,0 +1,3 @@
+[pytest]  
+filterwarnings =  
+    ignore::DeprecationWarning:polyrat.*
\ No newline at end of file
diff --git a/setup.py b/setup.py
index ec1ceea7..76188fa2 100644
--- a/setup.py
+++ b/setup.py
@@ -31,7 +31,8 @@ def get_version():
         'h5py',
         'lalsuite',
         'sxs',
-        'scri'
+        'scri',
+        'polyrat'
     ],
     classifiers=[
         "Intended Audience :: Science/Research",
diff --git a/test/generate_regression_data.py b/test/generate_regression_data.py
index e17cbd0d..6b309835 100644
--- a/test/generate_regression_data.py
+++ b/test/generate_regression_data.py
@@ -17,10 +17,15 @@
     type=str,
     required=True,
     help="EccDefinition method to save the regression data for.")
+parser.add_argument(
+    "--interp_method",
+    type=str,
+    required=True,
+    help="omega_gw_extrema_interpolation_method to save the regression data for.")
 args = parser.parse_args()
 
 
-def generate_regression_data(method):
+def generate_regression_data(method, interp_method):
     """Generate data for regression test using a method."""
     # Load test waveform
     lal_kwargs = {"approximant": "EccentricTD",
@@ -42,7 +47,7 @@ def generate_regression_data(method):
         raise Exception(f"method {method} is not available. Must be one of "
                         f"{available_methods}")
 
-    extra_kwargs = {}
+    extra_kwargs = {"omega_gw_extrema_interpolation_method": interp_method}
     user_kwargs = extra_kwargs.copy()
     regression_data.update({"extra_kwargs": extra_kwargs})
     # Try evaluating at an array of times
@@ -56,9 +61,9 @@ def generate_regression_data(method):
     meanano_ref = gwecc_dict["mean_anomaly"]
     # We save the measured data 3 reference times
     n = len(tref_out)
-    dict_tref =  {"time": [tref_out[0], tref_out[n//4], tref_out[n//2]],
-                  "eccentricity": [ecc_ref[0], ecc_ref[n//4], ecc_ref[n//2]],
-                  "mean_anomaly": [meanano_ref[0], meanano_ref[n//4], meanano_ref[n//2]]}
+    dict_tref =  {"time": [tref_out[n//8], tref_out[n//4], tref_out[n//2]],
+                  "eccentricity": [ecc_ref[n//8], ecc_ref[n//4], ecc_ref[n//2]],
+                  "mean_anomaly": [meanano_ref[n//8], meanano_ref[n//4], meanano_ref[n//2]]}
         
     # Try evaluating at an array of frequencies
     gwecc_dict = measure_eccentricity(
@@ -70,18 +75,18 @@ def generate_regression_data(method):
     ecc_ref = gwecc_dict["eccentricity"]
     meanano_ref = gwecc_dict["mean_anomaly"]
     n = len(fref_out)
-    dict_fref = {"frequency": [fref_out[0], fref_out[n//4], fref_out[n//2]],
-                 "eccentricity": [ecc_ref[0], ecc_ref[n//4], ecc_ref[n//2]],
-                 "mean_anomaly": [meanano_ref[0], meanano_ref[n//4], meanano_ref[n//2]]}
+    dict_fref = {"frequency": [fref_out[n//8], fref_out[n//4], fref_out[n//2]],
+                 "eccentricity": [ecc_ref[n//8], ecc_ref[n//4], ecc_ref[n//2]],
+                 "mean_anomaly": [meanano_ref[n//8], meanano_ref[n//4], meanano_ref[n//2]]}
     regression_data.update({"tref": dict_tref,
                             "fref": dict_fref})
 
     if not os.path.exists(data_dir):
         os.mkdir(data_dir)
     # save to a json file
-    fl = open(f"{data_dir}/{method}_regression_data.json", "w")
+    fl = open(f"{data_dir}/{method}_{interp_method}_regression_data.json", "w")
     json.dump(regression_data, fl)
     fl.close()
 
 # generate regression data
-generate_regression_data(args.method)
+generate_regression_data(args.method, args.interp_method)
diff --git a/test/regression_data/AmplitudeFits_regression_data.json b/test/regression_data/AmplitudeFits_regression_data.json
deleted file mode 100644
index faeab019..00000000
--- a/test/regression_data/AmplitudeFits_regression_data.json
+++ /dev/null
@@ -1 +0,0 @@
-{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {}, "tref": {"time": [-14769.419163230703, -11192.719163230704, -7616.019163230703], "eccentricity": [0.14715861688347398, 0.1347862810647561, 0.11943002062852026], "mean_anomaly": [0.0, 0.9679656489638546, 5.4841245849525535]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12338650648955707, 0.11542822840300593, 0.10530937888307923], "mean_anomaly": [2.9766378492619907, 1.494348360501732, 2.830869757710353]}}
\ No newline at end of file
diff --git a/test/regression_data/AmplitudeFits_spline_regression_data.json b/test/regression_data/AmplitudeFits_spline_regression_data.json
new file mode 100644
index 00000000..002b5202
--- /dev/null
+++ b/test/regression_data/AmplitudeFits_spline_regression_data.json
@@ -0,0 +1 @@
+{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {"omega_gw_extrema_interpolation_method": "spline"}, "tref": {"time": [-12981.119163243457, -11192.719163243455, -7616.019163243456], "eccentricity": [0.14125503147618124, 0.1347862810667858, 0.11943002063131647], "mean_anomaly": [3.271739184693624, 0.9679656489638688, 5.4841245849525535]}, "fref": {"frequency": [0.00413802852038928, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.11926882656303317, 0.11542822840106648, 0.10530937888140646], "mean_anomaly": [5.832592108883006, 1.494348360501732, 2.830869757710083]}}
\ No newline at end of file
diff --git a/test/regression_data/Amplitude_regression_data.json b/test/regression_data/Amplitude_regression_data.json
deleted file mode 100644
index ec603f0b..00000000
--- a/test/regression_data/Amplitude_regression_data.json
+++ /dev/null
@@ -1 +0,0 @@
-{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {}, "tref": {"time": [-14768.619163230704, -11225.819163230703, -7683.019163230703], "eccentricity": [0.14714463101268582, 0.1348939080003978, 0.11972197672186036], "mean_anomaly": [0.0, 0.6424466289474253, 4.757738663295697]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12336087182171862, 0.11538856688498489, 0.10524430745756652], "mean_anomaly": [2.955433309853852, 1.4699121058753946, 2.7932679992191325]}}
\ No newline at end of file
diff --git a/test/regression_data/Amplitude_spline_regression_data.json b/test/regression_data/Amplitude_spline_regression_data.json
new file mode 100644
index 00000000..38be1630
--- /dev/null
+++ b/test/regression_data/Amplitude_spline_regression_data.json
@@ -0,0 +1 @@
+{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {"omega_gw_extrema_interpolation_method": "spline"}, "tref": {"time": [-12997.219163243457, -11225.819163243457, -7683.019163243456], "eccentricity": [0.1412963090542808, 0.13489390800343737, 0.11972197672456475], "mean_anomaly": [3.118093974362708, 0.642446628947404, 4.757738663295697]}, "fref": {"frequency": [0.00413802852038928, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.11923643778091675, 0.11538856688636778, 0.10524430745877145], "mean_anomaly": [5.8083670411936055, 1.4699121058751956, 2.7932679992189193]}}
\ No newline at end of file
diff --git a/test/regression_data/FrequencyFits_regression_data.json b/test/regression_data/FrequencyFits_regression_data.json
deleted file mode 100644
index 3a0e061e..00000000
--- a/test/regression_data/FrequencyFits_regression_data.json
+++ /dev/null
@@ -1 +0,0 @@
-{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {}, "tref": {"time": [-14769.519163230703, -11193.019163230703, -7616.419163230703], "eccentricity": [0.14715912645792362, 0.13478783801420957, 0.11943256094553112], "mean_anomaly": [0.0, 0.9671579875661607, 5.4821375349876575]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12338317931464393, 0.11542730894258924, 0.1053075001590753], "mean_anomaly": [2.988070979817415, 1.5003684452927502, 2.8390552703477496]}}
\ No newline at end of file
diff --git a/test/regression_data/FrequencyFits_spline_regression_data.json b/test/regression_data/FrequencyFits_spline_regression_data.json
new file mode 100644
index 00000000..d532f681
--- /dev/null
+++ b/test/regression_data/FrequencyFits_spline_regression_data.json
@@ -0,0 +1 @@
+{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {"omega_gw_extrema_interpolation_method": "spline"}, "tref": {"time": [-12981.319163243457, -11193.019163243456, -7616.419163243456], "eccentricity": [0.1412559990968395, 0.13478783801299066, 0.11943256094180654], "mean_anomaly": [3.271739184693624, 0.9671579875661678, 5.4821375349876575]}, "fref": {"frequency": [0.00413802852038928, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.11926984282942188, 0.1154273089414749, 0.10530750015791723], "mean_anomaly": [5.833828675147657, 1.5003684452927075, 2.839055270348034]}}
\ No newline at end of file
diff --git a/test/regression_data/Frequency_regression_data.json b/test/regression_data/Frequency_regression_data.json
deleted file mode 100644
index 8eab7988..00000000
--- a/test/regression_data/Frequency_regression_data.json
+++ /dev/null
@@ -1 +0,0 @@
-{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {}, "tref": {"time": [-14768.719163230704, -11226.819163230703, -7684.919163230703], "eccentricity": [0.14714131326030366, 0.13489405451806347, 0.11972728556646606], "mean_anomaly": [0.0, 0.634809819093654, 4.7407656825570825]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12335278321610421, 0.11538428023363356, 0.10524044651268605], "mean_anomaly": [2.9664207187406504, 1.4731955044661476, 2.7980616912547234]}}
\ No newline at end of file
diff --git a/test/regression_data/Frequency_spline_regression_data.json b/test/regression_data/Frequency_spline_regression_data.json
new file mode 100644
index 00000000..9025b76f
--- /dev/null
+++ b/test/regression_data/Frequency_spline_regression_data.json
@@ -0,0 +1 @@
+{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {"omega_gw_extrema_interpolation_method": "spline"}, "tref": {"time": [-12997.819163243457, -11226.819163243457, -7684.919163243456], "eccentricity": [0.14129511303046816, 0.13489405452014136, 0.11972728556745593], "mean_anomaly": [3.114478792943153, 0.634809819093654, 4.7407656825570825]}, "fref": {"frequency": [0.00413802852038928, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.11923247746634258, 0.11538428023179614, 0.10524044651877096], "mean_anomaly": [5.811428452219403, 1.473195504466105, 2.798061691253885]}}
\ No newline at end of file
diff --git a/test/regression_data/ResidualAmplitude_regression_data.json b/test/regression_data/ResidualAmplitude_regression_data.json
deleted file mode 100644
index 75ef4a34..00000000
--- a/test/regression_data/ResidualAmplitude_regression_data.json
+++ /dev/null
@@ -1 +0,0 @@
-{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {}, "tref": {"time": [-14769.519163230703, -11193.019163230703, -7616.519163230703], "eccentricity": [0.14715912350724414, 0.134787764326868, 0.11943299066186941], "mean_anomaly": [0.0, 0.9662032412803185, 5.4800155541341695]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.123383193388651, 0.11542720219146674, 0.10530996061674747], "mean_anomaly": [2.9867814403803123, 1.499264733634135, 2.8323037712234154]}}
\ No newline at end of file
diff --git a/test/regression_data/ResidualAmplitude_spline_regression_data.json b/test/regression_data/ResidualAmplitude_spline_regression_data.json
new file mode 100644
index 00000000..9caccb28
--- /dev/null
+++ b/test/regression_data/ResidualAmplitude_spline_regression_data.json
@@ -0,0 +1 @@
+{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {"omega_gw_extrema_interpolation_method": "spline"}, "tref": {"time": [-12981.319163243457, -11193.019163243456, -7616.519163243456], "eccentricity": [0.14125596891213044, 0.13478776432643125, 0.11943299066201574], "mean_anomaly": [3.2708353893387336, 0.9662032412803043, 5.4800155541341695]}, "fref": {"frequency": [0.00413802852038928, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.11926975624561775, 0.1154272021945365, 0.10530996061768227], "mean_anomaly": [5.832818126814786, 1.4992647336341633, 2.8323037712231454]}}
\ No newline at end of file
diff --git a/test/regression_data/ResidualFrequency_regression_data.json b/test/regression_data/ResidualFrequency_regression_data.json
deleted file mode 100644
index b6df16de..00000000
--- a/test/regression_data/ResidualFrequency_regression_data.json
+++ /dev/null
@@ -1 +0,0 @@
-{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {}, "tref": {"time": [-14769.519163230703, -11193.019163230703, -7616.519163230703], "eccentricity": [0.14715912219591287, 0.13478782247732968, 0.11943303944129302], "mean_anomaly": [0.0, 0.9670110477321288, 5.481076544560906]}, "fref": {"frequency": [0.0039788735772973835, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.12338313657709132, 0.11542727032631284, 0.10531052353292425], "mean_anomaly": [2.9880704763720587, 1.500384961620071, 2.8327681616288416]}}
\ No newline at end of file
diff --git a/test/regression_data/ResidualFrequency_spline_regression_data.json b/test/regression_data/ResidualFrequency_spline_regression_data.json
new file mode 100644
index 00000000..1dc9abd0
--- /dev/null
+++ b/test/regression_data/ResidualFrequency_spline_regression_data.json
@@ -0,0 +1 @@
+{"waveform_kwargs": {"approximant": "EccentricTD", "q": 1.0, "chi1": [0.0, 0.0, 0.0], "chi2": [0.0, 0.0, 0.0], "Momega0": 0.01, "ecc": 0.1, "mean_ano": 0, "include_zero_ecc": true}, "extra_kwargs": {"omega_gw_extrema_interpolation_method": "spline"}, "tref": {"time": [-12981.319163243457, -11193.019163243456, -7616.519163243456], "eccentricity": [0.14125598209018864, 0.1347878224760185, 0.1194330394418408], "mean_anomaly": [3.271739184693624, 0.967011047732143, 5.481076544560906]}, "fref": {"frequency": [0.00413802852038928, 0.004297183463481175, 0.0047746482927568615], "eccentricity": [0.11926981413331783, 0.11542727032728273, 0.10531052353279924], "mean_anomaly": [5.833866315156342, 1.5003849616201421, 2.8327681616289127]}}
\ No newline at end of file
diff --git a/test/test_mks_vs_dimless_units.py b/test/test_mks_vs_dimless_units.py
index c48a4fd6..769abf2e 100644
--- a/test/test_mks_vs_dimless_units.py
+++ b/test/test_mks_vs_dimless_units.py
@@ -38,169 +38,184 @@ def test_mks_vs_dimless_units():
                        "D": 1})
     dataDictMKS = load_data.load_waveform(**lal_kwargs)
 
+    # use different omega_gw_extrema_interpolation methods
+    # TODO: Add rational_fit to the list
+    omega_gw_extrema_interpolation_methods = ["spline"]
+
     # List of all available methods
     available_methods = gw_eccentricity.get_available_methods()
     for method in available_methods:
-        # Try evaluating at a single dimless time t = -12000
-        idx = np.argmin(np.abs(dataDict["t"] - (-12000)))
-        gwecc_dict = measure_eccentricity(
-            tref_in=dataDict["t"][idx],
-            method=method,
-            dataDict=dataDict)
-        tref_out = gwecc_dict["tref_out"]
-        ecc_ref = gwecc_dict["eccentricity"]
-        meanano_ref = gwecc_dict["mean_anomaly"]
-        gwecc_object = gwecc_dict["gwecc_object"]
-        # Try evaluating at a single MKS time
-        gwecc_dict_MKS = measure_eccentricity(
-            tref_in=dataDictMKS["t"][idx],
-            method=method,
-            dataDict=dataDictMKS)
-        tref_out_MKS = gwecc_dict_MKS["tref_out"]
-        ecc_ref_MKS = gwecc_dict_MKS["eccentricity"]
-        meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"]
-        gwecc_object_MKS = gwecc_dict_MKS["gwecc_object"]
-        # Check that the tref_out is the same as rescaled (to dimless) tref_out_MKS
-        sec_to_dimless = 1/lal_kwargs["M"]/MTSUN_SI
-        np.testing.assert_allclose(
-            [tref_out],
-            [tref_out_MKS * sec_to_dimless],
-            err_msg=("tref_out at a single dimensionless and MKS"
-                     " time are inconsistent.\n"
-                     "x = Dimensionless, y = MKS converted to dimless"))
-        # Check if the measured ecc an mean ano are the same from the two units
-        np.testing.assert_allclose(
-            [ecc_ref],
-            [ecc_ref_MKS],
-            err_msg=("Eccentricity at a single dimensionless and MKS"
-                     " time gives different results.\n"
-                     "x = Dimensionless, y = MKS"))
-        np.testing.assert_allclose(
-            [meanano_ref],
-            [meanano_ref_MKS],
-            err_msg=("Mean anomaly at a single dimensionless and MKS"
-                     " time gives different results.\n"
-                     "x = Dimensionless, y = MKS"))
+        for omega_gw_extrema_interpolation_method in omega_gw_extrema_interpolation_methods:
+            extra_kwargs = {
+                "omega_gw_extrema_interpolation_method": omega_gw_extrema_interpolation_method}
+            # Try evaluating at a single dimless time t = -12000
+            idx = np.argmin(np.abs(dataDict["t"] - (-12000)))
+            gwecc_dict = measure_eccentricity(
+                tref_in=dataDict["t"][idx],
+                method=method,
+                dataDict=dataDict,
+                extra_kwargs=extra_kwargs)
+            tref_out = gwecc_dict["tref_out"]
+            ecc_ref = gwecc_dict["eccentricity"]
+            meanano_ref = gwecc_dict["mean_anomaly"]
+            gwecc_object = gwecc_dict["gwecc_object"]
+            # Try evaluating at a single MKS time
+            gwecc_dict_MKS = measure_eccentricity(
+                tref_in=dataDictMKS["t"][idx],
+                method=method,
+                dataDict=dataDictMKS,
+                extra_kwargs=extra_kwargs)
+            tref_out_MKS = gwecc_dict_MKS["tref_out"]
+            ecc_ref_MKS = gwecc_dict_MKS["eccentricity"]
+            meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"]
+            gwecc_object_MKS = gwecc_dict_MKS["gwecc_object"]
+            # Check that the tref_out is the same as rescaled (to dimless) tref_out_MKS
+            sec_to_dimless = 1/lal_kwargs["M"]/MTSUN_SI
+            np.testing.assert_allclose(
+                [tref_out],
+                [tref_out_MKS * sec_to_dimless],
+                err_msg=("tref_out at a single dimensionless and MKS"
+                         " time are inconsistent.\n"
+                         "x = Dimensionless, y = MKS converted to dimless"))
+            # Check if the measured ecc an mean ano are the same from the two units
+            np.testing.assert_allclose(
+                [ecc_ref],
+                [ecc_ref_MKS],
+                err_msg=("Eccentricity at a single dimensionless and MKS"
+                         " time gives different results.\n"
+                         "x = Dimensionless, y = MKS"))
+            np.testing.assert_allclose(
+                [meanano_ref],
+                [meanano_ref_MKS],
+                err_msg=("Mean anomaly at a single dimensionless and MKS"
+                         " time gives different results.\n"
+                         "x = Dimensionless, y = MKS"))
 
-        # Try evaluating at an array of dimless times
-        idx_start = np.argmin(np.abs(dataDict["t"] - (-12000)))
-        idx_end = np.argmin(np.abs(dataDict["t"] - (-5000)))
-        gwecc_dict = measure_eccentricity(
-            tref_in=dataDict["t"][idx_start: idx_end],
-            method=method,
-            dataDict=dataDict)
-        tref_out = gwecc_dict["tref_out"]
-        ecc_ref = gwecc_dict["eccentricity"]
-        meanano_ref = gwecc_dict["mean_anomaly"]
-        # Try evaluating at an array of MKS times
-        gwecc_dict_MKS = measure_eccentricity(
-            tref_in=dataDictMKS["t"][idx_start: idx_end],
-            method=method,
-            dataDict=dataDictMKS)
-        tref_out_MKS = gwecc_dict_MKS["tref_out"]
-        ecc_ref_MKS = gwecc_dict_MKS["eccentricity"]
-        meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"]
-        # check that the tref_out times are consistent
-        np.testing.assert_allclose(
-            [tref_out],
-            [tref_out_MKS * sec_to_dimless],
-            err_msg=("tref_out array for dimensionless and MKS"
-                     " tref_in are inconsistent.\n"
-                     "x = Dimensionless, y = MKS converted to dimless"))
-        # Check if the measured ecc an mean ano are the same from the two units
-        np.testing.assert_allclose(
-            ecc_ref,
-            ecc_ref_MKS,
-            err_msg=("Eccentricity at dimensionless and MKS array of"
-                     " times are different\n."
-                     "x = Dimensionless, y = MKS"))
-        # using unwrapped mean anomaly since 0 and 2pi should be treated as
-        # the same
-        np.testing.assert_allclose(
-            np.unwrap(meanano_ref),
-            np.unwrap(meanano_ref_MKS),
-            err_msg=("Mean anomaly at dimensionless and MKS array of"
-                     " times are different.\n"
-                     "x = Dimensionless, y = MKS"))
+            # Try evaluating at an array of dimless times
+            idx_start = np.argmin(np.abs(dataDict["t"] - (-12000)))
+            idx_end = np.argmin(np.abs(dataDict["t"] - (-5000)))
+            gwecc_dict = measure_eccentricity(
+                tref_in=dataDict["t"][idx_start: idx_end],
+                method=method,
+                dataDict=dataDict,
+                extra_kwargs=extra_kwargs)
+            tref_out = gwecc_dict["tref_out"]
+            ecc_ref = gwecc_dict["eccentricity"]
+            meanano_ref = gwecc_dict["mean_anomaly"]
+            # Try evaluating at an array of MKS times
+            gwecc_dict_MKS = measure_eccentricity(
+                tref_in=dataDictMKS["t"][idx_start: idx_end],
+                method=method,
+                dataDict=dataDictMKS,
+                extra_kwargs=extra_kwargs)
+            tref_out_MKS = gwecc_dict_MKS["tref_out"]
+            ecc_ref_MKS = gwecc_dict_MKS["eccentricity"]
+            meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"]
+            # check that the tref_out times are consistent
+            np.testing.assert_allclose(
+                [tref_out],
+                [tref_out_MKS * sec_to_dimless],
+                err_msg=("tref_out array for dimensionless and MKS"
+                         " tref_in are inconsistent.\n"
+                         "x = Dimensionless, y = MKS converted to dimless"))
+            # Check if the measured ecc an mean ano are the same from the two units
+            np.testing.assert_allclose(
+                ecc_ref,
+                ecc_ref_MKS,
+                err_msg=("Eccentricity at dimensionless and MKS array of"
+                         " times are different\n."
+                         "x = Dimensionless, y = MKS"))
+            # using unwrapped mean anomaly since 0 and 2pi should be treated as
+            # the same
+            np.testing.assert_allclose(
+                np.unwrap(meanano_ref),
+                np.unwrap(meanano_ref_MKS),
+                err_msg=("Mean anomaly at dimensionless and MKS array of"
+                         " times are different.\n"
+                         "x = Dimensionless, y = MKS"))
 
-        # Try evaluating at single dimensionless frequency
-        fref_in = gwecc_object.compute_orbit_averaged_omega_gw_between_extrema(
-            dataDict["t"][idx]) / (2 * np.pi)
-        gwecc_dict = measure_eccentricity(
-            fref_in=fref_in,
-            method=method,
-            dataDict=dataDict)
-        fref_out = gwecc_dict["fref_out"]
-        ecc_ref = gwecc_dict["eccentricity"]
-        meanano_ref = gwecc_dict["mean_anomaly"]
-        # Try evaluating at single MKS frequency
-        fref_in = gwecc_object_MKS.compute_orbit_averaged_omega_gw_between_extrema(
-            dataDictMKS["t"][idx]) / (2 * np.pi)
-        gwecc_dict_MKS = measure_eccentricity(
-            fref_in=fref_in,
-            method=method,
-            dataDict=dataDictMKS)
-        fref_out_MKS = gwecc_dict_MKS["fref_out"]
-        ecc_ref_MKS = gwecc_dict_MKS["eccentricity"]
-        meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"]
-        # Check the fref_out frequencies are consistent
-        np.testing.assert_allclose(
-            [fref_out],
-            [fref_out_MKS / sec_to_dimless],
-            err_msg=("fref_out for a single dimensionless and MKS"
-                     " fref_in are inconsistent.\n"
-                     "x = Dimensionless, y = MKS converted to dimless"))
-        # Check if the measured ecc an mean ano are the same from the two units
-        np.testing.assert_allclose(
-            [ecc_ref],
-            [ecc_ref_MKS],
-            err_msg=("Eccentricity at a single dimensionless and MKS"
-                     " frequency gives different results.\n"
-                     "x = Dimensionless, y = MKS"))
-        np.testing.assert_allclose(
-            [meanano_ref],
-            [meanano_ref_MKS],
-            err_msg=("Mean anomaly at a single dimensionless and MKS"
-                     " frequency gives different results.\n"
-                     "x = Dimensionless, y = MKS"))
+            # Try evaluating at single dimensionless frequency
+            fref_in = gwecc_object.compute_orbit_averaged_omega_gw_between_extrema(
+                dataDict["t"][idx]) / (2 * np.pi)
+            gwecc_dict = measure_eccentricity(
+                fref_in=fref_in,
+                method=method,
+                dataDict=dataDict,
+                extra_kwargs=extra_kwargs)
+            fref_out = gwecc_dict["fref_out"]
+            ecc_ref = gwecc_dict["eccentricity"]
+            meanano_ref = gwecc_dict["mean_anomaly"]
+            # Try evaluating at single MKS frequency
+            fref_in = gwecc_object_MKS.compute_orbit_averaged_omega_gw_between_extrema(
+                dataDictMKS["t"][idx]) / (2 * np.pi)
+            gwecc_dict_MKS = measure_eccentricity(
+                fref_in=fref_in,
+                method=method,
+                dataDict=dataDictMKS,
+                extra_kwargs=extra_kwargs)
+            fref_out_MKS = gwecc_dict_MKS["fref_out"]
+            ecc_ref_MKS = gwecc_dict_MKS["eccentricity"]
+            meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"]
+            # Check the fref_out frequencies are consistent
+            np.testing.assert_allclose(
+                [fref_out],
+                [fref_out_MKS / sec_to_dimless],
+                err_msg=("fref_out for a single dimensionless and MKS"
+                         " fref_in are inconsistent.\n"
+                         "x = Dimensionless, y = MKS converted to dimless"))
+            # Check if the measured ecc an mean ano are the same from the two units
+            np.testing.assert_allclose(
+                [ecc_ref],
+                [ecc_ref_MKS],
+                err_msg=("Eccentricity at a single dimensionless and MKS"
+                         " frequency gives different results.\n"
+                         "x = Dimensionless, y = MKS"))
+            np.testing.assert_allclose(
+                [meanano_ref],
+                [meanano_ref_MKS],
+                err_msg=("Mean anomaly at a single dimensionless and MKS"
+                         " frequency gives different results.\n"
+                         "x = Dimensionless, y = MKS"))
 
-        # Try evaluating at an array of dimensionless frequencies
-        fref_in = gwecc_object.compute_orbit_averaged_omega_gw_between_extrema(
-            dataDict["t"][idx: idx+500]) / (2 * np.pi)
-        gwecc_dict = measure_eccentricity(
-            fref_in=fref_in,
-            method=method,
-            dataDict=dataDict)
-        fref_out = gwecc_dict["fref_out"]
-        ecc_ref = gwecc_dict["eccentricity"]
-        meanano_ref = gwecc_dict["mean_anomaly"]
-        # Try evaluating at an array of MKS frequencies
-        fref_in = gwecc_object_MKS.compute_orbit_averaged_omega_gw_between_extrema(
-            dataDictMKS["t"][idx: idx+500]) / (2 * np.pi)
-        gwecc_dict_MKS = measure_eccentricity(
-            fref_in=fref_in,
-            method=method,
-            dataDict=dataDictMKS)
-        fref_out_MKS = gwecc_dict_MKS["fref_out"]
-        ecc_ref_MKS = gwecc_dict_MKS["eccentricity"]
-        meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"]
-        # Check fref_out
-        np.testing.assert_allclose(
-            [fref_out],
-            [fref_out_MKS / sec_to_dimless],
-            err_msg=("fref_out for an array of dimensionless and MKS"
-                     " fref_in are inconsistent.\n"
-                     "x = Dimensionless, y = MKS converted to dimless"))
-        # Check if the measured ecc an mean ano are the same from the two units
-        np.testing.assert_allclose(
-            ecc_ref,
-            ecc_ref_MKS,
-            err_msg=("Eccentricity at dimensionless and MKS array of"
-                     " frequencies are different.\n"
-                     "x = Dimensionless, y = MKS"))
-        np.testing.assert_allclose(
-            np.unwrap(meanano_ref),
-            np.unwrap(meanano_ref_MKS),
-            err_msg=("Mean anomaly at dimensionless and MKS array of"
-                     " frequencies are different.\n"
-                     "x = Dimensionless, y = MKS"))
+            # Try evaluating at an array of dimensionless frequencies
+            fref_in = gwecc_object.compute_orbit_averaged_omega_gw_between_extrema(
+                dataDict["t"][idx: idx+500]) / (2 * np.pi)
+            gwecc_dict = measure_eccentricity(
+                fref_in=fref_in,
+                method=method,
+                dataDict=dataDict,
+                extra_kwargs=extra_kwargs)
+            fref_out = gwecc_dict["fref_out"]
+            ecc_ref = gwecc_dict["eccentricity"]
+            meanano_ref = gwecc_dict["mean_anomaly"]
+            # Try evaluating at an array of MKS frequencies
+            fref_in = gwecc_object_MKS.compute_orbit_averaged_omega_gw_between_extrema(
+                dataDictMKS["t"][idx: idx+500]) / (2 * np.pi)
+            gwecc_dict_MKS = measure_eccentricity(
+                fref_in=fref_in,
+                method=method,
+                dataDict=dataDictMKS,
+                extra_kwargs=extra_kwargs)
+            fref_out_MKS = gwecc_dict_MKS["fref_out"]
+            ecc_ref_MKS = gwecc_dict_MKS["eccentricity"]
+            meanano_ref_MKS = gwecc_dict_MKS["mean_anomaly"]
+            # Check fref_out
+            np.testing.assert_allclose(
+                [fref_out],
+                [fref_out_MKS / sec_to_dimless],
+                err_msg=("fref_out for an array of dimensionless and MKS"
+                         " fref_in are inconsistent.\n"
+                         "x = Dimensionless, y = MKS converted to dimless"))
+            # Check if the measured ecc an mean ano are the same from the two units
+            np.testing.assert_allclose(
+                ecc_ref,
+                ecc_ref_MKS,
+                err_msg=("Eccentricity at dimensionless and MKS array of"
+                         " frequencies are different.\n"
+                         "x = Dimensionless, y = MKS"))
+            np.testing.assert_allclose(
+                np.unwrap(meanano_ref),
+                np.unwrap(meanano_ref_MKS),
+                err_msg=("Mean anomaly at dimensionless and MKS array of"
+                         " frequencies are different.\n"
+                         "x = Dimensionless, y = MKS"))
diff --git a/test/test_regression.py b/test/test_regression.py
index ca4e3345..2de08489 100644
--- a/test/test_regression.py
+++ b/test/test_regression.py
@@ -4,6 +4,7 @@
 import gw_eccentricity
 from gw_eccentricity import load_data
 from gw_eccentricity import measure_eccentricity
+import glob
 
 # locally it passs without problem but on github action, we need to set this tolerance
 atol = 1e-9
@@ -14,55 +15,56 @@ def test_regression():
     available_methods = gw_eccentricity.get_available_methods()
     for method in available_methods:
         # Load regression data
-        regression_data_file = f"test/regression_data/{method}_regression_data.json"
-        # Load the regression data
-        fl = open(regression_data_file, "r")
-        regression_data = json.load(fl)
-        fl.close()
-        # waveform kwargs
-        lal_kwargs = regression_data["waveform_kwargs"]
-        # load waveform data
-        dataDict = load_data.load_waveform(**lal_kwargs)
-        extra_kwargs = regression_data["extra_kwargs"]
+        regression_data_files = glob.glob(f"test/regression_data/{method}_*_regression_data.json")
+        for regression_data_file in regression_data_files:
+            # Load the regression data
+            fl = open(regression_data_file, "r")
+            regression_data = json.load(fl)
+            fl.close()
+            # waveform kwargs
+            lal_kwargs = regression_data["waveform_kwargs"]
+            # load waveform data
+            dataDict = load_data.load_waveform(**lal_kwargs)
+            extra_kwargs = regression_data["extra_kwargs"]
 
-        # Try evaluating at times where regression data are saved
-        regression_data_at_tref = regression_data["tref"]
-        tref_in = regression_data_at_tref["time"]
-        gwecc_dict = measure_eccentricity(
-            tref_in=tref_in,
-            method=method,
-            dataDict=dataDict,
-            extra_kwargs=extra_kwargs)
-        tref_out = gwecc_dict["tref_out"]
-        ecc_ref = gwecc_dict["eccentricity"]
-        meanano_ref = gwecc_dict["mean_anomaly"]
-        # Compare the measured data with the saved data
-        np.testing.assert_allclose(
-            ecc_ref, regression_data_at_tref["eccentricity"],
-            atol=atol,
-            err_msg="measured and saved eccentricity at saved times do not match.")
-        np.testing.assert_allclose(
-            meanano_ref, regression_data_at_tref["mean_anomaly"],
-            atol=atol,
-            err_msg="measured and saved mean anomaly at saved times do not match.")
+            # Try evaluating at times where regression data are saved
+            regression_data_at_tref = regression_data["tref"]
+            tref_in = regression_data_at_tref["time"]
+            gwecc_dict = measure_eccentricity(
+                tref_in=tref_in,
+                method=method,
+                dataDict=dataDict,
+                extra_kwargs=extra_kwargs)
+            tref_out = gwecc_dict["tref_out"]
+            ecc_ref = gwecc_dict["eccentricity"]
+            meanano_ref = gwecc_dict["mean_anomaly"]
+            # Compare the measured data with the saved data
+            np.testing.assert_allclose(
+                ecc_ref, regression_data_at_tref["eccentricity"],
+                atol=atol,
+                err_msg="measured and saved eccentricity at saved times do not match.")
+            np.testing.assert_allclose(
+                meanano_ref, regression_data_at_tref["mean_anomaly"],
+                atol=atol,
+                err_msg="measured and saved mean anomaly at saved times do not match.")
 
-        # Try evaluating at frequencies where regression data are saved
-        regression_data_at_fref = regression_data["fref"]
-        fref_in = regression_data_at_fref["frequency"]
-        gwecc_dict = measure_eccentricity(
-            fref_in=fref_in,
-            method=method,
-            dataDict=dataDict,
-            extra_kwargs=extra_kwargs)
-        fref_out = gwecc_dict["fref_out"]
-        ecc_ref = gwecc_dict["eccentricity"]
-        meanano_ref = gwecc_dict["mean_anomaly"]
-        # Compare the measured data with the saved data
-        np.testing.assert_allclose(
-            ecc_ref, regression_data_at_fref["eccentricity"],
-            atol=atol,
-            err_msg="measured and saved eccentricity at saved frequencies do not match.")
-        np.testing.assert_allclose(
-            meanano_ref, regression_data_at_fref["mean_anomaly"],
-            atol=atol,
-            err_msg="measured and saved mean anomaly at saved frequencies do not match.")
+            # Try evaluating at frequencies where regression data are saved
+            regression_data_at_fref = regression_data["fref"]
+            fref_in = regression_data_at_fref["frequency"]
+            gwecc_dict = measure_eccentricity(
+                fref_in=fref_in,
+                method=method,
+                dataDict=dataDict,
+                extra_kwargs=extra_kwargs)
+            fref_out = gwecc_dict["fref_out"]
+            ecc_ref = gwecc_dict["eccentricity"]
+            meanano_ref = gwecc_dict["mean_anomaly"]
+            # Compare the measured data with the saved data
+            np.testing.assert_allclose(
+                ecc_ref, regression_data_at_fref["eccentricity"],
+                atol=atol,
+                err_msg="measured and saved eccentricity at saved frequencies do not match.")
+            np.testing.assert_allclose(
+                meanano_ref, regression_data_at_fref["mean_anomaly"],
+                atol=atol,
+                err_msg="measured and saved mean anomaly at saved frequencies do not match.")