diff --git a/.github/workflows/ci_tests.yml b/.github/workflows/ci_tests.yml
index 5e24a17..ea369de 100644
--- a/.github/workflows/ci_tests.yml
+++ b/.github/workflows/ci_tests.yml
@@ -27,40 +27,25 @@ jobs:
             python: 3.x
             toxenv: codestyle
 
-          - name: Python 3.7 with minimal dependencies
+          - name: Python 3.9 with minimal dependencies
             os: ubuntu-latest
-            python: 3.7
-            toxenv: py37-test
+            python: 3.9
+            toxenv: py39-test
 
-          - name: Python 3.8 with all optional dependencies and coverage checking
+          - name: Python 3.11 with all optional dependencies and coverage checking
             os: ubuntu-latest
-            python: 3.8
-            toxenv: py38-test-alldeps-cov
+            python: 3.11
+            toxenv: py311-test-alldeps-cov
 
-          - name: macOS - Python 3.8 with all optional dependencies
+          - name: macOS - Python 3.11 with all optional dependencies
             os: macos-latest
-            python: 3.8
-            toxenv: py38-test-alldeps
+            python: 3.11
+            toxenv: py311-test-alldeps
 
-          - name: Windows - Python 3.8 with all optional dependencies
+          - name: Windows - Python 3.11 with all optional dependencies
             os: windows-latest
-            python: 3.8
-            toxenv: py38-test-alldeps
-
-          # - name: Python 3.7 with oldest supported version of all dependencies
-          #   os: ubuntu-16.04
-          #   python: 3.7
-          #   toxenv: py37-test-oldestdeps
-
-          # - name: Python 3.8 with latest dev versions of key dependencies
-          #   os: ubuntu-latest
-          #   python: 3.8
-          #   toxenv: py38-test-devdeps
-
-          # - name: Test building of Sphinx docs
-          #   os: ubuntu-latest
-          #   python: 3.x
-          #   toxenv: build_docs
+            python: 3.11
+            toxenv: py311-test-alldeps
 
     steps:
     - name: Checkout code
diff --git a/ell1fit/create_parfile.py b/ell1fit/create_parfile.py
index 949c20e..f348139 100644
--- a/ell1fit/create_parfile.py
+++ b/ell1fit/create_parfile.py
@@ -40,10 +40,13 @@ def update_model(model, value_dict):
         value = mean * factor + initial
         err = max(neg, pos) * factor
         if par == "Phase":
-            new_model.TZRMJD.value = -value / new_model.F0.value / 86400 + PEPOCH
-            new_model.TZRMJD.uncertainty_value = err / new_model.F0.value / 86400
-            new_model.TZRMJD.frozen = False
-            # new_model.TZRMJD.value =  PEPOCH
+            try:
+                new_model.TZRMJD.value = -value / new_model.F0.value / 86400 + PEPOCH
+                new_model.TZRMJD.uncertainty_value = err / new_model.F0.value / 86400
+                new_model.TZRMJD.frozen = False
+                # new_model.TZRMJD.value =  PEPOCH
+            except ValueError:
+                pass
             continue
         if par == "PB":
             value /= 86400
@@ -56,7 +59,12 @@ def update_model(model, value_dict):
         getattr(new_model, par).uncertainty_value = err
         getattr(new_model, par).frozen = False
 
-    logging.info(new_model.as_parfile())
+    try:
+        # This fails on windows
+        logging.info(new_model.as_parfile())
+    except Exception as e:
+        print(e)
+        pass
     return new_model
 
 
diff --git a/ell1fit/ell1fit.py b/ell1fit/ell1fit.py
index 9c21481..39b5138 100755
--- a/ell1fit/ell1fit.py
+++ b/ell1fit/ell1fit.py
@@ -8,6 +8,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 from hendrics.io import load_events
+from stingray.pulse.pulsar import z_n_events, z_n_binned_events
 from pint.models import get_model
 
 import matplotlib as mpl
@@ -17,7 +18,7 @@
 import emcee
 import corner
 from numba import njit, vectorize, int64, float32, float64, prange
-from astropy.table import Table, vstack
+from astropy.table import Table, vstack, TableMergeError
 from astropy.time import Time
 from scipy.ndimage import gaussian_filter
 from scipy.optimize import minimize
@@ -41,7 +42,7 @@
     "grid.color": "grey",
     "grid.linewidth": 0.3,
     "grid.linestyle": ":",
-    "axes.grid.axis": "y",
+    "axes.grid.axis": "both",
     "axes.grid.which": "both",
     "axes.axisbelow": False,
     "axes.labelsize": 8,
@@ -188,16 +189,11 @@ def phases_from_zero_to_one(phase):
     """Normalize pulse phases from 0 to 1
     Examples
     --------
-    >>> phases_from_zero_to_one(0.1)
-    0.1
-    >>> phases_from_zero_to_one(-0.9)
-    0.1
-    >>> phases_from_zero_to_one(0.9)
-    0.9
-    >>> phases_from_zero_to_one(3.1)
-    0.1
+    >>> assert np.isclose(phases_from_zero_to_one(0.1), 0.1)
+    >>> assert np.isclose(phases_from_zero_to_one(-0.9), 0.1)
+    >>> assert np.isclose(phases_from_zero_to_one(0.9), 0.9)
+    >>> assert np.isclose(phases_from_zero_to_one(3.1), 0.1)
     >>> assert np.allclose(phases_from_zero_to_one([0.1, 3.1, -0.9]), 0.1)
-    True
     """
 
     return phase - np.floor(phase)
@@ -208,14 +204,10 @@ def phases_around_zero(phase):
     """Normalize pulse phases from -0.5 to 0.5
     Examples
     --------
-    >>> phases_around_zero(0.6)
-    -0.4
-    >>> phases_around_zero(-0.9)
-    0.1
-    >>> phases_around_zero(3.9)
-    -0.1
-    >>> assert np.allclose(phases_from_zero_to_one([0.6, -0.4]), -0.4)
-    True
+    >>> assert np.isclose(phases_around_zero(0.6), -0.4)
+    >>> assert np.isclose(phases_around_zero(-0.9), 0.1)
+    >>> assert np.isclose(phases_around_zero(3.9), -0.1)
+    >>> assert np.allclose(phases_around_zero([0.6, -0.4]), -0.4)
     """
     ph = phase - np.floor(phase)
     while ph >= 0.5:
@@ -288,7 +280,9 @@ def template_func(x):
 
         phases_fine = np.arange(0.5 * dph_fine, 3, dph_fine)
 
-        templ_func_fine = interp1d(phases_fine, template_fine, kind="cubic", assume_sorted=True)
+        templ_func_fine = interp1d(
+            phases_fine, template_fine, kind="cubic", assume_sorted=True
+        )
 
         additional_phase = (
             np.argmax(template_fine[: final_nbin * oversample_factor])
@@ -311,10 +305,16 @@ def template_func(x):
     template = template[:final_nbin].real
 
     fig = plt.figure(figsize=(3.5, 2.65))
-    plt.plot(np.arange(0.5 / nbin, 1, 1 / nbin), profile, drawstyle="steps-mid", label="data")
+    plt.plot(
+        np.arange(0.5 / nbin, 1, 1 / nbin), profile, drawstyle="steps-mid", label="data"
+    )
     plt.plot(phas[:final_nbin], template, label="template values", ls="--", lw=2)
     plt.plot(
-        phas[:final_nbin], template_func(phas[:final_nbin]), label="template func", ls=":", lw=2
+        phas[:final_nbin],
+        template_func(phas[:final_nbin]),
+        label="template func",
+        ls=":",
+        lw=2,
     )
     plt.plot(
         phas[:final_nbin],
@@ -329,7 +329,7 @@ def template_func(x):
     return template * final_nbin / nbin, additional_phase
 
 
-def likelihood(phases, template_func, weights=None):
+def pletsch_clarke_likelihood(phases, template_func, weights=None):
     probs = template_func(phases)
     if weights is None:
         return np.log(probs).sum()
@@ -337,6 +337,11 @@ def likelihood(phases, template_func, weights=None):
         return np.log(weights * probs + 1.0 - weights).sum()
 
 
+def rayleigh_as_likelihood(phases, *args, **kwargs):
+    prob = z_n_events(phases, 1)
+    return prob
+
+
 def get_template_func(template):
     """Get a cubic interpolation function of a pulse template.
     Parameters
@@ -405,7 +410,9 @@ def simple_ell1_deorbit_numba(times, PB, A1, TASC, EPS1, EPS2, tolerance=1e-8):
             old_out = out_times[i]
             phase = omega * out_times[i]
             twophase = 2 * phase
-            out_times[i] = t - A1 * (np.sin(phase) + k1 * np.sin(twophase) + k2 * np.cos(twophase))
+            out_times[i] = t - A1 * (
+                np.sin(phase) + k1 * np.sin(twophase) + k2 * np.cos(twophase)
+            )
         out_times[i] += TASC
 
         # out_times[i] = times[i] - A1 * np.sin(omega * (out_times[i] - TASC))
@@ -457,7 +464,12 @@ def calculate_result_array_from_samples(sampler, labels):
 
 
 def plot_mcmc_results(
-    sampler=None, backend=None, flat_samples=None, labels=None, fname="results.jpg", **plot_kwargs
+    sampler=None,
+    backend=None,
+    flat_samples=None,
+    labels=None,
+    fname="results.jpg",
+    **plot_kwargs,
 ):
     assert np.any([a is not None for a in [sampler, backend, flat_samples]]), (
         "At least one between backend, sampler, or flat_samples, should be specified, in",
@@ -472,7 +484,9 @@ def plot_mcmc_results(
 
         flat_samples, _ = get_flat_samples(sampler)
 
-    fig = corner.corner(flat_samples, labels=labels, quantiles=[0.16, 0.5, 0.84], **plot_kwargs)
+    fig = corner.corner(
+        flat_samples, labels=labels, quantiles=[0.16, 0.5, 0.84], **plot_kwargs
+    )
     fig.savefig(fname, dpi=300)
 
 
@@ -485,7 +499,6 @@ def safe_run_sampler(
     corner_labels=None,
     n_autocorr=200,
 ):
-
     # https://emcee.readthedocs.io/en/stable/tutorials/monitor/?highlight=run_mcmc#saving-monitoring-progress
     # We'll track how the average autocorrelation time estimate changes
     starting_pars = np.asarray(starting_pars)
@@ -561,7 +574,10 @@ def safe_run_sampler(
 
     result_dict, flat_samples = calculate_result_array_from_samples(sampler, labels)
     plot_mcmc_results(
-        flat_samples=flat_samples, labels=labels, fname=outroot + "_corner.jpg", backend=backend
+        flat_samples=flat_samples,
+        labels=labels,
+        fname=outroot + "_corner.jpg",
+        backend=backend,
     )
 
     return result_dict
@@ -716,17 +732,21 @@ def fast_phase(times, frequency_derivatives):
     if len(frequency_derivatives) == 1:
         return _fast_phase(times, frequency_derivatives[0])
     elif len(frequency_derivatives) == 2:
-        return _fast_phase_fdot(times, frequency_derivatives[0], frequency_derivatives[1])
+        return _fast_phase_fdot(
+            times, frequency_derivatives[0], frequency_derivatives[1]
+        )
     elif len(frequency_derivatives) == 3:
         return _fast_phase_fddot(
-            times, frequency_derivatives[0], frequency_derivatives[1], frequency_derivatives[2]
+            times,
+            frequency_derivatives[0],
+            frequency_derivatives[1],
+            frequency_derivatives[2],
         )
 
     return _fast_phase_generic(times, np.array(frequency_derivatives))
 
 
 def _calculate_phases(times_from_pepoch, pars_dict, tolerance=1e-8):
-
     n_files = len(times_from_pepoch)
     list_phases_from_zero_to_one = []
     pb = pars_dict["PB"]
@@ -802,7 +822,10 @@ def return_unc(param):
         "EPS1": [model.EPS1.value.astype(float), return_unc(model.EPS1)],
         "EPS2": [model.EPS2.value.astype(float), return_unc(model.EPS2)],
         "PBDOT": [model.PBDOT.value.astype(float), return_unc(model.PBDOT)],
-        "PEPOCH": [model.PEPOCH.value.astype(float), return_unc(model.PEPOCH)],  # I added Pepoch
+        "PEPOCH": [
+            model.PEPOCH.value.astype(float),
+            return_unc(model.PEPOCH),
+        ],  # I added Pepoch
     }
 
     count = 0
@@ -819,16 +842,13 @@ def _load_and_format_events(
     event_file, energy_range, pepoch, plotlc=True, plotfile="lightcurve.jpg"
 ):
     events = load_events(event_file)
+    events.apply_gtis(inplace=True)
+
     if plotlc:
         lc = events.to_lc(100)
 
         fig = plt.figure("LC", figsize=(3.5, 2.65))
-        plt.plot(_sec_to_mjd(lc.time, events.mjdref), lc.counts / lc.dt)
-        GTI = _sec_to_mjd(events.gti, events.mjdref)
-        for g0, g1 in zip(GTI[:, 1], GTI[:, 0]):
-            plt.axvspan(g0, g1, color="r", alpha=0.5)
-        plt.xlabel("MJD")
-        plt.ylabel("Count rate")
+        lc.plot(ax=plt.gca())
         plt.savefig(plotfile)
         plt.close(fig)
 
@@ -854,6 +874,7 @@ def optimize_solution(
     nharm=1,
     outroot="out",
     tolerance=1e-8,
+    likelihood_func=pletsch_clarke_likelihood,
 ):
     def logprior(pars):
         if np.any(np.isnan(pars)):
@@ -883,7 +904,7 @@ def func_to_maximize(pars):
 
         ll = 0
         for i in range(len(phases)):
-            ll += likelihood(phases[i], template_func[i])
+            ll += likelihood_func(phases[i], template_func[i])
 
         return ll + lp
 
@@ -966,7 +987,6 @@ def func(x):
 def assign_logpriors(
     parnames, parvalunc
 ):  # parvalunc is a dictionary with mean values ([0]) and uncertainties ([1])of the parameters.
-
     logps = []
     logging.info("Setting up priors")
     for par in parnames:
@@ -990,7 +1010,9 @@ def assign_logpriors(
             logps.append(_flat_logprior(-np.inf, np.inf))
         else:
             log_line += f"normal with mean {parvalunc[par][0]} and std {abs(parvalunc[par][1]):.2e}"
-            logps.append(norm(loc=parvalunc[par][0], scale=abs(parvalunc[par][1])).logpdf)
+            logps.append(
+                norm(loc=parvalunc[par][0], scale=abs(parvalunc[par][1])).logpdf
+            )
         logging.info(log_line)
 
     return logps
@@ -1001,7 +1023,6 @@ def order_of_magnitude(value):
 
 
 def get_factors(parnames, model, observation_length):
-
     n_files = len(observation_length)
     zoom = []
     P = model[0].PB.value * 86400
@@ -1015,7 +1036,9 @@ def get_factors(parnames, model, observation_length):
         if matchobj:
             order = int(matchobj.group(1))
             file_n = int(matchobj.group(2))
-            zoom.append(order_of_magnitude(1 / observation_length[file_n] ** (order + 1)))
+            zoom.append(
+                order_of_magnitude(1 / observation_length[file_n] ** (order + 1))
+            )
         elif par == "A1":
             zoom.append(min(1, order_of_magnitude(1 / np.pi / 2 / F)))
         elif par == "PB":
@@ -1086,7 +1109,9 @@ def split_output_results(result_table, n_files, fit_parameters):
         for i in list(range(n_files))[::-1]:
             par_to_test = f"{par}_{i}"
 
-            cols = look_for_string_in_list_of_strings(common_table.colnames, par_to_test)
+            cols = look_for_string_in_list_of_strings(
+                common_table.colnames, par_to_test
+            )
             for colname in cols:
                 clean_colname = colname.replace(f"{par}_{i}", f"{par}")
 
@@ -1100,6 +1125,35 @@ def split_output_results(result_table, n_files, fit_parameters):
     return output_tables
 
 
+def safe_save(results, output_file, **write_kwargs):
+    """
+    Examples
+    --------
+    >>> results = Table({"a": [2]})
+    >>> results_2 = Table({"a": ["3"]})
+    >>> output_file = "blabla.csv"
+    >>> safe_save(results, output_file)
+    >>> safe_save(results, output_file)
+    >>> out = Table.read(output_file)
+    >>> len(out)
+    2
+    >>> os.unlink(output_file)
+    >>> os.unlink("old_" + output_file)
+    """
+    if os.path.exists(output_file):
+        old = Table.read(output_file)
+        old.write("old_" + output_file, overwrite=True)
+        try:
+            results = vstack([old, results])
+        except TableMergeError:
+            warnings.warn(
+                "Merging old and new results failed. Old results were saved in a separate file."
+            )
+
+    results.write(output_file, overwrite=True, **write_kwargs)
+    return
+
+
 def ell1fit(
     files,
     parfiles,
@@ -1110,6 +1164,7 @@ def ell1fit(
     fit_parameters=["F0"],
     minimize_first=False,
     general_outroot=None,
+    likelihood_func=pletsch_clarke_likelihood,
 ):
     n_files = len(files)
     assert len(parfiles) == len(
@@ -1123,13 +1178,19 @@ def ell1fit(
         pepoch.append(model[i].PEPOCH.value)
 
         if hasattr(model[i], "T0") or model[i].BINARY.value != "ELL1":
-            raise ValueError("This script wants an ELL1 model, with TASC, not T0, defined")
+            raise ValueError(
+                "This script wants an ELL1 model, with TASC, not T0, defined"
+            )
 
         model[i].change_binary_epoch(pepoch[i])
 
     nbin = max(16, nharm * 8)
 
     energy_str = _format_energy_string(energy_range)
+    likelihood_str = ""
+    if likelihood_func == rayleigh_as_likelihood:
+        likelihood_str = "_rayleigh"
+
     nharm_str = ""
     if nharm > 1:
         nharm_str = f"_N{nharm}"
@@ -1172,7 +1233,7 @@ def ell1fit(
     list_parameter_names = sorted(fit_parameters)
 
     for f in parameters:
-        if f.startswith("Phase"):
+        if f.startswith("Phase") and likelihood_func == pletsch_clarke_likelihood:
             parameter_names.append(f)
             continue
         for g in list_parameter_names:
@@ -1188,7 +1249,14 @@ def get_outroot(file_n=None):
         else:
             initial_outroot = "out"
 
-        outroot = initial_outroot + "_" + "_".join(list_parameter_names) + energy_str + nharm_str
+        outroot = (
+            initial_outroot
+            + "_"
+            + "_".join(list_parameter_names)
+            + energy_str
+            + nharm_str
+            + likelihood_str
+        )
         return outroot
 
     times_from_pepoch = [[] for _ in range(n_files)]
@@ -1197,7 +1265,10 @@ def get_outroot(file_n=None):
     for i in range(n_files):
         fname = files[i]
         times_from_pepoch[i], gtis = _load_and_format_events(
-            fname, energy_range, pepoch[i], plotfile=get_outroot(i) + f"_lightcurve_{i}.jpg"
+            fname,
+            energy_range,
+            pepoch[i],
+            plotfile=get_outroot(i) + f"_lightcurve_{i}.jpg",
         )
         expo[i] += np.sum(np.diff(gtis, axis=1))
 
@@ -1206,14 +1277,19 @@ def get_outroot(file_n=None):
     logprior_funcs = assign_logpriors(parameter_names, parameters_with_unc)
     factors = get_factors(parameter_names, model, observation_length)
 
-    profile = folded_profile(times_from_pepoch, parameters, nbin=nbin, tolerance=tolerance)
+    profile = folded_profile(
+        times_from_pepoch, parameters, nbin=nbin, tolerance=tolerance
+    )
 
     template_func = []
     pulsed_frac = []
 
     for i in range(n_files):
         template, additional_phase = create_template_from_profile_harm(
-            profile[i], nharm=nharm, final_nbin=200, imagefile=get_outroot(i) + "_template.jpg"
+            profile[i],
+            nharm=nharm,
+            final_nbin=200,
+            imagefile=get_outroot(i) + "_template.jpg",
         )
 
         template_func.append(get_template_func(template))
@@ -1247,6 +1323,7 @@ def get_outroot(file_n=None):
         nharm=nharm,
         outroot=[get_outroot(i) for i in range(n_files)] + [get_outroot(None)],
         tolerance=tolerance,
+        likelihood_func=likelihood_func,
     )
 
     for i in range(n_files):
@@ -1272,6 +1349,7 @@ def get_outroot(file_n=None):
 
     for i in range(n_files):
         results[f"pf_{i}"] = pulsed_frac[i]
+        results[f"Z11_{i}"] = z_n_binned_events(profile[i], nharm)
 
     for i in range(n_files):
         results[f"ctrate_{i}"] = times_from_pepoch[i].size / expo[i]
@@ -1286,12 +1364,7 @@ def get_outroot(file_n=None):
 
     output_file = get_outroot(None) + "_results.ecsv"
 
-    if os.path.exists(output_file):
-        old = Table.read(output_file)
-        old.write("old_" + output_file, overwrite=True)
-        results = vstack([old, results])
-
-    results.write(output_file, overwrite=True)
+    safe_save(results, output_file)
 
     list_result = split_output_results(results, n_files, list_parameter_names)
 
@@ -1308,7 +1381,9 @@ def main(args=None):
     """Main function called by the `ell1fit` script"""
     import argparse
 
-    description = "Fit an ELL1 model and frequency derivatives to an X-ray " "pulsar observation."
+    description = (
+        "Fit an ELL1 model and frequency derivatives to an X-ray " "pulsar observation."
+    )
     parser = argparse.ArgumentParser(description=description)
     parser.add_argument("files", help="List of files", nargs="+")
     parser.add_argument(
@@ -1323,7 +1398,9 @@ def main(args=None):
             "All other models will be ignored."
         ),
     )
-    parser.add_argument("-o", "--outroot", type=str, default=None, help="Root of output file names")
+    parser.add_argument(
+        "-o", "--outroot", type=str, default=None, help="Root of output file names"
+    )
     parser.add_argument(
         "-N",
         "--nharm",
@@ -1358,12 +1435,22 @@ def main(args=None):
         help="Comma-separated list of parameters to fit",
         default="F0,F1",
     )
+    parser.add_argument(
+        "--likelihood",
+        type=str,
+        help="Can be PC (Pletsch & Clarke, default) or Rayleigh",
+        default="PC",
+    )
     parser.add_argument("--minimize-first", action="store_true", default=False)
 
     args = parser.parse_args(args)
     files = args.files
     parfiles = args.parfile
 
+    like = pletsch_clarke_likelihood
+    if args.likelihood.lower() == "rayleigh":
+        like = rayleigh_as_likelihood
+
     ell1fit(
         files,
         parfiles,
@@ -1374,4 +1461,5 @@ def main(args=None):
         fit_parameters=args.parameters.split(","),
         minimize_first=args.minimize_first,
         general_outroot=args.outroot,
+        likelihood_func=like,
     )
diff --git a/ell1fit/tests/test_ell1fit.py b/ell1fit/tests/test_ell1fit.py
index 04840c3..6fb206f 100644
--- a/ell1fit/tests/test_ell1fit.py
+++ b/ell1fit/tests/test_ell1fit.py
@@ -1,6 +1,9 @@
+import os
+
 import pytest
 import numpy as np
-from ..ell1fit import add_circular_orbit_numba, add_ell1_orbit_numba
+from astropy.table import Table
+from ..ell1fit import add_circular_orbit_numba, add_ell1_orbit_numba, safe_save
 from ..ell1fit import simple_circular_deorbit_numba, simple_ell1_deorbit_numba
 
 
@@ -28,3 +31,15 @@ def test_ell1_orbit(PB, A1, E1, E2):
     orbited = add_ell1_orbit_numba(times, PB, A1, TASC, E1, E2)
     deorbited = simple_ell1_deorbit_numba(orbited, PB, A1, TASC, E1, E2, tolerance=1e-8)
     assert np.all(np.abs(deorbited - times) < 1e-8)
+
+
+def test_safe_save():
+    results = Table({"a": [2]})
+    results_2 = Table({"a": ["3"]})
+    output_file = "blabla.csv"
+    safe_save(results, output_file)
+    safe_save(results_2, output_file)
+    with pytest.warns(UserWarning, match="Merging old and new"):
+        safe_save(results_2, output_file)
+    os.unlink(output_file)
+    os.unlink("old_" + output_file)
diff --git a/ell1fit/tests/test_execution.py b/ell1fit/tests/test_execution.py
index 271fb5c..9eb1c18 100644
--- a/ell1fit/tests/test_execution.py
+++ b/ell1fit/tests/test_execution.py
@@ -1,21 +1,33 @@
 import os
+import sys
 import glob
+import pytest
 from ell1fit.ell1fit import main as main_ell1fit
 from ell1fit.create_parfile import main as main_ell1par
 
 
 curdir = os.path.abspath(os.path.dirname(__file__))
 datadir = os.path.join(curdir, "data")
+if sys.platform.startswith("win"):
+    pytest.skip(
+        "skipping tests that are known to fail on windows", allow_module_level=True
+    )
 
 
-class TestExecution():
+class TestExecution:
     @classmethod
     def setup_class(cls):
         cls.event_files = sorted(glob.glob(os.path.join(datadir, "events[01].nc")))
         cls.param_files = sorted(glob.glob(os.path.join(datadir, "events[01].par")))
 
-    def test_ell1fit_and_ell1par(self):
-        cmdlines = self.event_files + ["-p"] + self.param_files + ["-P", "F0,PB,A1,TASC"]
+    @pytest.mark.parametrize("likelihood", ["PC", "Rayleigh"])
+    def test_ell1fit_and_ell1par(self, likelihood):
+        cmdlines = (
+            self.event_files
+            + ["-p"]
+            + self.param_files
+            + ["-P", "F0,PB,A1,TASC", "--likelihood", likelihood]
+        )
 
         cmdline1 = cmdlines + ["--nsteps", "100"]
         cmdline2 = cmdlines + ["--nsteps", "200"]
@@ -23,20 +35,27 @@ def test_ell1fit_and_ell1par(self):
         # Get to 100, then continue up to 200
         main_ell1fit(cmdline1)
         main_ell1fit(cmdline2)
+        label = "_A1_F0_PB_TASC"
+        if likelihood == "Rayleigh":
+            label += "_rayleigh"
 
-        outputs = sorted(glob.glob(os.path.join(datadir, "events[01]_A1_F0_PB_TASC_results.ecsv")))
+        outputs = sorted(
+            glob.glob(os.path.join(datadir, f"events[01]{label}_results.ecsv"))
+        )
         for out in outputs:
             assert os.path.exists(out)
 
         main_ell1par(f"{outputs[0]} -p {self.param_files[0]}".split())
         main_ell1par(f"{outputs[1]} -p {self.param_files[1]}".split())
 
-        out_param = sorted(glob.glob(os.path.join(datadir, "events[01]_A1_F0_PB_TASC_results.par")))
+        out_param = sorted(
+            glob.glob(os.path.join(datadir, "events[01]_A1_F0_PB_TASC_results.par"))
+        )
         for out in out_param:
             assert os.path.exists(out)
 
     @classmethod
     def teardown_class(cls):
-        outs = glob.glob(os.path.join(datadir, '*A1_*TASC*'))
+        outs = glob.glob(os.path.join(datadir, "*A1_*TASC*"))
         for out in outs:
             os.remove(out)
diff --git a/setup.cfg b/setup.cfg
index 28525e2..ce784d9 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -17,7 +17,7 @@ packages = find:
 python_requires = >=3.7
 setup_requires = setuptools_scm
 install_requires =
-    astropy
+    astropy<6
     emcee
     matplotlib
     pint-pulsar
diff --git a/tox.ini b/tox.ini
index 041c7db..14ac5ac 100644
--- a/tox.ini
+++ b/tox.ini
@@ -1,8 +1,6 @@
 [tox]
 envlist =
-    py{36,37,38}-test{,-alldeps,-devdeps}{,-cov}
-    py{36,37,38}-test-numpy{116,117,118}
-    py{36,37,38}-test-astropy{30,40,lts}
+    py{39,310,311,312}-test{,-alldeps,-devdeps}{,-cov}
     build_docs
     linkcheck
     codestyle
@@ -67,6 +65,7 @@ extras =
 
 commands =
     pip freeze
+    pip install git+https://github.com/stingraysoftware/stingray.git
     !cov: pytest --pyargs ell1fit {toxinidir}/docs {posargs}
     cov: pytest --pyargs ell1fit {toxinidir}/docs --cov ell1fit --cov-config={toxinidir}/setup.cfg {posargs}
     cov: coverage xml -o {toxinidir}/coverage.xml