diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 7b9c62f..785760c 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,6 +1,6 @@
 repos:
 - repo: https://github.com/pre-commit/pre-commit-hooks
-  rev: v4.3.0
+  rev: v4.4.0
   hooks:
   - id: trailing-whitespace
   - id: end-of-file-fixer
@@ -10,7 +10,7 @@ repos:
   - id: check-ast
   - id: check-merge-conflict
 - repo: https://github.com/psf/black
-  rev: 22.10.0
+  rev: 23.1.0
   hooks:
   - id: black
     exclude: venv*/|dustmaps/|grids/
diff --git a/basta/_version.py b/basta/_version.py
index c68196d..a955fda 100644
--- a/basta/_version.py
+++ b/basta/_version.py
@@ -1 +1 @@
-__version__ = "1.2.0"
+__version__ = "1.2.1"
diff --git a/basta/bastamain.py b/basta/bastamain.py
index a0bb657..c97b9fe 100644
--- a/basta/bastamain.py
+++ b/basta/bastamain.py
@@ -11,7 +11,7 @@
 import numpy as np
 from tqdm import tqdm
 
-from basta import freq_fit, stats, process_output, priors, distances, plot_seismic
+from basta import freq_fit, stats, process_output, priors, distances, plot_driver
 from basta import utils_seismic as su
 from basta import utils_general as util
 from basta._version import __version__
@@ -264,6 +264,15 @@ def BASTA(
             obsintervals,
         ) = su.prepare_obs(inputparams, verbose=verbose, debug=debug)
 
+        # Apply prior on dnufit to mimick the range defined by dnufrac
+        if fitfreqs["dnuprior"] and ("dnufit" not in limits):
+            dnufit_frac = fitfreqs["dnufrac"] * fitfreqs["dnufit"]
+            dnuerr = max(3 * fitfreqs["dnufit_err"], dnufit_frac)
+            limits["dnufit"] = [
+                fitfreqs["dnufit"] - dnuerr,
+                fitfreqs["dnufit"] + dnuerr,
+            ]
+
     # Check if grid is interpolated
     try:
         Grid["header/interpolation_time"][()]
@@ -348,6 +357,7 @@ def BASTA(
 
         # Translate True/False to Yes/No
         strmap = ("No", "Yes")
+        print("  - Automatic prior on dnu: {0}".format(strmap[fitfreqs["dnuprior"]]))
         print(
             "  - Constraining lowest l = 0 (n = {0}) with f = {1:.3f} +/-".format(
                 obskey[1, 0], obs[0, 0]
@@ -383,7 +393,14 @@ def BASTA(
                 strmap[fitfreqs["threepoint"]]
             )
         )
-        print("  - Value of dnu: {0:.3f} microHz".format(fitfreqs["dnufit"]))
+        if fitfreqs["dnufit_err"]:
+            print(
+                "  - Value of dnu: {0:.3f} +/- {1:.3f} microHz".format(
+                    fitfreqs["dnufit"], fitfreqs["dnufit_err"]
+                )
+            )
+        else:
+            print("  - Value of dnu: {0:.3f} microHz".format(fitfreqs["dnufit"]))
         print("  - Value of numax: {0:.3f} microHz".format(fitfreqs["numax"]))
 
         weightcomment = ""
@@ -503,7 +520,7 @@ def BASTA(
         trackcounter += len(group.items())
 
     # Prepare the main loop
-    shapewarn = False
+    shapewarn = 0
     warn = True
     selectedmodels = {}
     noofind = 0
@@ -527,7 +544,7 @@ def BASTA(
 
             # For grid with interpolated tracks, skip tracks flagged as empty
             if grid_is_intpol:
-                if libitem["FeHini_weight"][()] < 0:
+                if libitem["IntStatus"][()] < 0:
                     continue
 
             # Check for diffusion
@@ -576,7 +593,6 @@ def BASTA(
 
             # Check which models have phases as specified
             if "phase" in inputparams:
-
                 # Mapping of verbose input phases to internal numbers
                 pmap = {
                     "pre-ms": 1,
@@ -754,7 +770,7 @@ def BASTA(
         )
 
     # Raise possible warnings
-    if shapewarn:
+    if shapewarn == 1:
         print(
             "Warning: Found models with fewer frequencies than observed!",
             "These were set to zero likelihood!",
@@ -764,6 +780,11 @@ def BASTA(
                 "This is probably due to the interpolation scheme. Lookup",
                 "`interpolate_frequencies` for more details.",
             )
+    if shapewarn == 2:
+        print(
+            "Warning: Models without frequencies overlapping with observed",
+            "ignored due to interpolation of ratios being impossible.",
+        )
     if noofposind == 0:
         fio.no_models(starid, inputparams, "No models found")
         return
@@ -777,9 +798,9 @@ def BASTA(
 
     # Find and print highest likelihood model info
     maxPDF_path, maxPDF_ind = stats.get_highest_likelihood(
-        Grid, selectedmodels, outparams
+        Grid, selectedmodels, inputparams
     )
-    stats.get_lowest_chi2(Grid, selectedmodels, outparams)
+    stats.get_lowest_chi2(Grid, selectedmodels, inputparams)
 
     # Generate posteriors of ascii- and plotparams
     # --> Print posteriors to console and log
@@ -798,219 +819,24 @@ def BASTA(
         experimental=experimental,
         validationmode=validationmode,
     )
-
     # Make frequency-related plots
     freqplots = inputparams.get("freqplots")
     if fitfreqs["active"] and len(freqplots):
-        # Check which plots to create
-        allfplots = freqplots[0] == True
-        if any(x == "allechelle" for x in freqplots):
-            freqplots += ["dupechelle", "echelle", "pairechelle"]
-        if any(x in freqtypes.rtypes for x in freqplots):
-            freqplots += ["ratios"]
-
-        # Naming of plots preparation
-        plotfmt = inputparams["plotfmt"]
-        plotfname = outfilename + "_{0}." + plotfmt
-
-        rawmaxmod = Grid[maxPDF_path + "/osc"][maxPDF_ind]
-        rawmaxmodkey = Grid[maxPDF_path + "/osckey"][maxPDF_ind]
-        maxmod = su.transform_obj_array(rawmaxmod)
-        maxmodkey = su.transform_obj_array(rawmaxmodkey)
-        maxmod = maxmod[:, maxmodkey[0, :] < 2.5]
-        maxmodkey = maxmodkey[:, maxmodkey[0, :] < 2.5]
-        maxjoins = freq_fit.calc_join(
-            mod=maxmod,
-            modkey=maxmodkey,
-            obs=obs,
+        plot_driver.plot_all_seismic(
+            freqplots,
+            Grid=Grid,
+            fitfreqs=fitfreqs,
+            obsfreqmeta=obsfreqmeta,
+            obsfreqdata=obsfreqdata,
             obskey=obskey,
+            obs=obs,
             obsintervals=obsintervals,
+            selectedmodels=selectedmodels,
+            path=maxPDF_path,
+            ind=maxPDF_ind,
+            plotfname=outfilename + "_{0}." + inputparams["plotfmt"],
+            debug=debug,
         )
-        maxjoinkeys, maxjoin = maxjoins
-        maxmoddnu = Grid[maxPDF_path + "/dnufit"][maxPDF_ind]
-
-        if allfplots or "echelle" in freqplots:
-            plot_seismic.echelle(
-                selectedmodels,
-                Grid,
-                obs,
-                obskey,
-                mod=maxmod,
-                modkey=maxmodkey,
-                dnu=fitfreqs["dnufit"],
-                join=maxjoin,
-                joinkeys=maxjoinkeys,
-                pair=False,
-                duplicate=False,
-                output=plotfname.format("echelle_uncorrected"),
-            )
-        if allfplots or "pairechelle" in freqplots:
-            plot_seismic.echelle(
-                selectedmodels,
-                Grid,
-                obs,
-                obskey,
-                mod=maxmod,
-                modkey=maxmodkey,
-                dnu=fitfreqs["dnufit"],
-                join=maxjoin,
-                joinkeys=maxjoinkeys,
-                pair=True,
-                duplicate=False,
-                output=plotfname.format("pairechelle_uncorrected"),
-            )
-        if allfplots or "dupechelle" in freqplots:
-            plot_seismic.echelle(
-                selectedmodels,
-                Grid,
-                obs,
-                obskey,
-                mod=maxmod,
-                modkey=maxmodkey,
-                dnu=fitfreqs["dnufit"],
-                join=maxjoin,
-                joinkeys=maxjoinkeys,
-                duplicate=True,
-                pair=True,
-                output=plotfname.format("dupechelle_uncorrected"),
-            )
-
-        if fitfreqs["fcor"] == "None":
-            corjoin = maxjoin
-            coeffs = [1]
-        elif fitfreqs["fcor"] == "HK08":
-            corjoin, coeffs = freq_fit.HK08(
-                joinkeys=maxjoinkeys,
-                join=maxjoin,
-                nuref=fitfreqs["numax"],
-                bcor=fitfreqs["bexp"],
-            )
-        elif fitfreqs["fcor"] == "BG14":
-            corjoin, coeffs = freq_fit.BG14(
-                joinkeys=maxjoinkeys, join=maxjoin, scalnu=fitfreqs["numax"]
-            )
-        elif fitfreqs["fcor"] == "cubicBG14":
-            corjoin, coeffs = freq_fit.cubicBG14(
-                joinkeys=maxjoinkeys, join=maxjoin, scalnu=fitfreqs["numax"]
-            )
-
-        if len(coeffs) > 1:
-            print("The surface correction coefficients are", *coeffs)
-        else:
-            print("The surface correction coefficient is", *coeffs)
-
-        if allfplots or "echelle" in freqplots:
-            plot_seismic.echelle(
-                selectedmodels,
-                Grid,
-                obs,
-                obskey,
-                mod=maxmod,
-                modkey=maxmodkey,
-                dnu=fitfreqs["dnufit"],
-                join=corjoin,
-                joinkeys=maxjoinkeys,
-                freqcor=fitfreqs["fcor"],
-                coeffs=coeffs,
-                scalnu=fitfreqs["numax"],
-                pair=False,
-                duplicate=False,
-                output=plotfname.format("echelle"),
-            )
-        if allfplots or "pairechelle" in freqplots:
-            plot_seismic.echelle(
-                selectedmodels,
-                Grid,
-                obs,
-                obskey,
-                mod=maxmod,
-                modkey=maxmodkey,
-                dnu=fitfreqs["dnufit"],
-                join=corjoin,
-                joinkeys=maxjoinkeys,
-                freqcor=fitfreqs["fcor"],
-                coeffs=coeffs,
-                scalnu=fitfreqs["numax"],
-                pair=True,
-                duplicate=False,
-                output=plotfname.format("pairechelle"),
-            )
-        if allfplots or "dupechelle" in freqplots:
-            plot_seismic.echelle(
-                selectedmodels,
-                Grid,
-                obs,
-                obskey,
-                mod=maxmod,
-                modkey=maxmodkey,
-                dnu=fitfreqs["dnufit"],
-                join=corjoin,
-                joinkeys=maxjoinkeys,
-                freqcor=fitfreqs["fcor"],
-                coeffs=coeffs,
-                scalnu=fitfreqs["numax"],
-                duplicate=True,
-                pair=True,
-                output=plotfname.format("dupechelle"),
-            )
-        if "freqcormap" in freqplots or debug:
-            plot_seismic.correlation_map(
-                "freqs",
-                obsfreqdata,
-                plotfname.format("freqs_cormap"),
-                obskey=obskey,
-            )
-        if obsfreqmeta["getratios"]:
-            for ratseq in obsfreqmeta["ratios"]["plot"]:
-                ratnamestr = "ratios_{0}".format(ratseq)
-                plot_seismic.ratioplot(
-                    obsfreqdata,
-                    maxjoinkeys,
-                    maxjoin,
-                    maxmodkey,
-                    maxmod,
-                    ratseq,
-                    output=plotfname.format(ratnamestr),
-                    threepoint=fitfreqs["threepoint"],
-                    interp_ratios=fitfreqs["interp_ratios"],
-                )
-                if fitfreqs["correlations"]:
-                    plot_seismic.correlation_map(
-                        ratseq,
-                        obsfreqdata,
-                        output=plotfname.format(ratnamestr + "_cormap"),
-                    )
-
-        if obsfreqmeta["getepsdiff"]:
-            for epsseq in obsfreqmeta["epsdiff"]["plot"]:
-                epsnamestr = "epsdiff_{0}".format(epsseq)
-                plot_seismic.epsilon_difference_diagram(
-                    mod=maxmod,
-                    modkey=maxmodkey,
-                    moddnu=maxmoddnu,
-                    sequence=epsseq,
-                    obsfreqdata=obsfreqdata,
-                    output=plotfname.format(epsnamestr),
-                )
-                if fitfreqs["correlations"]:
-                    plot_seismic.correlation_map(
-                        epsseq,
-                        obsfreqdata,
-                        output=plotfname.format(epsnamestr + "_cormap"),
-                    )
-        if obsfreqmeta["getepsdiff"] and debug:
-            if len(obsfreqmeta["epsdiff"]["plot"]) > 0:
-                plot_seismic.epsilon_difference_components_diagram(
-                    mod=maxmod,
-                    modkey=maxmodkey,
-                    moddnu=maxmoddnu,
-                    obs=obs,
-                    obskey=obskey,
-                    dnudata=obsfreqdata["freqs"]["dnudata"],
-                    obsfreqdata=obsfreqdata,
-                    obsfreqmeta=obsfreqmeta,
-                    output=plotfname.format("DEBUG_epsdiff_components"),
-                )
     else:
         print(
             "Did not get any frequency file input, skipping ratios and echelle plots."
diff --git a/basta/constants.py b/basta/constants.py
index e01bf7e..98a6e12 100644
--- a/basta/constants.py
+++ b/basta/constants.py
@@ -32,6 +32,27 @@ class freqtypes:
     surfeffcorrs = ["HK08", "BG14", "cubicBG14"]
 
 
+@dataclass
+class statdata:
+    """
+    Constant values for statistics, to ensure consistensy across code
+
+    Contains
+    --------
+    quantiles : list
+        Median, lower and upper percentiles of Bayesian posterior
+        distributions to draw
+    nsamples : int
+        Number of samples to draw when sampling
+    nsigma : float
+        Fractional standard deviation used for smoothing
+    """
+
+    quantiles = [0.5, 0.158655, 0.841345]
+    nsamples = 100000
+    nsigma = 0.25
+
+
 @dataclass
 class parameters:
     """
@@ -254,7 +275,6 @@ def exclude_params(excludeparams):
         return parnames
 
     def get_keys(inputparams):
-
         """
         Takes a list of input parameters (or a
         single parameter) as strings and returns
diff --git a/basta/distances.py b/basta/distances.py
index e95bde4..e4cb97a 100644
--- a/basta/distances.py
+++ b/basta/distances.py
@@ -60,8 +60,8 @@ def LOS_reddening(distanceparams):
     if "EBV" in distanceparams:
         return lambda x: np.asarray(
             np.random.normal(
-                distanceparams["EBV"][1],
-                distanceparams["EBV"][2] - distanceparams["EBV"][1],
+                distanceparams["EBV"][0],
+                distanceparams["EBV"][2] - distanceparams["EBV"][0],
                 size=[
                     len(i) if isinstance(i, collections.Iterable) else 1 for i in [x]
                 ][0],
@@ -275,7 +275,6 @@ def add_absolute_magnitudes(
 
     print("\nPreparing distance/parallax/magnitude input ...", flush=True)
 
-    qs = [0.158655, 0.5, 0.841345]
     fitparams = inputparams["fitparams"]
     distanceparams = inputparams["distanceparams"]
 
@@ -387,22 +386,28 @@ def add_absolute_magnitudes(
             plt.savefig(outfilename + "_DEBUG_absms_%s.png" % filt)
             plt.close()
 
-        absms_qs = stats.quantile_1D(absms, labsms, qs)
-        distanceparams["As"][filt] = list(stats.quantile_1D(As, labsms, qs))
+        absms_qs = stats.quantile_1D(absms, labsms, cnsts.statdata.quantiles)
+        distanceparams["As"][filt] = list(
+            stats.quantile_1D(As, labsms, cnsts.statdata.quantiles)
+        )
 
         # Extended dictionary for use in Kiel diagram, and clarifying it as a prior
         inputparams["magnitudes"][filt] = {
             "prior": M_prior,
-            "median": absms_qs[1],
-            "errp": np.abs(absms_qs[2] - absms_qs[1]),
-            "errm": np.abs(absms_qs[1] - absms_qs[0]),
+            "median": absms_qs[0],
+            "errp": np.abs(absms_qs[2] - absms_qs[0]),
+            "errm": np.abs(absms_qs[0] - absms_qs[1]),
         }
 
     # Get an estimate from all filters
     labsms_joined = np.exp(llabsms_joined - np.log(np.sum(np.exp(llabsms_joined))))
 
-    distanceparams["EBV"] = list(stats.quantile_1D(EBVs, labsms_joined, qs))
-    distanceparams["priordistance"] = list(stats.quantile_1D(dists, labsms_joined, qs))
+    distanceparams["EBV"] = list(
+        stats.quantile_1D(EBVs, labsms_joined, cnsts.statdata.quantiles)
+    )
+    distanceparams["priordistance"] = list(
+        stats.quantile_1D(dists, labsms_joined, cnsts.statdata.quantiles)
+    )
 
     # Constrain metallicity within the limits of the color transformations
     metal = "MeH" if "MeH" in inputparams["limits"] else "FeH"
diff --git a/basta/fileio.py b/basta/fileio.py
index bc0a94c..e5f60cf 100644
--- a/basta/fileio.py
+++ b/basta/fileio.py
@@ -705,7 +705,7 @@ def read_freq(filename, nottrustedfile=None, covarfre=False):
         f = np.concatenate([f, frecu[given_l][incrn]])
         e = np.concatenate([e, errors[given_l][incrn]])
     assert len(f) == len(n) == len(e) == len(l)
-    obskey = np.asarray([l, n], dtype=np.int)
+    obskey = np.asarray([l, n], dtype=int)
     obs = np.array([f, e])
 
     # Remove untrusted frequencies
@@ -722,7 +722,7 @@ def read_freq(filename, nottrustedfile=None, covarfre=False):
             print("Not-trusted frequencies found")
             if nottrustedobskey.shape == (2,):
                 nottrustedobskey = [nottrustedobskey]
-            for (l, n) in nottrustedobskey:
+            for l, n in nottrustedobskey:
                 nottrustedmask = (obskey[0] == l) & (obskey[1] == n)
                 print("Removed mode at", obs[:, nottrustedmask][0][0], "µHz")
                 obskey = obskey[:, ~nottrustedmask]
@@ -1117,9 +1117,8 @@ def read_glitch(filename):
     tmpFit[:, 1] = glitchfit[:, 4]
     tmpFit[:, 2] = glitchfit[:, 5]
     cov = MinCovDet().fit(tmpFit).covariance_
-    covinv = np.linalg.pinv(cov, rcond=1e-8)
 
-    return glitchparams, cov, covinv
+    return glitchparams, cov
 
 
 def read_precomputedratios(
@@ -1443,7 +1442,6 @@ def read_allseismic(
             if datos is not None:
                 obsfreqdata[ratiotype]["data"] = datos[0]
                 obsfreqdata[ratiotype]["cov"] = datos[1]
-                obsfreqdata[ratiotype]["covinv"] = datos[2]
             else:
                 if ratiotype in obsfreqmeta["ratios"]["fit"]:
                     # Fail
@@ -1463,7 +1461,6 @@ def read_allseismic(
         datos = read_glitch(fitfreqs["glhfile"])
         obsfreqdata["glitches"]["data"] = datos[0]
         obsfreqdata["glitches"]["cov"] = datos[1]
-        obsfreqdata["glitches"]["covinv"] = datos[2]
 
     # Get epsilon differences
     if obsfreqmeta["getepsdiff"]:
@@ -1482,7 +1479,6 @@ def read_allseismic(
                 )
                 obsfreqdata[epsdifffit]["data"] = datos[0]
                 obsfreqdata[epsdifffit]["cov"] = datos[1]
-                obsfreqdata[epsdifffit]["covinv"] = datos[2]
             elif epsdifffit in obsfreqmeta["epsdiff"]["plot"]:
                 datos = freq_fit.compute_epsilondiff(
                     obskey,
@@ -1495,14 +1491,18 @@ def read_allseismic(
                 )
                 obsfreqdata[epsdifffit]["data"] = datos[0]
                 obsfreqdata[epsdifffit]["cov"] = datos[1]
-                obsfreqdata[epsdifffit]["covinv"] = datos[2]
 
+    # As the inverse covariance matrix is actually what is used, compute it once
     # Diagonalise covariance matrices if correlations is set to False
-    if not fitfreqs["correlations"]:
-        for key in obsfreqdata.keys():
-            for mat in ["cov", "covinv"]:
-                full = obsfreqdata[key][mat]
-                obsfreqdata[key][mat] = np.identity(full.shape[0]) * full
+    for key in obsfreqdata.keys():
+        if not fitfreqs["correlations"]:
+            obsfreqdata[key]["cov"] = (
+                np.identity(obsfreqdata[key]["cov"].shape[0]) * obsfreqdata[key]["cov"]
+            )
+        if "covinv" not in obsfreqdata[key].keys():
+            obsfreqdata[key]["covinv"] = np.linalg.pinv(
+                obsfreqdata[key]["cov"], rcond=1e-8
+            )
 
     return obskey, obs, obsfreqdata, obsfreqmeta
 
diff --git a/basta/freq_fit.py b/basta/freq_fit.py
index a8e77b1..b3fe6e4 100644
--- a/basta/freq_fit.py
+++ b/basta/freq_fit.py
@@ -189,7 +189,7 @@ def calc_join(mod, modkey, obs, obskey, obsintervals=None, dnu=None):
 
                 joinkeys.append(
                     np.transpose(
-                        [l * np.ones(osum, dtype=np.int), nmod[mfilter], nobs[ofilter]]
+                        [l * np.ones(osum, dtype=int), nmod[mfilter], nobs[ofilter]]
                     )
                 )
                 join.append(
@@ -662,12 +662,14 @@ def compute_ratios(obskey, obs, ratiotype, nrealisations=10000, threepoint=False
         Ratio requested from `ratiotype`.
     ratio_cov : array
         Covariance matrix of the requested ratio.
-    ratio_covinv : array
-        The inverse of the covariance matrix of the requested ratio.
     """
     ratio = compute_ratioseqs(obskey, obs, ratiotype, threepoint=threepoint)
 
-    ratio_cov, ratio_covinv = su.compute_cov_from_mc(
+    # Check for valid return
+    if ratio is None:
+        return None
+
+    ratio_cov = su.compute_cov_from_mc(
         ratio.shape[1],
         obskey,
         obs,
@@ -675,7 +677,7 @@ def compute_ratios(obskey, obs, ratiotype, nrealisations=10000, threepoint=False
         args={"threepoint": threepoint},
         nrealisations=nrealisations,
     )
-    return ratio, ratio_cov, ratio_covinv
+    return ratio, ratio_cov
 
 
 def compute_ratioseqs(obskey, obs, sequence, threepoint=False):
@@ -720,15 +722,15 @@ def compute_ratioseqs(obskey, obs, sequence, threepoint=False):
     n2 = obskey[1, obskey[0, :] == 2]
 
     if (len(f0) == 0) or (len(f0) != (n0[-1] - n0[0] + 1)):
-        r01 = False
-        r10 = False
+        r01 = None
+        r10 = None
 
     if (len(f1) == 0) or (len(f1) != (n1[-1] - n1[0] + 1)):
-        r01 = False
-        r10 = False
+        r01 = None
+        r10 = None
 
     if (len(f2) == 0) or (len(f2) != (n2[-1] - n2[0] + 1)):
-        r02 = False
+        r02 = None
 
     # Two-point frequency ratio R02
     # -----------------------------
@@ -886,18 +888,24 @@ def compute_ratioseqs(obskey, obs, sequence, threepoint=False):
         return r10
 
     elif sequence == "r012":
+        if r01 is None or r02 is None:
+            return None
         # R012 (R01 followed by R02) ordered by n (R01 first for identical n)
         mask = np.argsort(np.append(r01[3, :], r02[3, :] + 0.1))
         r012 = np.hstack((r01, r02))[:, mask]
         return r012
 
     elif sequence == "r102":
+        if r10 is None or r02 is None:
+            return None
         # R102 (R10 followed by R02) ordered by n (R10 first for identical n)
         mask = np.argsort(np.append(r10[3, :], r02[3, :] + 0.1))
         r102 = np.hstack((r10, r02))[:, mask]
         return r102
 
     elif sequence == "r010":
+        if r01 is None or r10 is None:
+            return None
         # R010 (R01 followed by R10) ordered by n (R01 first for identical n)
         mask = np.argsort(np.append(r01[3, :], r10[3, :] + 0.1))
         r010 = np.hstack((r01, r10))[:, mask]
@@ -969,8 +977,6 @@ def compute_epsilondiff(
         Array containing the modes in the observed data.
     epsdiff_cov : array
         Covariances matrix.
-    epsdiff_covinv : array
-        The inverse of the covariance matrix.
     """
 
     # Remove modes outside of l=0 range
@@ -995,7 +1001,7 @@ def compute_epsilondiff(
     epsdiff = compute_epsilondiffseqs(
         osckey, osc, avgdnu, sequence=sequence, nsorting=nsorting
     )
-    epsdiff_cov, epsdiff_covinv = su.compute_cov_from_mc(
+    epsdiff_cov = su.compute_cov_from_mc(
         epsdiff.shape[1],
         osckey,
         osc,
@@ -1004,7 +1010,7 @@ def compute_epsilondiff(
         nrealisations=nrealisations,
     )
 
-    return epsdiff, epsdiff_cov, epsdiff_covinv
+    return epsdiff, epsdiff_cov
 
 
 def compute_epsilondiffseqs(
diff --git a/basta/interpolation_across.py b/basta/interpolation_across.py
index 5de7b83..d2c9ff7 100644
--- a/basta/interpolation_across.py
+++ b/basta/interpolation_across.py
@@ -15,6 +15,9 @@
 from basta import interpolation_helpers as ih
 from basta import plot_interp as ip
 
+import traceback
+
+
 # ======================================================================================
 # Interpolation helper routines
 # ======================================================================================
@@ -75,8 +78,8 @@ def _check_sobol(grid, res):
     return sobol
 
 
-def _calc_across_points(
-    base, baseparams, tri, sobol, outbasename, debug=False, verbose=False
+def _calc_cartesian_points(
+    base, baseparams, tri, outbasename, debug=False, verbose=False
 ):
     """
     Determine the new points for tracks in the base parameters, for either Cartesian
@@ -100,9 +103,6 @@ def _calc_across_points(
     tri : object
         Triangulation of the base.
 
-    sobol : float/bool
-        Scale resolution for Sobol-sampled interpolation, False for Cartesian.
-
     outbasename : str
         Name of the outputted plot of the base.
 
@@ -117,70 +117,125 @@ def _calc_across_points(
 
     """
 
-    if not sobol:
-        # Cartesian routine. Stores arrays of points to be added and all points
-        newbase = None
-        wholebase = copy.deepcopy(base)
-
-        # For each interpolation parameter, add the desired number of points
-        for i, (par, res) in enumerate(baseparams.items()):
-            newpoints = None
-            # Unique values of the parameter
-            uniq = np.unique(wholebase[:, i])
-            # New spacing in parameter
-            diff = np.mean(np.diff(uniq)) / (res + 1)
-            # For each requested new point, add an offsetted copy of the base
-            for j in range(res):
-                points = wholebase.copy()
-                points[:, i] += diff * (j + 1)
-                if type(newpoints) != np.ndarray:
-                    newpoints = points
-                else:
-                    newpoints = np.vstack((newpoints, points))
-            # Update the arrays
-            wholebase = np.vstack((wholebase, newpoints))
-            if type(newbase) != np.ndarray:
-                newbase = newpoints
+    # Stores arrays of points to be added and all points
+    newbase = None
+    wholebase = copy.deepcopy(base)
+
+    # For each interpolation parameter, add the desired number of points
+    for i, (par, res) in enumerate(baseparams.items()):
+        newpoints = None
+        # Unique values of the parameter
+        uniq = np.unique(wholebase[:, i])
+        # New spacing in parameter
+        diff = np.mean(np.diff(uniq)) / (res + 1)
+        # For each requested new point, add an offsetted copy of the base
+        for j in range(res):
+            points = wholebase.copy()
+            points[:, i] += diff * (j + 1)
+            if type(newpoints) != np.ndarray:
+                newpoints = points
             else:
-                newbase = np.vstack((newbase, newpoints))
+                newpoints = np.vstack((newpoints, points))
+        # Update the arrays
+        wholebase = np.vstack((wholebase, newpoints))
+        if type(newbase) != np.ndarray:
+            newbase = newpoints
+        else:
+            newbase = np.vstack((newbase, newpoints))
+
+    # Find all points within triangulation
+    mask = tri.find_simplex(newbase)
+    newbase = newbase[mask != -1]
+    trindex = tri.find_simplex(newbase)
+
+    # Plot of old vs. new base of subgrid
+    if len(baseparams) > 1 and (debug or verbose):
+        outname = outbasename.split(".")[-2] + "_all"
+        outname += "." + outbasename.split(".")[-1]
+        success = ip.base_corner(baseparams, base, newbase, tri, False, outname)
+        if success:
+            print(
+                "Initial across interpolation base has been plotted in",
+                "figure",
+                outname,
+            )
+    return newbase, trindex
+
+
+def _calc_sobol_points(
+    base, baseparams, tri, sobol, outbasename, debug=False, verbose=False
+):
+    """
+    Determine the new points for tracks in the base parameters, for Sobol sampling.
+    It determines a new Sobol sampling which satisfy an increase in number of tracks,
+    given a scale value.
+
+    Also plots the old vs new base of the interpolation, no plot for dim(base) = 1,
+    corner plot for dim(base) > 2.
+
+    Parameters
+    ----------
+    base : array
+        The current base of the grid, formed as (number of tracks, parameters in base).
+
+    baseparams : dict
+        Dictionary of the parameters forming the grid, with the required resolution of
+        the parameters.
+
+    tri : object
+        Triangulation of the base.
+
+    sobol : float
+        Scale resolution for Sobol-sampled interpolation
+
+    outbasename : str
+        Name of the outputted plot of the base.
+
+    Returns
+    -------
+    newbase : array
+        A base of the new points in the base, same structure as input base.
+
+    trindex : array
+        List of simplexes of the new points, for determination of the enveloping
+        tracks.
+
+    sob_nums : array
+        Sobol numbers used to generate new base, which is needed for determining
+        volume weights of the tracks
+    """
 
-        # Find all points within triangulation
+    # Check that we increase the number of tracks
+    lorgbase = len(base)
+    lnewbase = int(lorgbase * sobol)
+    assert lnewbase > lorgbase
+    ndim = len(baseparams)
+    l_trim = 1
+
+    # Try sampling the parameter space, and retry until increase met
+    while l_trim / sobol < lorgbase:
+        # Extract Sobol sequences
+        lnewbase = int(lnewbase * 1.2)
+        sob_nums = np.zeros((lnewbase, ndim))
+        iseed = 1
+        for i in range(lnewbase):
+            iseed, sob_nums[i, :] = sobol_numbers.i8_sobol(ndim, iseed)
+
+        # Assign parameter values by sequence
+        newbase = []
+        for npar in range(ndim):
+            Cmin = min(base[:, npar])
+            Cmax = max(base[:, npar])
+            newbase.append((Cmax - Cmin) * sob_nums[:, npar] + Cmin)
+
+        # Remove points outside subgrid
+        newbase = np.asarray(newbase).T
         mask = tri.find_simplex(newbase)
         newbase = newbase[mask != -1]
-        trindex = tri.find_simplex(newbase)
-
-    elif sobol:
-        # Check that we increase the number of tracks
-        lorgbase = len(base)
-        lnewbase = int(lorgbase * sobol)
-        assert lnewbase > lorgbase
-        ndim = len(baseparams)
-        l_trim = 1
-
-        # Try sampling the parameter space, and retry until increase met
-        while l_trim / sobol < lorgbase:
-            # Extract Sobol sequences
-            lnewbase = int(lnewbase * 1.2)
-            sob_nums = np.zeros((lnewbase, ndim))
-            iseed = 1
-            for i in range(lnewbase):
-                iseed, sob_nums[i, :] = sobol_numbers.i8_sobol(ndim, iseed)
-
-            # Assign parameter values by sequence
-            newbase = []
-            for npar in range(ndim):
-                Cmin = min(base[:, npar])
-                Cmax = max(base[:, npar])
-                newbase.append((Cmax - Cmin) * sob_nums[:, npar] + Cmin)
-
-            # Remove points outside subgrid
-            newbase = np.asarray(newbase).T
-            mask = tri.find_simplex(newbase)
-            newbase = newbase[mask != -1]
-            l_trim = len(newbase[:, 0])
-
-        # Compute new simplex list
-        trindex = tri.find_simplex(newbase)
+        l_trim = len(newbase[:, 0])
+
+    # Compute new simplex list
+    trindex = tri.find_simplex(newbase)
 
     # Plot of old vs. new base of subgrid
     if len(baseparams) > 1 and (debug or verbose):
@@ -193,7 +248,7 @@ def _calc_across_points(
                 "figure",
                 outname,
             )
-    return newbase, trindex
+    return newbase, trindex, sob_nums[mask != -1]
 
 
 # ======================================================================================
@@ -299,7 +354,7 @@ def _interpolate_across(
         numfmt = len(tracklist[0][0].split("track")[-1])
 
         # Form basis of varied parameters
-        bpars = [par.decode("UTF-8") for par in grid["header/active_weights"]]
+        bpars = [par.decode("UTF-8") for par in grid["header/pars_sampled"]]
         baseparams = {par: resolution[par] for par in bpars}
         const_vars = {}
         for par in headvars:
@@ -317,7 +372,7 @@ def _interpolate_across(
         newnum = 0
 
         # Parameters for forming basis
-        bpars = [par.decode("UTF-8") for par in grid["header/active_weights"]]
+        bpars = [par.decode("UTF-8") for par in grid["header/pars_sampled"]]
         baseparams = {par: resolution[par] for par in bpars}
         const_vars = {}
         isochhead = os.path.join("header", basepath)
@@ -344,7 +399,7 @@ def _interpolate_across(
             if not isomode:
                 index2d = np.array(np.transpose([index, index]))
             if not (any(index) and sum(index) > 2):
-                outfile[os.path.join(name, "FeHini_weight")] = -1
+                outfile[os.path.join(name, "IntStatus")] = -1
             else:
                 # Write everything from the old grid to the new in the region
                 for key in grid[name].keys():
@@ -367,9 +422,14 @@ def _interpolate_across(
     # Determine the base params for new tracks
     print("\nBuilding triangulation ... ", end="", flush=True)
     triangulation = spatial.Delaunay(base)
-    new_points, trindex = _calc_across_points(
-        base, baseparams, triangulation, sobol, outbasename, debug, verbose
-    )
+    if sobol:
+        new_points, trindex, sobnums = _calc_sobol_points(
+            base, baseparams, triangulation, sobol, outbasename, debug, verbose
+        )
+    else:
+        new_points, trindex = _calc_cartesian_points(
+            base, baseparams, triangulation, outbasename, debug, verbose
+        )
     print("done!")
 
     # List of tracknames for accessing grid
@@ -433,7 +493,7 @@ def _interpolate_across(
             warstr += " of the enveloping {0}!".format(modestr)
             print(warstr)
             success[tracknum] = False
-            outfile[os.path.join(libname, "FeHini_weight")] = -1
+            outfile[os.path.join(libname, "IntStatus")] = -1
             continue
 
         # Assume equal spacing, but approximately the same number of points
@@ -510,19 +570,18 @@ def _interpolate_across(
                 dsetosc = outfile.create_dataset(
                     name=os.path.join(libname, "osc"),
                     shape=(Npoints, 2),
-                    dtype=h5py.special_dtype(vlen=np.float),
+                    dtype=h5py.special_dtype(vlen=float),
                 )
                 dsetosckey = outfile.create_dataset(
                     name=os.path.join(libname, "osckey"),
                     shape=(Npoints, 2),
-                    dtype=h5py.special_dtype(vlen=np.int),
+                    dtype=h5py.special_dtype(vlen=int),
                 )
                 for i in range(Npoints):
                     dsetosc[i] = newosc[i]
                     dsetosckey[i] = newosckey[i]
 
             # Dealing with constants of the track
-            outfile[os.path.join(libname, "FeHini_weight")] = 1.0
             for par, parval in zip(baseparams, point):
                 keypath = os.path.join(libname, par)
                 try:
@@ -544,6 +603,9 @@ def _interpolate_across(
             keypath = os.path.join(libname, dname)
             outfile[keypath] = ih.bay_weights(outfile[parpath])
 
+            # Successfully interpolated, mark it as such
+            outfile[os.path.join(libname, "IntStatus")] = 0
+
             if debug:
                 debugnum = str(int(newnum + tracknum)).zfill(numfmt)
                 plotpath = os.path.join(
@@ -584,7 +646,7 @@ def _interpolate_across(
                 None
             success[tracknum] = False
             print("Error:", sys.exc_info()[1])
-            outfile[os.path.join(libname, "FeHini_weight")] = -1
+            outfile[os.path.join(libname, "IntStatus")] = -1
             print("Interpolation failed for {0}".format(libname))
             if debug:
                 print("Point at:")
@@ -616,5 +678,8 @@ def _interpolate_across(
 
     # Write the new tracks to the header, and recalculate the weights
     outfile = ih._extend_header(outfile, basepath, headvars)
-    outfile = ih._recalculate_weights(outfile, basepath, headvars)
+    if "volume" in grid["header/active_weights"] or sobol:
+        outfile = ih._recalculate_volume_weights(outfile, basepath, sobnums)
+    else:
+        outfile = ih._recalculate_param_weights(outfile, basepath)
     return grid, outfile, fail
diff --git a/basta/interpolation_along.py b/basta/interpolation_along.py
index 275eec6..2be53b1 100644
--- a/basta/interpolation_along.py
+++ b/basta/interpolation_along.py
@@ -240,6 +240,14 @@ def _interpolate_along(
 
     overwrite = grid == outfile
 
+    # Check if grid is interpolated
+    try:
+        grid["header/interpolation_time"][()]
+    except KeyError:
+        grid_is_intpol = False
+    else:
+        grid_is_intpol = True
+
     if debug:
         verbose = True
     if verbose:
@@ -317,8 +325,11 @@ def _interpolate_along(
             # Update progress bar in the start of the loop to count skipped tracks
             pbar.update(1)
 
-            if libitem["FeHini_weight"][()] < 0:
-                continue
+            # Skip failed input tracks
+            if grid_is_intpol:
+                if libitem["IntStatus"][()] < 0:
+                    continue
+
             # Make sure the user provides a valid parameter as resolution requirement
             # --> In case of failure, assume the user wanted a dnu-related parameter
             if resolution["param"].lower() not in freqres:
@@ -468,28 +479,31 @@ def _interpolate_along(
                     dsetosc = outfile.create_dataset(
                         name=os.path.join(libitem.name, "osc"),
                         shape=(Npoints, 2),
-                        dtype=h5py.special_dtype(vlen=np.float),
+                        dtype=h5py.special_dtype(vlen=float),
                     )
                     dsetosckey = outfile.create_dataset(
                         name=os.path.join(libitem.name, "osckey"),
                         shape=(Npoints, 2),
-                        dtype=h5py.special_dtype(vlen=np.int),
+                        dtype=h5py.special_dtype(vlen=int),
                     )
                     for i in range(Npoints):
                         dsetosc[i] = osclist[i]
                         dsetosckey[i] = osckeylist[i]
 
+                # Successfully interpolated, mark it as such
+                if "IntStatus" in outfile[libitem.name]:
+                    del outfile[os.path.join(libitem.name, "IntStatus")]
+                outfile[os.path.join(libitem.name, "IntStatus")] = 0
                 trackcounter += 1
 
             #
             # *** BLOCK 1b: Handle tracks without any models inside selection ***
             #
             else:
-                # Flag the empty tracks with a single entry: A negative weight
-                # (note: requires at least BASTA v0.28 to be useful)
-                if overwrite:
-                    del outfile[os.path.join(libitem.name, "FeHini_weight")]
-                outfile[os.path.join(libitem.name, "FeHini_weight")] = -1
+                # Flag the empty track as "failed", to avoid errors in the main BASTA loop
+                if "IntStatus" in outfile[libitem.name]:
+                    del outfile[os.path.join(libitem.name, "IntStatus")]
+                outfile[os.path.join(libitem.name, "IntStatus")] = -1
                 if verbose:
                     print()
 
diff --git a/basta/interpolation_helpers.py b/basta/interpolation_helpers.py
index 4f495ec..2d3b0b6 100644
--- a/basta/interpolation_helpers.py
+++ b/basta/interpolation_helpers.py
@@ -149,7 +149,10 @@ def get_selectedmodels(grid, basepath, limits, cut=True, show_progress=True):
             if show_progress:
                 pbar.update(1)
             # If empty, skip
-            if not len(libitem) or libitem["FeHini_weight"][()] == -1:
+            if not len(libitem):
+                continue
+            # If previously interpolated, and track failed
+            if "IntStatus" in libitem and libitem["IntStatus"][()] < 0:
                 continue
 
             # Full list of indexes, set False if model outside limits
@@ -493,8 +496,8 @@ def interpolate_frequencies(
             ]
 
         # Add the stacked ".obs-style" lists with all freq information to storage arrays
-        osckeylist.append(np.array([lf, nf], dtype=np.int))
-        osclist.append(np.array([ff, ef], dtype=np.float))
+        osckeylist.append(np.array([lf, nf], dtype=int))
+        osclist.append(np.array([ff, ef], dtype=float))
 
     return osckeylist, osclist
 
@@ -538,7 +541,7 @@ def _extend_header(outfile, basepath, headvars):
             values = np.zeros(ltracks)
             for _, group in outfile[basepath].items():
                 for n, (_, libitem) in enumerate(group.items()):
-                    if libitem["FeHini_weight"][()] != -1:
+                    if libitem["IntStatus"][()] != -1:
                         values[n] = libitem[var][0]
             del outfile[headpath]
             outfile[headpath] = values.tolist()
@@ -603,10 +606,10 @@ def _write_header(grid, outfile, basepath):
     return grid, outfile
 
 
-def _recalculate_weights(outfile, basepath, headvars):
+def _recalculate_param_weights(outfile, basepath):
     """
     Recalculates the weights of the tracks/isochrones, for the new grid.
-    Tracks not transferred from old grid has FeHini_weight = -1.
+    Tracks not transferred from old grid has IntStatus = -1.
 
     Parameters
     ----------
@@ -627,6 +630,7 @@ def _recalculate_weights(outfile, basepath, headvars):
 
     """
     isomode = False if "grid" in basepath else True
+    headvars = outfile["header/active_weights"][()]
 
     # Collect the relevant tracks/isochrones
     mask = []
@@ -634,9 +638,9 @@ def _recalculate_weights(outfile, basepath, headvars):
     for nogroup, (gname, group) in enumerate(outfile[basepath].items()):
         # Determine which tracks are actually present
         for name, libitem in group.items():
-            mask.append(libitem["FeHini_weight"][()])
+            mask.append(libitem["IntStatus"][()])
             names.append(os.path.join(gname, name))
-    mask = np.where(np.array(mask) > 0)[0]
+    mask = np.where(np.array(mask) >= 0)[0]
     active = np.asarray(names)[mask]
 
     # For each parameter, collect values, recalculate weights, and replace old weight
@@ -659,3 +663,86 @@ def _recalculate_weights(outfile, basepath, headvars):
                 del outfile[weight_path]
                 outfile[weight_path] = weights[i]
     return outfile
+
+
+def _recalculate_volume_weights(outfile, basepath, sobnums):
+    """
+    Recalculates the weights of the tracks/isochrones, for the new grid.
+    Tracks not transferred from old grid has IntStatus = -1.
+
+    Parameters
+    ----------
+    outfile : hdf5 file
+        New grid file to write to.
+
+    basepath : str
+        Path in the grid where the tracks are stored. The default value given in
+        parent functions applies to standard grids of tracks.
+
+    sobnums : array
+        Sobol numbers used to generate new base
+
+    extend : bool
+        Whether the old tracks have been preserved, for which the Sobol numbers
+        needs to be recovered.
+
+    Returns
+    -------
+    outfile : hdf5 file
+        New grid file to write to.
+
+    """
+    # Collect the relevant tracks/isochrones
+    IntStatus = []
+    names = []
+    for gname, group in outfile[basepath].items():
+        # Determine which tracks are actually present
+        for name, libitem in group.items():
+            IntStatus.append(libitem["IntStatus"][()])
+            names.append(os.path.join(gname, name))
+    mask = np.where(np.array(IntStatus) >= 0)[0]
+    active = np.asarray(names)[mask]
+
+    # Import necessary packages only used here
+    import bottleneck as bn
+    from basta import sobol_numbers
+
+    # Use IntStatus mask, read key numbers
+    sobnums = sobnums[mask, :]
+    ntracks = sobnums.shape[0]
+    ndim = sobnums.shape[1]
+
+    # Generate oversampled grid
+    iseed = 2
+    osfactor = 100
+    osntracks = osfactor * ntracks
+    osbase = np.zeros((osntracks, ndim))
+    for i in range(osntracks):
+        iseed, osbase[i, :] = sobol_numbers.i8_sobol(ndim, iseed)
+
+    # For every track in the grid, gather all the points from the
+    # oversampled grid which are closest to the track at hand
+    npoints = np.zeros(ntracks, dtype=int)
+    for i in range(osntracks):
+        diff = osbase[i, :] - sobnums
+        distance2 = bn.nansum(diff**2, axis=-1)
+        nclosest = bn.nanargmin(distance2)
+        npoints[nclosest] += 1
+
+    # Transform to volume weight per track
+    weights = npoints / bn.nansum(npoints)
+
+    # Write the weights
+    for i, name in enumerate(active):
+        weight_path = os.path.join(basepath, name, "volume_weight")
+        try:
+            outfile[weight_path] = weights[i]
+        except:
+            del outfile[weight_path]
+            outfile[weight_path] = weights[i]
+
+    # Write the active weights as only volume
+    del outfile["header/active_weights"]
+    outfile["header/active_weights"] = ["volume"]
+
+    return outfile
diff --git a/basta/plot_corner.py b/basta/plot_corner.py
index 0a4bc92..8beee69 100644
--- a/basta/plot_corner.py
+++ b/basta/plot_corner.py
@@ -188,6 +188,13 @@ def corner(
         for i, x in enumerate(xs):
             try:
                 xbin = np.histogram_bin_edges(x, bins="auto", range=np.sort(prange[i]))
+                if len(xbin) > 1000:
+                    print(
+                        "Parameter {0} resulted in {1} bins, raising MemoryError".format(
+                            labels[i], len(xbin)
+                        )
+                    )
+                    raise MemoryError
             except MemoryError:
                 print(
                     "WARNING! Using 'auto' as bin-rule causes a memory crash!"
@@ -353,6 +360,7 @@ def corner(
                     ax.set_title(labels[i], y=1.25, **label_kwargs)
                 else:
                     ax.set_xlabel(labels[i], **label_kwargs)
+                    ax.xaxis.set_label_coords(0.5, -0.35)
 
             # use MathText for axes ticks
             ax.xaxis.set_major_formatter(ScalarFormatter(useMathText=use_math_text))
diff --git a/basta/plot_driver.py b/basta/plot_driver.py
new file mode 100644
index 0000000..87b948f
--- /dev/null
+++ b/basta/plot_driver.py
@@ -0,0 +1,342 @@
+import numpy as np
+
+from basta import freq_fit, plot_seismic
+from basta import utils_seismic as su
+from basta.constants import freqtypes
+
+
+def plot_all_seismic(
+    freqplots,
+    Grid,
+    fitfreqs,
+    obsfreqmeta,
+    obsfreqdata,
+    obskey,
+    obs,
+    obsintervals,
+    selectedmodels,
+    path,
+    ind,
+    plotfname,
+    debug=False,
+):
+    """
+    Driver for producing all seismic related plots
+
+    Parameters
+    ----------
+    freqplots : list
+        Plots to be produced
+    Grid : hdf5 file
+        Stellar models, as tracks or isochrones
+    fitfreqs : dict
+        Input frequency fitting options/controls
+    obsfreqmeta : dict
+        The requested information about which frequency products to fit or
+        plot, unpacked for easier access later.
+    obsfreqdata : dict
+        Requested frequency-dependent data such as glitches, ratios, and
+        epsilon difference. It also contains the covariance matrix and its
+        inverse of the individual frequency modes.
+        The keys correspond to the science case, e.g. `r01`, `glitch`, or
+        `e012`.
+        Inside each case, you find the data (`data`), the covariance matrix
+        (`cov`), and its inverse (`covinv`).
+    obskey : array
+        Array containing the angular degrees and radial orders of obs
+    obs : array
+        Individual frequencies and uncertainties.
+    obsintervals : array
+        Array containing the endpoints of the intervals used in the frequency
+        fitting routine in :func:'freq_fit.calc_join'.
+    selectedmodels : dict
+        Contains information on all models with a non-zero likelihood.
+    path : str
+        Path to the highest likelihood track/isocohrone in the grid
+    ind : int
+        Index of the highest likelihood model in the track
+    plotfname : str
+        Output plotname format
+    debug : bool
+        Whether to produce debugging output
+
+    """
+
+    # Check which plots to create
+    allfplots = freqplots[0] == True
+    if any(x == "allechelle" for x in freqplots):
+        freqplots += ["dupechelle", "echelle", "pairechelle"]
+    if any(x in freqtypes.rtypes for x in freqplots):
+        freqplots += ["ratios"]
+    try:
+        rawmaxmod = Grid[path + "/osc"][ind]
+        rawmaxmodkey = Grid[path + "/osckey"][ind]
+        maxmod = su.transform_obj_array(rawmaxmod)
+        maxmodkey = su.transform_obj_array(rawmaxmodkey)
+        maxmod = maxmod[:, maxmodkey[0, :] < 2.5]
+        maxmodkey = maxmodkey[:, maxmodkey[0, :] < 2.5]
+        maxjoins = freq_fit.calc_join(
+            mod=maxmod,
+            modkey=maxmodkey,
+            obs=obs,
+            obskey=obskey,
+            obsintervals=obsintervals,
+        )
+        maxjoinkeys, maxjoin = maxjoins
+        maxmoddnu = Grid[path + "/dnufit"][ind]
+    except Exception as e:
+        print("\nFrequency plots initialisation failed with the error:", e)
+        return None
+
+    # Extract the original observed dnu for use on the echelle diagrams
+    # --> (equivalent to re-scaling if solar scaling activated)
+    plotdnu = fitfreqs["dnu_obs"]
+
+    if allfplots or "echelle" in freqplots:
+        try:
+            plot_seismic.echelle(
+                selectedmodels,
+                Grid,
+                obs,
+                obskey,
+                mod=maxmod,
+                modkey=maxmodkey,
+                dnu=plotdnu,
+                join=maxjoin,
+                joinkeys=maxjoinkeys,
+                pair=False,
+                duplicate=False,
+                output=plotfname.format("echelle_uncorrected"),
+            )
+        except Exception as e:
+            print("\nUncorrected echelle failed with the error:", e)
+
+    if allfplots or "pairechelle" in freqplots:
+        try:
+            plot_seismic.echelle(
+                selectedmodels,
+                Grid,
+                obs,
+                obskey,
+                mod=maxmod,
+                modkey=maxmodkey,
+                dnu=plotdnu,
+                join=maxjoin,
+                joinkeys=maxjoinkeys,
+                pair=True,
+                duplicate=False,
+                output=plotfname.format("pairechelle_uncorrected"),
+            )
+        except Exception as e:
+            print("\nUncorrected pairechelle failed with the error:", e)
+
+    if allfplots or "dupechelle" in freqplots:
+        try:
+            plot_seismic.echelle(
+                selectedmodels,
+                Grid,
+                obs,
+                obskey,
+                mod=maxmod,
+                modkey=maxmodkey,
+                dnu=plotdnu,
+                join=maxjoin,
+                joinkeys=maxjoinkeys,
+                duplicate=True,
+                pair=True,
+                output=plotfname.format("dupechelle_uncorrected"),
+            )
+        except Exception as e:
+            print("\nUncorrected dupechelle failed with the error:", e)
+
+    if fitfreqs["fcor"] == "None":
+        corjoin = maxjoin
+        coeffs = [1]
+    elif fitfreqs["fcor"] == "HK08":
+        corjoin, coeffs = freq_fit.HK08(
+            joinkeys=maxjoinkeys,
+            join=maxjoin,
+            nuref=fitfreqs["numax"],
+            bcor=fitfreqs["bexp"],
+        )
+    elif fitfreqs["fcor"] == "BG14":
+        corjoin, coeffs = freq_fit.BG14(
+            joinkeys=maxjoinkeys, join=maxjoin, scalnu=fitfreqs["numax"]
+        )
+    elif fitfreqs["fcor"] == "cubicBG14":
+        corjoin, coeffs = freq_fit.cubicBG14(
+            joinkeys=maxjoinkeys, join=maxjoin, scalnu=fitfreqs["numax"]
+        )
+
+    if len(coeffs) > 1:
+        print("The surface correction coefficients are", *coeffs)
+    else:
+        print("The surface correction coefficient is", *coeffs)
+
+    if allfplots or "echelle" in freqplots:
+        try:
+            plot_seismic.echelle(
+                selectedmodels,
+                Grid,
+                obs,
+                obskey,
+                mod=maxmod,
+                modkey=maxmodkey,
+                dnu=plotdnu,
+                join=corjoin,
+                joinkeys=maxjoinkeys,
+                freqcor=fitfreqs["fcor"],
+                coeffs=coeffs,
+                scalnu=fitfreqs["numax"],
+                pair=False,
+                duplicate=False,
+                output=plotfname.format("echelle"),
+            )
+        except Exception as e:
+            print("\nEchelle failed with the error:", e)
+
+    if allfplots or "pairechelle" in freqplots:
+        try:
+            plot_seismic.echelle(
+                selectedmodels,
+                Grid,
+                obs,
+                obskey,
+                mod=maxmod,
+                modkey=maxmodkey,
+                dnu=plotdnu,
+                join=corjoin,
+                joinkeys=maxjoinkeys,
+                freqcor=fitfreqs["fcor"],
+                coeffs=coeffs,
+                scalnu=fitfreqs["numax"],
+                pair=True,
+                duplicate=False,
+                output=plotfname.format("pairechelle"),
+            )
+        except Exception as e:
+            print("\nPairechelle failed with the error:", e)
+
+    if allfplots or "dupechelle" in freqplots:
+        try:
+            plot_seismic.echelle(
+                selectedmodels,
+                Grid,
+                obs,
+                obskey,
+                mod=maxmod,
+                modkey=maxmodkey,
+                dnu=plotdnu,
+                join=corjoin,
+                joinkeys=maxjoinkeys,
+                freqcor=fitfreqs["fcor"],
+                coeffs=coeffs,
+                scalnu=fitfreqs["numax"],
+                duplicate=True,
+                pair=True,
+                output=plotfname.format("dupechelle"),
+            )
+        except Exception as e:
+            print("\nDupechelle failed with the error:", e)
+
+    if "freqcormap" in freqplots or debug:
+        try:
+            plot_seismic.correlation_map(
+                "freqs",
+                obsfreqdata,
+                plotfname.format("freqs_cormap"),
+                obskey=obskey,
+            )
+        except Exception as e:
+            print("\nFrequencies correlation map failed with the error:", e)
+
+    if obsfreqmeta["getratios"]:
+        for ratseq in obsfreqmeta["ratios"]["plot"]:
+            try:
+                ratnamestr = "ratios_{0}".format(ratseq)
+                plot_seismic.ratioplot(
+                    obsfreqdata,
+                    maxjoinkeys,
+                    maxjoin,
+                    maxmodkey,
+                    maxmod,
+                    ratseq,
+                    output=plotfname.format(ratnamestr),
+                    threepoint=fitfreqs["threepoint"],
+                    interp_ratios=fitfreqs["interp_ratios"],
+                )
+            except Exception as e:
+                print(
+                    "\nRatio plot for {} sequence failed with the error:".format(
+                        ratseq
+                    ),
+                    e,
+                )
+
+            if fitfreqs["correlations"]:
+                try:
+                    plot_seismic.correlation_map(
+                        ratseq,
+                        obsfreqdata,
+                        output=plotfname.format(ratnamestr + "_cormap"),
+                    )
+                except Exception as e:
+                    print(
+                        "\nRatio correlation map for {} sequence failed with the error:".format(
+                            ratseq
+                        ),
+                        e,
+                    )
+
+    if obsfreqmeta["getepsdiff"]:
+        for epsseq in obsfreqmeta["epsdiff"]["plot"]:
+            try:
+                epsnamestr = "epsdiff_{0}".format(epsseq)
+                plot_seismic.epsilon_difference_diagram(
+                    mod=maxmod,
+                    modkey=maxmodkey,
+                    moddnu=maxmoddnu,
+                    sequence=epsseq,
+                    obsfreqdata=obsfreqdata,
+                    output=plotfname.format(epsnamestr),
+                )
+            except Exception as e:
+                print(
+                    "\nEpsilon difference plot for {} sequence failed with the error:".format(
+                        epsseq
+                    ),
+                    e,
+                )
+
+            if fitfreqs["correlations"]:
+                try:
+                    plot_seismic.correlation_map(
+                        epsseq,
+                        obsfreqdata,
+                        output=plotfname.format(epsnamestr + "_cormap"),
+                    )
+                except Exception as e:
+                    print(
+                        "\nEpsilon difference correlation map for {} sequence failed with the error:".format(
+                            epsseq
+                        ),
+                        e,
+                    )
+
+    if obsfreqmeta["getepsdiff"] and debug:
+        if len(obsfreqmeta["epsdiff"]["plot"]) > 0:
+            try:
+                plot_seismic.epsilon_difference_components_diagram(
+                    mod=maxmod,
+                    modkey=maxmodkey,
+                    moddnu=maxmoddnu,
+                    obs=obs,
+                    obskey=obskey,
+                    dnudata=obsfreqdata["freqs"]["dnudata"],
+                    obsfreqdata=obsfreqdata,
+                    obsfreqmeta=obsfreqmeta,
+                    output=plotfname.format("DEBUG_epsdiff_components"),
+                )
+            except Exception as e:
+                print("\nEpsilon difference compoenent plot failed with the error:", e)
+    return None
diff --git a/basta/plot_kiel.py b/basta/plot_kiel.py
index 5568ec2..b823222 100644
--- a/basta/plot_kiel.py
+++ b/basta/plot_kiel.py
@@ -312,14 +312,13 @@ def kiel(
         )
         iteration = zip(axis, tefflim, logglim)
     else:
-        fig, axis = plt.subplots(1, 1)
+        fig, axis = plt.subplots(1, 1, figsize=(8.47, 6))
         iteration = zip([axis], tefflim, logglim)
 
     # Iteration over the subplots, the same is done with different limits
     for ax, tlim, glim in iteration:
         max_logPDF = selectedmodels[maxPDF_path].logPDF.max()
         for track in tracks:
-
             # Make a copy to allow manipulation
             xs = gu.h5py_to_array(Grid[track + "/Teff"])
             ys = gu.h5py_to_array(Grid[track + "/logg"])
diff --git a/basta/plot_seismic.py b/basta/plot_seismic.py
index 550fe4c..2d86e0c 100644
--- a/basta/plot_seismic.py
+++ b/basta/plot_seismic.py
@@ -308,11 +308,13 @@ def echelle(
                         if ((fm[i] % modx) > linelimit) & (
                             (fo[i] % modx) < (modx - linelimit)
                         ):
-                            a = (fm[i] - fo[i]) / ((fm[i] - fo[i]) % modx)
                             if duplicate:
                                 x0 = -1
                             else:
                                 x0 = 0
+                            a = (fo[i] - fm[i]) / abs(
+                                fo[i] % modx - (fm[i] % modx - modx)
+                            )
                             ax.plot(
                                 (fm[i] % modx - modx, fo[i] % modx),
                                 (fm[i], fo[i]),
@@ -322,14 +324,14 @@ def echelle(
                             )
                             ax.plot(
                                 (fm[i] % modx, modx),
-                                (fm[i], a * np.ceil(fm[i] / modx) * modx),
+                                (fm[i], fm[i] + a * (modx - fm[i] % modx)),
                                 c=colors["l" + str(l)],
                                 alpha=0.7,
                                 zorder=10,
                             )
                             ax.plot(
                                 (x0, fo[i] % modx - modx),
-                                (a * np.ceil(fm[i] / modx) * modx, fo[i]),
+                                (fm[i] + a * (modx - fm[i] % modx), fo[i]),
                                 c=colors["l" + str(l)],
                                 alpha=0.7,
                                 zorder=10,
@@ -341,7 +343,9 @@ def echelle(
                                 x0 = -1
                             else:
                                 x0 = 0
-                            a = (fm[i] - fo[i]) / ((fm[i] - fo[i]) % modx)
+                            a = (fm[i] - fo[i]) / abs(
+                                fm[i] % modx - (fo[i] % modx - modx)
+                            )
                             ax.plot(
                                 (fo[i] % modx - modx, fm[i] % modx),
                                 (fo[i], fm[i]),
@@ -351,14 +355,14 @@ def echelle(
                             )
                             ax.plot(
                                 (fo[i] % modx, modx),
-                                (fo[i], a * np.ceil(fo[i] / modx) * modx),
+                                (fo[i], fo[i] + a * (modx - fo[i] % modx)),
                                 c=colors["l" + str(l)],
                                 alpha=0.7,
                                 zorder=10,
                             )
                             ax.plot(
                                 (x0, fm[i] % modx - modx),
-                                (a * np.ceil(fo[i] / modx) * modx, fm[i]),
+                                (fo[i] + a * (modx - fo[i] % modx), fm[i]),
                                 c=colors["l" + str(l)],
                                 alpha=0.7,
                                 zorder=10,
@@ -429,6 +433,7 @@ def echelle(
     if output is not None:
         plt.savefig(output, bbox_inches="tight")
         print("Saved figure to " + output)
+        plt.close(fig)
 
 
 def ratioplot(
@@ -482,13 +487,17 @@ def ratioplot(
     fig, ax = plt.subplots(1, 1)
 
     obsratio = obsfreqdata[ratiotype]["data"]
+    # Exit if there are no ratios to plot
+    if obsratio is None:
+        plt.close(fig)
+        return
+
     obsratio_cov = obsfreqdata[ratiotype]["cov"]
     obsratio_err = np.sqrt(np.diag(obsratio_cov))
 
     if interp_ratios:
-        broad_key, broad_osc = su.extend_modjoin(joinkeys, join, modkey, mod)
         modratio = freq_fit.compute_ratioseqs(
-            broad_key, broad_osc, ratiotype, threepoint=threepoint
+            modkey, mod, ratiotype, threepoint=threepoint
         )
     else:
         modratio = freq_fit.compute_ratioseqs(
@@ -542,10 +551,14 @@ def ratioplot(
             intfunc = interp1d(
                 modratio[1, modmask], modratio[0, modmask], kind="linear"
             )
-            newmod = intfunc(obsratio[1, obsmask])
+            # When only plotting, not fitting, model freqs can be outside observed range
+            rangemask = np.ones(sum(obsmask), dtype=bool)
+            rangemask &= obsratio[1, obsmask] > min(modratio[1, modmask])
+            rangemask &= obsratio[1, obsmask] < max(modratio[1, modmask])
+            newmod = intfunc(obsratio[1, obsmask][rangemask])
             marker = splinemarkers[1] if "1" in str(rtype) else splinemarkers[2]
             (intp,) = ax.plot(
-                obsratio[1, obsmask],
+                obsratio[1, obsmask][rangemask],
                 newmod,
                 marker=marker,
                 color="k",
@@ -574,6 +587,8 @@ def ratioplot(
 
     ax.set_xlabel(r"Frequency ($\mu$Hz)")
     ax.set_ylabel(f"Frequency ratio ({ratiotype})")
+    ylim = ax.get_ylim()
+    ax.set_ylim(max(ylim[0], 0), ylim[1])
 
     if output is not None:
         fig.savefig(output, bbox_inches="tight")
@@ -711,11 +726,14 @@ def epsilon_difference_diagram(
 
     ax.set_xlabel(r"Frequency ($\mu$Hz)")
     ax.set_ylabel(r"Epsilon difference $\delta\epsilon_{0\ell}$")
+    ylim = ax.get_ylim()
+    ax.set_ylim(max(ylim[0], 0), ylim[1])
 
     fig.tight_layout()
     if output is not None:
         print("Saved figure to " + output)
         fig.savefig(output)
+        plt.close(fig)
     else:
         return fig
 
@@ -745,11 +763,15 @@ def correlation_map(fittype, obsfreqdata, output, obskey=None):
 
     elif fittype in freqtypes.rtypes:
         data = obsfreqdata[fittype]["data"]
+        if data is None:
+            return
         fmtstr = r"$r_{{{:02d}}}({{{:d}}})$"
         ln_zip = zip(data[2, :], data[3, :])
 
     elif fittype in freqtypes.epsdiff:
         data = obsfreqdata[fittype]["data"]
+        if data is None:
+            return
         fmtstr = r"$\delta\epsilon_{{{:02d}}}({{{:d}}})$"
         ln_zip = zip(data[2, :], data[3, :])
 
@@ -777,8 +799,8 @@ def correlation_map(fittype, obsfreqdata, output, obskey=None):
 
     if output is not None:
         fig.savefig(output, bbox_inches="tight")
-        plt.close(fig)
         print("Saved figure to " + output)
+        plt.close(fig)
 
 
 ###############
@@ -867,11 +889,10 @@ def epsilon_difference_components_diagram(
     expol = np.where(np.logical_or(nu12 < min(nu0), nu12 > max(nu0)))[0]
 
     # Definition of figure
-    fig, ax = plt.subplots(2, 2, figsize=(11, 9), sharex="col", sharey="row")
+    fig, ax = plt.subplots(2, 2, figsize=(8.47, 6), sharex="col", sharey="row")
 
     # Epsilon differences
     for ll in l_available:
-
         indobs = obsepsdiff[2] == ll
         indmod = modepsdiff[2] == ll
 
@@ -905,7 +926,7 @@ def epsilon_difference_components_diagram(
             obsepsdiff[1][indobs],
             obsepsdiff[0][indobs],
             yerr=obsepsdiff_err[indobs],
-            marker=obsmarker,
+            marker=modmarkers["l" + str(ll)],
             color=colors["l" + str(ll)],
             markeredgewidth=0.5,
             markeredgecolor="k",
@@ -944,7 +965,7 @@ def epsilon_difference_components_diagram(
             eps,
             yerr=err,
             linestyle=None,
-            fmt=obsmarker,
+            fmt=modmarkers["l" + str(ll)],
             color=colors["l" + str(ll)],
             zorder=3,
             label=elab % ll,
@@ -997,3 +1018,4 @@ def epsilon_difference_components_diagram(
 
     fig.savefig(output)
     print("Saved figure to " + output)
+    plt.close(fig)
diff --git a/basta/process_output.py b/basta/process_output.py
index 24167e6..9ee7763 100644
--- a/basta/process_output.py
+++ b/basta/process_output.py
@@ -11,7 +11,7 @@
 
 import basta.fileio as fio
 from basta.constants import sydsun as sydc
-from basta.constants import parameters
+from basta.constants import parameters, statdata
 from basta.utils_distances import distance_from_mag
 from basta.distances import get_absorption, LOS_reddening
 from basta import utils_general as util
@@ -92,17 +92,12 @@ def compute_posterior(
     hout_dist.append("starid")
     out_dist.append(starid)
 
-    # Settings for the quantiles and samples
-    qs = [0.5, 0.158655, 0.841345]
-    nsigma = 0.25
-    nsamples = 100000
-
     # Generate PDF values
     logy = np.concatenate([ts.logPDF for ts in selectedmodels.values()])
     noofind = len(logy)
     nonzeroprop = np.isfinite(logy)
     logy = logy[nonzeroprop]
-    nsamples = min(nsamples, noofind)
+    nsamples = min(statdata.nsamples, noofind)
 
     # Likelihood is only defined up to a multiplicative constant of
     # proportionality, therefore we subtract max(logy) from logy to make sure
@@ -124,7 +119,7 @@ def compute_posterior(
     # Corner plot kwargs
     ckwargs = {
         "show_titles": True,
-        "quantiles": qs,
+        "quantiles": statdata.quantiles,
         "smooth": 1,
         "smooth1d": "kde",
         "title_kwargs": {"fontsize": 10},
@@ -158,10 +153,14 @@ def compute_posterior(
 
             # Create posteriors from weighted histograms
             dinterp.append(
-                stats.posterior(d_all, nonzeroprop, sampled_indices, nsigma=nsigma)
+                stats.posterior(
+                    d_all, nonzeroprop, sampled_indices, nsigma=statdata.nsigma
+                )
             )
             EBVinterp.append(
-                stats.posterior(EBV_all, nonzeroprop, sampled_indices, nsigma=nsigma)
+                stats.posterior(
+                    EBV_all, nonzeroprop, sampled_indices, nsigma=statdata.nsigma
+                )
             )
 
             # Compute centroid and uncertainties and print them
@@ -469,16 +468,24 @@ def compute_posterior(
                 x = util.get_parameter_values(
                     library_param, Grid, selectedmodels, noofind
                 )
-                lp_interval = np.quantile(x[nonzeroprop][sampled_indices], qs[1:])
+                lp_interval = np.quantile(
+                    x[nonzeroprop][sampled_indices], statdata.quantiles[1:]
+                )
 
                 x = util.get_parameter_values("FeH", Grid, selectedmodels, noofind)
-                feh_interval = np.quantile(x[nonzeroprop][sampled_indices], qs[1:])
+                feh_interval = np.quantile(
+                    x[nonzeroprop][sampled_indices], statdata.quantiles[1:]
+                )
 
                 x = util.get_parameter_values("Teff", Grid, selectedmodels, noofind)
-                Teffout = np.quantile(x[nonzeroprop][sampled_indices], qs)
+                Teffout = np.quantile(
+                    x[nonzeroprop][sampled_indices], statdata.quantiles
+                )
 
                 x = util.get_parameter_values("logg", Grid, selectedmodels, noofind)
-                loggout = np.quantile(x[nonzeroprop][sampled_indices], qs)
+                loggout = np.quantile(
+                    x[nonzeroprop][sampled_indices], statdata.quantiles
+                )
 
                 try:
                     fig = plot_kiel.kiel(
@@ -556,18 +563,22 @@ def compute_posterior(
         # Find quantiles of massini/age and FeH to determine what tracks to plot
         library_param = "massini" if "tracks" in gridtype.lower() else "age"
         x = util.get_parameter_values(library_param, Grid, selectedmodels, noofind)
-        lp_interval = np.quantile(x[nonzeroprop][sampled_indices], qs[1:])
+        lp_interval = np.quantile(
+            x[nonzeroprop][sampled_indices], statdata.quantiles[1:]
+        )
 
         # Use correct metallicity (only important for alpha enhancement)
         metalname = "MeH" if "MeH" in fitparams else "FeH"
         x = util.get_parameter_values(metalname, Grid, selectedmodels, noofind)
-        feh_interval = np.quantile(x[nonzeroprop][sampled_indices], qs[1:])
+        feh_interval = np.quantile(
+            x[nonzeroprop][sampled_indices], statdata.quantiles[1:]
+        )
 
         x = util.get_parameter_values("Teff", Grid, selectedmodels, noofind)
-        Teffout = np.quantile(x[nonzeroprop][sampled_indices], qs)
+        Teffout = np.quantile(x[nonzeroprop][sampled_indices], statdata.quantiles)
 
         x = util.get_parameter_values("logg", Grid, selectedmodels, noofind)
-        loggout = np.quantile(x[nonzeroprop][sampled_indices], qs)
+        loggout = np.quantile(x[nonzeroprop][sampled_indices], statdata.quantiles)
 
         try:
             fig = plot_kiel.kiel(
diff --git a/basta/stats.py b/basta/stats.py
index ee51085..9b80fc9 100644
--- a/basta/stats.py
+++ b/basta/stats.py
@@ -1,6 +1,7 @@
 """
 Key statistics functions
 """
+import os
 import copy
 import math
 import collections
@@ -12,7 +13,8 @@
 
 from basta import freq_fit, glitch
 from basta import utils_seismic as su
-from basta.constants import freqtypes
+from basta.constants import sydsun as sydc
+from basta.constants import freqtypes, statdata
 
 # Define named tuple used in selectedmodels
 Trackstats = collections.namedtuple("Trackstats", "index logPDF chi2")
@@ -64,7 +66,7 @@ def chi2_astero(
     ind,
     fitfreqs,
     warnings=True,
-    shapewarn=False,
+    shapewarn=0,
     debug=False,
     verbose=False,
 ):
@@ -99,10 +101,10 @@ def chi2_astero(
         Contains all user inputted frequency fitting options.
     warnings : bool
         If True, print something when it fails.
-    shapewarn : bool
-        If a mismatch in array dimensions of the fitted parameters is
-        encountered, this is set to True in order to warn the user at the end
-        of the run.
+    shapewarn : int
+        If a mismatch in array dimensions of the fitted parameters, or range
+        of frequencies is encountered, this is set to corresponding integer
+        in order to warn the user at the end of the run.
     debug : str
         Flag for print control.
     verbose : str
@@ -216,30 +218,31 @@ def chi2_astero(
 
         # Interpolate model ratios to observed frequencies
         if fitfreqs["interp_ratios"]:
-            # Get extended frequency modes to provide interpolation range of ratios
-            broad_key, broad_osc = su.extend_modjoin(joinkeys, join, modkey, mod)
-            if broad_key is None:
-                shapewarn = True
-                chi2rut = np.inf
-                return chi2rut, warnings, shapewarn
-            # Get model ratios
+            # Get all available model ratios
             broadratio = freq_fit.compute_ratioseqs(
-                broad_key,
-                broad_osc,
+                modkey,
+                mod,
                 ratiotype,
                 threepoint=fitfreqs["threepoint"],
             )
             modratio = copy.deepcopy(obsfreqdata[ratiotype]["data"])
+
             # Seperate and interpolate within the separate r01, r10 and r02 sequences
-            for rtype in set(obsfreqdata[ratiotype]["data"][2, :]):
-                obsmask = obsfreqdata[ratiotype]["data"][2, :] == rtype
+            for rtype in set(modratio[2, :]):
+                obsmask = modratio[2, :] == rtype
                 modmask = broadratio[2, :] == rtype
+                # Check we have the range to interpolate
+                if (
+                    modratio[1, obsmask][0] < broadratio[1, modmask][0]
+                    or modratio[1, obsmask][-1] > broadratio[1, modmask][-1]
+                ):
+                    chi2rut = np.inf
+                    shapewarn = 2
+                    return chi2rut, warnings, shapewarn
                 intfunc = interp1d(
                     broadratio[1, modmask], broadratio[0, modmask], kind="linear"
                 )
-                modratio[0, obsmask] = intfunc(
-                    obsfreqdata[ratiotype]["data"][1, obsmask]
-                )
+                modratio[0, obsmask] = intfunc(modratio[1, obsmask])
 
         else:
             modratio = freq_fit.compute_ratioseqs(
@@ -253,7 +256,7 @@ def chi2_astero(
         if x.shape[0] == covinv.shape[0]:
             chi2rut += (x.T.dot(covinv).dot(x)) / w
         else:
-            shapewarn = True
+            shapewarn = 1
             chi2rut = np.inf
 
     # Add contribution from glitches
@@ -286,14 +289,13 @@ def chi2_astero(
             if x.shape[0] == covinv.shape[0]:
                 chi2rut += (x.T.dot(covinv).dot(x)) / w
             else:
-                shapewarn = True
-                chirut = np.inf
+                shapewarn = 1
+                chi2rut = np.inf
         else:
             chi2rut = np.inf
             return chi2rut, warnings, shapewarn
 
     if any([x in freqtypes.epsdiff for x in fitfreqs["fittypes"]]):
-
         epsdifftype = list(set(fitfreqs["fittypes"]).intersection(freqtypes.epsdiff))[0]
         obsepsdiff = obsfreqdata[epsdifftype]["data"]
         # Purge model freqs of unused modes
@@ -415,7 +417,7 @@ def chi_for_plot(selectedmodels):
     return maxPDFchi2, minchi2
 
 
-def get_highest_likelihood(Grid, selectedmodels, outparams):
+def get_highest_likelihood(Grid, selectedmodels, inputparams):
     """
     Find highest likelihood model and print info.
 
@@ -425,8 +427,9 @@ def get_highest_likelihood(Grid, selectedmodels, outparams):
         The already loaded grid, containing the tracks/isochrones.
     selectedmodels : dict
         Contains information on all models with a non-zero likelihood.
-    outparams : dict
-        Dict containing the wanted output of the run.
+    inputparams : dict
+        The standard bundle of all fitting information
+
 
     Returns
     -------
@@ -448,18 +451,33 @@ def get_highest_likelihood(Grid, selectedmodels, outparams):
         print("  - Name:", Grid[maxPDF_path + "/name"][maxPDF_ind].decode("utf-8"))
 
     # Print parameters
+    outparams = inputparams["asciiparams"]
+    dnu_scales = inputparams.get("dnu_scales", {})
     for param in outparams:
         if param == "distance":
             continue
-        print(
-            "  - {0:10}: {1:12.6f}".format(
-                param, Grid[maxPDF_path + "/" + param][maxPDF_ind]
-            )
-        )
+        paramval = Grid[os.path.join(maxPDF_path, param)][maxPDF_ind]
+
+        # Handle the scaled asteroseismic parameters
+        if param.startswith("dnu") and param not in ["dnufit", "dnufitMos12"]:
+            dnu_rescal = dnu_scales.get(param, 1.00)
+            scaleval = paramval * inputparams.get("dnusun", sydc.SUNdnu) / dnu_rescal
+        elif param.startswith("numax"):
+            scaleval = paramval * inputparams.get("numsun", sydc.SUNnumax)
+        elif param in ["dnufit", "dnufitMos12"]:
+            scaleval = paramval / dnu_scales.get(param, 1.00)
+        else:
+            scaleval = None
+
+        if scaleval:
+            scaleprt = f"(after rescaling: {scaleval:12.6f})"
+        else:
+            scaleprt = ""
+        print("  - {0:10}: {1:12.6f} {2}".format(param, paramval, scaleprt))
     return maxPDF_path, maxPDF_ind
 
 
-def get_lowest_chi2(Grid, selectedmodels, outparams):
+def get_lowest_chi2(Grid, selectedmodels, inputparams):
     """
     Find model with lowest chi2 value and print info.
 
@@ -469,8 +487,8 @@ def get_lowest_chi2(Grid, selectedmodels, outparams):
         The already loaded grid, containing the tracks/isochrones.
     selectedmodels : dict
         Contains information on all models with a non-zero likelihood.
-    outparams : dict
-        Dict containing the wanted output of the run.
+    inputparams : dict
+        The standard bundle of all fitting information
 
     Returns
     -------
@@ -489,14 +507,30 @@ def get_lowest_chi2(Grid, selectedmodels, outparams):
         print("  - Name:", Grid[minchi2_path + "/name"][minchi2_ind].decode("utf-8"))
 
     # Print parameters
+    outparams = inputparams["asciiparams"]
+    dnu_scales = inputparams.get("dnu_scales", {})
     for param in outparams:
         if param == "distance":
             continue
-        print(
-            "  - {0:10}: {1:12.6f}".format(
-                param, Grid[minchi2_path + "/" + param][minchi2_ind]
-            )
-        )
+        paramval = Grid[os.path.join(minchi2_path, param)][minchi2_ind]
+
+        # Handle the scaled asteroseismic parameters
+        if param.startswith("dnu") and param not in ["dnufit", "dnufitMos12"]:
+            dnu_rescal = dnu_scales.get(param, 1.00)
+            scaleval = paramval * inputparams.get("dnusun", sydc.SUNdnu) / dnu_rescal
+        elif param.startswith("numax"):
+            scaleval = paramval * inputparams.get("numsun", sydc.SUNnumax)
+        elif param in ["dnufit", "dnufitMos12"]:
+            scaleval = paramval / dnu_scales.get(param, 1.00)
+        else:
+            scaleval = None
+
+        if scaleval:
+            scaleprt = f"(after rescaling: {scaleval:12.6f})"
+        else:
+            scaleprt = ""
+        print("  - {0:10}: {1:12.6f} {2}".format(param, paramval, scaleprt))
+
     return minchi2_path, minchi2_ind
 
 
@@ -616,15 +650,14 @@ def calc_key_stats(x, centroid, uncert, weights=None):
         84'th percentile if quantile selected, None for standard
         deviation.
     """
-    # Definition of Bayesian 16, 50, and 84 percentiles
-    qs = [0.5, 0.158655, 0.841345]
+
     xp = None
 
     # Handling af all different combinations of input
     if uncert == "quantiles" and not type(weights) == type(None):
-        xcen, xm, xp = quantile_1D(x, weights, qs)
+        xcen, xm, xp = quantile_1D(x, weights, statdata.quantiles)
     elif uncert == "quantiles":
-        xcen, xm, xp = np.quantile(x, qs)
+        xcen, xm, xp = np.quantile(x, statdata.quantiles)
     else:
         xm = np.std(x)
     if centroid == "mean" and not type(weights) == type(None):
diff --git a/basta/utils_general.py b/basta/utils_general.py
index e6da6a3..b8bf7b1 100644
--- a/basta/utils_general.py
+++ b/basta/utils_general.py
@@ -207,8 +207,8 @@ def compare_output_to_input(
         if "distance" in hout_dist:
             idx = np.nonzero([x == "distance_joint" for x in hout_dist])[0][0]
             priordistqs = inputparams["distanceparams"]["priordistance"]
-            priorerrm = priordistqs[1] - priordistqs[0]
-            priorerrp = priordistqs[2] - priordistqs[1]
+            priorerrm = priordistqs[0] - priordistqs[1]
+            priorerrp = priordistqs[2] - priordistqs[0]
             if uncert == "quantiles":
                 outerr = (out_dist[idx + 1] + out_dist[idx + 2]) / 2
             else:
diff --git a/basta/utils_seismic.py b/basta/utils_seismic.py
index 0ffd5f9..74d70ea 100644
--- a/basta/utils_seismic.py
+++ b/basta/utils_seismic.py
@@ -46,19 +46,27 @@ def solar_scaling(Grid, inputparams, diffusion=None):
     """
     print("\nTransforming solar-based asteroseismic quantities:", flush=True)
 
-    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    # BLOCK 1: Conversion into solar units
-    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    print("* Converting to solar units ...")
-
     # Check for solar values, if not set then use default
     dnusun = inputparams.get("dnusun", sydc.SUNdnu)
     numsun = inputparams.get("numsun", sydc.SUNnumax)
 
-    # ------------------------------
-    # TASK 1.1: Conversion of parameters
-    # ------------------------------
+    # Obtain parameter lists
     fitparams = inputparams.get("fitparams")
+    fitfreqs = inputparams.get("fitfreqs", {})
+    limits = inputparams.get("limits")
+
+    # If fitting frequencies, make sure to keep a copy of the original deltaNu
+    if fitfreqs["active"]:
+        fitfreqs["dnu_obs"] = fitfreqs["dnufit"]
+
+    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    # BLOCK 1: Conversion into solar units
+    # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    print("* Converting to solar units if needed...")
+
+    # ----------------------------------
+    # TASK 1.1: Conversion of parameters
+    # ----------------------------------
     fitpar_convert = [
         par for par in fitparams if (par.startswith("dnu") or par.startswith("numax"))
     ]
@@ -87,7 +95,6 @@ def solar_scaling(Grid, inputparams, diffusion=None):
     # TASK 1.2: Conversion of limits
     # ------------------------------
     # Duplicates the approach above)
-    limits = inputparams.get("limits")
     limits_convert = [
         par for par in limits if (par.startswith("dnu") or par.startswith("numax"))
     ]
@@ -159,37 +166,58 @@ def solar_scaling(Grid, inputparams, diffusion=None):
     # ------------------------------
     dnu_scales = {}
     for dnu in sunmoddnu:
-        if dnu in fitparams:
-            if dnu in ["dnufit", "dnufitMos12"]:
-                # Scaling factor is DNU_SUN_GRID / DNU_SUN_OBSERVED
+        if (dnu in fitparams) or (dnu in fitfreqs):
+            if dnu in fitparams:
+                if dnu in ["dnufit", "dnufitMos12"]:
+                    # Scaling factor is DNU_SUN_GRID / DNU_SUN_OBSERVED
+                    dnu_rescal = sunmoddnu[dnu] / dnusun
+                    print(
+                        "  - {0} scaled by {1:.4f} from {2:.2f} to {3:.2f} microHz".format(
+                            dnu,
+                            dnu_rescal,
+                            fitparams[dnu][0],
+                            fitparams[dnu][0] * dnu_rescal,
+                        ),
+                        "(grid Sun: {0:.2f} microHz, real Sun: {1:.2f} microHz)".format(
+                            sunmoddnu[dnu], dnusun
+                        ),
+                    )
+                else:
+                    # Using the scaling relations on the solar model in the grid generally
+                    # yields DNU_SUN_SCALING_GRID != DNU_SUN_SCALING_OBS . The exact value
+                    # of the solar model DNU from scaling relations is stored in the grid
+                    # in solar units (a number close to 1, but not exactly 1). This number
+                    # is used as the scaling factor of solar-unit input dnu's
+                    dnu_rescal = sunmoddnu[dnu]
+                    print(
+                        "  - {0} scaled by a factor {1:.8f} according to the".format(
+                            dnu, dnu_rescal
+                        ),
+                        "grid-solar-model value from scaling relations",
+                    )
+                fitparams[dnu] = [(dnu_rescal) * p for p in fitparams[dnu]]
+            else:
+                # If in frequency fitting, it is always dnufit (in microHz)
+                # --> Scaling factor is DNU_SUN_GRID / DNU_SUN_OBSERVED
                 dnu_rescal = sunmoddnu[dnu] / dnusun
                 print(
                     "  - {0} scaled by {1:.4f} from {2:.2f} to {3:.2f} microHz".format(
                         dnu,
                         dnu_rescal,
-                        fitparams[dnu][0],
-                        fitparams[dnu][0] * dnu_rescal,
+                        fitfreqs[dnu],
+                        fitfreqs[dnu] * dnu_rescal,
                     ),
                     "(grid Sun: {0:.2f} microHz, real Sun: {1:.2f} microHz)".format(
                         sunmoddnu[dnu], dnusun
                     ),
                 )
-            else:
-                # Using the scaling relations on the solar model in the grid generally
-                # yields DNU_SUN_SCALING_GRID != DNU_SUN_SCALING_OBS . The exact value
-                # of the solar model DNU from scaling relations is stored in the grid
-                # in solar units (a number close to 1, but not exactly 1). This number
-                # is used as the scaling factor of solar-unit input dnu's
-                dnu_rescal = sunmoddnu[dnu]
-                print(
-                    "  - {0} scaled by a factor {1:.8f} according to the".format(
-                        dnu, dnu_rescal
-                    ),
-                    "grid-solar-model value from scaling relations",
-                )
+                fitfreqs[dnu] = dnu_rescal * fitfreqs[dnu]
+                if fitfreqs[dnu + "_err"]:
+                    fitfreqs[dnu + "_err"] = dnu_rescal * fitfreqs[dnu + "_err"]
+
             print("    (Note: Will be scaled back before outputting results!)")
-            fitparams[dnu] = [(dnu_rescal) * p for p in fitparams[dnu]]
             dnu_scales[dnu] = dnu_rescal
+
     inputparams["dnu_scales"] = dnu_scales
 
     print("Done!")
@@ -532,9 +560,7 @@ def compute_cov_from_mc(nr, osckey, osc, fittype, args, nrealisations=10000):
         print(f"Frobenius norm {fnorm} > 1e-6")
         print("Warning: Covariance failed to converge")
 
-    n_covinv = np.linalg.pinv(n_cov, rcond=1e-8)
-
-    return n_cov, n_covinv
+    return n_cov
 
 
 def extend_modjoin(joinkey, join, modkey, mod):
diff --git a/basta/utils_xml.py b/basta/utils_xml.py
index 0f5582e..134f263 100644
--- a/basta/utils_xml.py
+++ b/basta/utils_xml.py
@@ -75,6 +75,7 @@ def create_xmltag(
     star = SubElement(main, "star", {"starid": str(starid)})
 
     # Special treatment of dnu and numax for frequency fitting
+    # ++ Make sure to only add them once
     nuset = {"dnu": False, "numax": False}
 
     # Loop over fitting parameters
@@ -104,25 +105,9 @@ def create_xmltag(
         if not np.isclose(paramerr, missingval):
             SubElement(star, param, {"value": str(paramval)})
 
-    # Handle interpolation parameters
-    for param in intpollim:
-        out = {}
-        gparam = "dnu" if "dnu" in param else param
-        nucheck = not nuset[gparam] if gparam in nuset else False
-        if "abstol" in intpollim[param] or nucheck:
-            paramval = _get_param(paramvals, params, gparam)
-            if not np.isclose(paramval, missingval):
-                out["value"] = str(paramval)
-            if nucheck:
-                nuset[gparam] = True
-        if "sigmacut" in intpollim[param]:
-            paramerr = _get_param(paramvals, params, gparam + "_err")
-            if not np.isclose(paramerr, missingval):
-                out["error"] = str(paramerr)
-        SubElement(star, gparam, out)
-
-    # Handle the nottrustedfile object
+    # Handle the special things for frequency fitting
     if freqparams and any(x in fitparams for x in freqtypes.alltypes):
+        # The nottrustedfile object
         if "nottrustedfile" in freqparams:
             if isinstance(freqparams["nottrustedfile"], dict):
                 if starid in freqparams["nottrustedfile"].keys():
@@ -137,12 +122,40 @@ def create_xmltag(
                 )
             else:
                 raise ValueError("Nottrustedfile is neither a dict or a str")
+
+        # Always add dnu and numax for frequency fitting
         if not nuset["dnu"]:
             dnu = _get_param(paramvals, params, "dnu")
-            SubElement(star, "dnu", {"value": str(dnu)})
+            try:
+                dnu_err = _get_param(paramvals, params, "dnu_err")
+                SubElement(star, "dnu", {"value": str(dnu), "error": str(dnu_err)})
+            except IndexError:
+                SubElement(star, "dnu", {"value": str(dnu)})
+            nuset["dnu"] = True
+
         if not nuset["numax"]:
             numax = _get_param(paramvals, params, "numax")
             SubElement(star, "numax", {"value": str(numax)})
+            nuset["numax"] = True
+
+    # Handle interpolation parameters (special treatment of dnu; not add if added)
+    for param in intpollim:
+        print(param)
+        out = {}
+        gparam = "dnu" if "dnu" in param else param
+        nucheck = not nuset[gparam] if gparam in nuset else False
+        if "abstol" in intpollim[param] or nucheck:
+            paramval = _get_param(paramvals, params, gparam)
+            if not np.isclose(paramval, missingval):
+                out["value"] = str(paramval)
+            if nucheck:
+                nuset[gparam] = True
+        if "sigmacut" in intpollim[param]:
+            paramerr = _get_param(paramvals, params, gparam + "_err")
+            if not np.isclose(paramerr, missingval):
+                out["error"] = str(paramerr)
+        if not nuset[gparam]:
+            SubElement(star, gparam, out)
 
     return main
 
diff --git a/basta/xml_create.py b/basta/xml_create.py
index 5bb82e9..23aadb2 100644
--- a/basta/xml_create.py
+++ b/basta/xml_create.py
@@ -309,7 +309,6 @@ def generate_xml(
         or ("distance" in outparams)
         or ("distance" in cornerplots)
     ):
-
         # Convert to tuple if the user specified only one filter as a sting
         if isinstance(filters, str):
             filters = (filters,)
diff --git a/basta/xml_run.py b/basta/xml_run.py
index abd9661..9c7d91d 100644
--- a/basta/xml_run.py
+++ b/basta/xml_run.py
@@ -520,6 +520,9 @@ def run_xml(
         fitfreqs["nsorting"] = strtobool(
             _find_get(root, "default/freqparams/nsorting", "value", "True")
         )
+        fitfreqs["dnuprior"] = strtobool(
+            _find_get(root, "default/freqparams/dnuprior", "value", "True")
+        )
 
         # Read seismic weight quantities
         dof = _find_get(root, "default/freqparams/dof", "value", None)
@@ -722,6 +725,10 @@ def run_xml(
                 # dnufit for prior, numax for scaling
                 fitfreqs["dnufit"] = float(_find_get(star, "dnu", "value"))
                 fitfreqs["numax"] = float(_find_get(star, "numax", "value"))
+                try:
+                    fitfreqs["dnufit_err"] = float(_find_get(star, "dnu", "error"))
+                except ValueError:
+                    pass
 
             # Add to inputparams
             inputparams["fitfreqs"] = fitfreqs
diff --git a/docs/conf.py b/docs/conf.py
index 6ae1be7..2cc750e 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -31,10 +31,10 @@
 author = "The BASTA Team"
 
 # The short X.Y version.
-version = "1.2.0"
+version = "1.2.1"
 
 # The full version, including alpha/beta/rc tags.
-release = "1.2.0"
+release = "1.2.1"
 
 
 # -- General configuration ------------------------------------------------
@@ -71,7 +71,7 @@
 #
 # This is also used if you do content translation via gettext catalogs.
 # Usually you set "language" from the command line for these cases.
-language = None
+language = "en"
 
 # List of patterns, relative to source directory, that match files and
 # directories to ignore when looking for source files.
diff --git a/docs/examples.rst b/docs/examples.rst
index 763db99..ba3670f 100644
--- a/docs/examples.rst
+++ b/docs/examples.rst
@@ -3,8 +3,6 @@
 Examples of fits
 ================
 
-*We are currently having issues with the examples pages (images are not showing in Firefox and Chrome). We are working on a release to fix it. Until then, please try with Safari. We are sorry for the inconvenience!*
-
 The following sections give practical examples of commonly used fits with BASTA, and explain the options in the code following the template given in :py:func:`define_input` in ``${BASTADIR}/examples/create_inputfile.py`` (in the following this will be denoted as :py:meth:`create_inputfile.define_input`).
 
 The fitting examples are located in ``${BASTADIR}/examples/xmlinput/`` and a reference set of output can be found in ``${BASTADIR}/examples/reference/``, which includes all fits described in the following sections. These can be compared to the output you obtain when running the code as a further sanity check. Overview of fitting examples:
diff --git a/docs/examples_freqs.rst b/docs/examples_freqs.rst
index 9147694..c0b3fdd 100644
--- a/docs/examples_freqs.rst
+++ b/docs/examples_freqs.rst
@@ -1,11 +1,8 @@
 .. _example_freqs:
 
-Individual frequencies, ratios, glitches
+Methods using individual frequencies
 ========================================
 
-*We are currently having issues with the examples pages (images are not showing in Firefox and Chrome). We are working on a release to fix it. Until then, please try with Safari. We are sorry for the inconvenience!*
-
-
 Using grids that include theoretically computed oscillation frequencies (see :ref:`grids`) BASTA can fit these
 individual frequencies with a surface correction, as well as combination of frequencies. In the following we show
 examples of the blocks that must be modified in :py:meth:`create_inputfile.define_input` to produce these types of fits.
@@ -93,12 +90,12 @@ The fit should take less than a minute and the output is stored in ``${BASTADIR}
 corner plot and Kiel diagrams, the code produces output of the fit to the individual frequencies in form of echelle
 diagrams for both corrected and uncorrected frequencies:
 
-.. figure:: ../examples/reference/freqs/16CygA_pairechelle_uncorrected.pdf
+.. figure:: figures/freqs/16CygA_pairechelle_uncorrected.png
    :alt: Echelle diagram showing the uncorrected frequencies of the best fit model to 16 Cyg A in the grid.
 
    Echelle diagram showing the uncorrected frequencies of the best fit model to 16 Cyg A in the grid.
 
-.. figure:: ../examples/reference/freqs/16CygA_pairechelle.pdf
+.. figure:: figures/freqs/16CygA_pairechelle.png
    :alt: Echelle diagram after the BG14 frequency correction to the best fit model to 16 Cyg A in the grid.
 
    Echelle diagram after the BG14 frequency correction to the best fit model to 16 Cyg A in the grid.
@@ -124,18 +121,59 @@ one simply adds the following ``fitparam`` (for the case of :math:`r_{012}` as a
 
 The variable ``freqplots`` can also be set to ``True``, which will produce plots of the ratios and corresponding echelle
 diagrams even though individual frequencies are not fitted. We provide an example to run this fit in
-``${BASTADIR}/examples/xmlinput/create_inputfiles_ratios.py`` which produces the file ``input_ratios.xml``. Running
+``${BASTADIR}/examples/xmlinput/create_inputfile_ratios.py`` which produces the file ``input_ratios.xml``. Running
 this file stores the results of the fit in ``${BASTADIR}/examples/output/ratios/``, and the resulting ratios should look
 as follows:
 
-.. figure:: ../examples/reference/ratios/16CygA_ratios.pdf
+.. figure:: figures/ratios/16CygA_ratios_r012.png
    :alt: Frequency ratios of the best fit model to 16 Cyg A in the grid.
 
    Frequency ratios of the best fit model to 16 Cyg A in the grid.
 
 BASTA uses by default the five-point small frequency separation formulation for computing the ratios, which is the
-recommended option. If the user should want to use the three-point formulation instead, this can be done by adding
-``"threepoint": True`` in the ``define_fit["freqparams"]`` dictionary.
+recommended option. Additionally, interpolation of the model ratios to the observed frequencies are applied in the fit.
+Finally, the correlations between the ratios are taken into account by using the full covariance matrix. Any of these
+settings can of cource be changed should the user wish to do so.
+
+
+Epsilon differences
+-------------------
+
+Similar to frequency ratios, BASTA can also fit the surface-independent frequency phase differences, commonly
+referred to as epsilon differences (Winther et. al, in preparation). The individual set of differences
+(:math:`\delta\epsilon_{01}, \delta\epsilon_{02}`) as well as the combined set can be fitted by adding the
+correpsonding keyword to ``fitparams`` (here for the case :math:`\delta\epsilon_{012}`):
+
+.. code-block:: python
+
+    # ==================================================================================
+    # BLOCK 2: Fitting control
+    # ==================================================================================
+    define_fit["fitparams"] = ("Teff", "FeH", "e012")
+
+    # ==================================================================================
+    # BLOCK 4: Plotting control
+    # ==================================================================================
+    define_plots["freqplots"] = "epsdiff"
+
+Adding ``epsdiff`` to ``freqplots`` produces the corresponding figure, which can also generally be produced when
+individual frequencies are available. An example of how to run this fit is provided in
+``${BASTADIR}/examples/xmlinput/create_inputfile_epsilondifference.py`` which produces the file ``input_epsilondifference.xml``.
+Running this file stores the results of the fit in ``${BASTADIR}/examples/output/epsilon/``, and the resulting
+epsilon differences should look as follows:
+
+.. figure:: figures/epsilon/16CygA_epsdiff_e012.png
+    :alt: Epsilon differences of the best fit model to 16 Cyg A in the grid.
+
+    Epsilon differences of the best fit model to 16 Cyg A in the grid.
+
+Note that since the determination of epsilon differences relies on interpolating the :math:`\ell=0` epsilons to the frequency locations
+of the :math:`\ell=1,2` modes, one would extrapolate the :math:`\ell=0` epsilons if the frequency locations of the
+:math:`\ell=1,2` goes outside the interval of the frequency locations of the :math:`\ell=0` modes. These are therefore
+excluded, and thus the number of epsilon differences may not be equal to the number of :math:`\ell=1,2` frequencies.
+
+As noted above for the ratios, correlations/covariances are taken into account in the fit.
+
 
 Frequency glitches
 ------------------
@@ -182,9 +220,9 @@ Since the ``.glh`` file is located in the same folder as the individual frequenc
     }
 
 You can find the corresponding python script to produce the input file for this fit in
-``${BASTADIR}/examples/xmlinput/create_inputfiles_glitches.py``. The output should look as follows:
+``${BASTADIR}/examples/xmlinput/create_inputfile_glitches.py``. The output should look as follows:
 
-.. figure:: ../examples/reference/glitches/16CygA_corner.pdf
+.. figure:: figures/glitches/16CygA_corner.png
    :alt: Corner plot of the 16 Cyg A fit using glitches.
 
    Corner plot of the 16 Cyg A fit using glitches.
@@ -264,7 +302,7 @@ simply be run as
 
 The resulting duplicated echelle diagram should look as like the following.
 
-.. figure:: ../examples/reference/subgiant/Valid_245_dupechelle.pdf
+.. figure:: figures/subgiant/Valid_245_dupechelle.png
    :alt: Echelle diagram after the BG14 frequency correction to the best fit model to Validation star 245.
 
    Echelle diagram after the BG14 frequency correction to the best fit model to Validation star 245.
diff --git a/docs/examples_global.rst b/docs/examples_global.rst
index cd461a0..40aea1a 100644
--- a/docs/examples_global.rst
+++ b/docs/examples_global.rst
@@ -3,8 +3,6 @@
 Spectroscopy and global asteroseismic parameters
 ================================================
 
-*We are currently having issues with the examples pages (images are not showing in Firefox and Chrome). We are working on a release to fix it. Until then, please try with Safari. We are sorry for the inconvenience!*
-
 This example describes the fitting of what we call global parameters of the grid, with a specific application to
 asteroseismology for illustration purposes. As a general recommendation, the user should explore the
 :py:meth:`constants.parameters` (in ``${BASTADIR}/basta/constants.py``) function for a complete list of available fitting parameters.
@@ -86,9 +84,9 @@ You should see the following output printed in your terminal:
 
     A total of 1 star(s) will be fitted with {Teff, FeH, dnufit, numax} to the grid 'BASTADIR/grids/Garstec_16CygA.hdf5'.
 
-    This will output {Teff, FeH, radPhot, massfin, age} to a results file.
+    This will output {Teff, FeH, dnufit, numax, radPhot, massfin, age} to a results file.
 
-    Corner plots include {Teff, FeH, radPhot, massfin, age} with observational bands on {Teff, FeH, dnufit, numax}.
+    Corner plots include {Teff, FeH, dnufit, numax, radPhot, massfin, age} with observational bands on {Teff, FeH, dnufit, numax}.
     Kiel diagrams will be made with observational bands on {Teff, FeH, dnufit, numax}.
 
     A restricted flat prior will be applied to: Teff, FeH.
@@ -106,7 +104,7 @@ Once the file is created, run BASTA as explained to perform the fit:
 The output of the fit can be found in ``${BASTADIR}/examples/output/global/``. It includes a Kiel diagram that should
 look like the following:
 
-.. figure:: ../examples/reference/global/16CygA_kiel.pdf
+.. figure:: figures/global/16CygA_kiel.png
    :alt: Kiel diagram plot of the 16 Cyg A fit using global asteroseismic quantities.
 
    Kiel diagram of the 16 Cyg A fit using global asteroseismic quantities.
@@ -117,7 +115,7 @@ and 84 percentiles mass and metallicity output of the solution, and are **not**
 
 Finally, a corner plot of the parameters included in ``cornerplots`` is also part of the output:
 
-.. figure:: ../examples/reference/global/16CygA_corner.pdf
+.. figure:: figures/global/16CygA_corner.png
    :alt: Corner plot of the 16 Cyg A fit using global asteroseismic quantities.
 
    Corner plot of the 16 Cyg A fit using global asteroseismic quantities.
diff --git a/docs/examples_interpolation.rst b/docs/examples_interpolation.rst
index 4a2a52f..d965fc6 100644
--- a/docs/examples_interpolation.rst
+++ b/docs/examples_interpolation.rst
@@ -3,8 +3,6 @@
 Interpolation of tracks and isochrones
 =======================================
 
-*We are currently having issues with the examples pages (images are not showing in Firefox and Chrome). We are working on a release to fix it. Until then, please try with Safari. We are sorry for the inconvenience!*
-
 If the resolutions of a particular grid is considered insufficient for a particular fit, BASTA includes the option of
 performing interpolation across and/or along tracks and isochrones. What the procedure does is to produce a new
 interpolated grid of models, store it in ``hdf5`` format, and perform a fit to the input given in ``fitparams``. Since
@@ -74,12 +72,12 @@ Finally, we define the settings for the interpolation `across/between` tracks:
 Since the input grid we are using has been constructed using Sobol sampling, we define a ``scale`` parameter of 5 in
 ``["resolution"]``, which increases the number of tracks with respect to the original by a factor of 6. The following figures shows the new sampling, and compares the desired resolution to the `current` distributions of parameter values in the grid:
 
-.. figure:: ../examples/reference/preview_interp_MS/16CygA_interp_preview_across_resolution.pdf
+.. figure:: figures/preview_interp_MS/16CygA_interp_preview_across_resolution.png
    :alt: Distribution in mass and metallicity of the current grid and desired interpolated grid.
 
    Distribution in mass and metallicity of the current grid (base) and desired interpolated grid.
 
-.. figure:: ../examples/reference/preview_interp_MS/16CygA_interp_preview_along_resolution.pdf
+.. figure:: figures/preview_interp_MS/16CygA_interp_preview_along_resolution.png
    :alt: Distribution in individual frequency resolution of the current grid.
 
    Distribution in individual frequency resolution of the current grid.
@@ -152,22 +150,22 @@ After the interpolation and fit are performed the results are stored in ``${BAST
 including the new interpolated grid. The following figures compare the Kiel diagrams of the grids with and without
 interpolation, as well as the corner plots.
 
-.. figure:: ../examples/reference/freqs/16CygA_kiel.pdf
+.. figure:: figures/freqs/16CygA_kiel.png
    :alt: Kiel diagram of the 16 Cyg A fit using the original grid.
 
    Kiel diagram of the 16 Cyg A fit using the original grid.
 
-.. figure:: ../examples/reference/interp_MS/16CygA_kiel.pdf
+.. figure:: figures/interp_MS/16CygA_kiel.png
    :alt: Kiel diagram of the 16 Cyg A fit using the interpolated grid.
 
    Kiel diagram of the 16 Cyg A fit using the interpolated grid.
 
-.. figure:: ../examples/reference/freqs/16CygA_corner.pdf
+.. figure:: figures/freqs/16CygA_corner.png
    :alt: Corner plot of the 16 Cyg A fit using the original grid.
 
    Corner plot of the 16 Cyg A fit using the original grid.
 
-.. figure:: ../examples/reference/interp_MS/16CygA_corner.pdf
+.. figure:: figures/interp_MS/16CygA_corner.png
    :alt: Corner plot of the 16 Cyg A fit using the interpolated grid.
 
    Corner plot of the 16 Cyg A fit using the interpolated grid.
diff --git a/docs/examples_isochrones.rst b/docs/examples_isochrones.rst
index 91df6f4..6d7c30f 100644
--- a/docs/examples_isochrones.rst
+++ b/docs/examples_isochrones.rst
@@ -3,8 +3,6 @@
 Isochrones, distances, and evolutionary phase
 =============================================
 
-*We are currently having issues with the examples pages (images are not showing in Firefox and Chrome). We are working on a release to fix it. Until then, please try with Safari. We are sorry for the inconvenience!*
-
 BASTA is shipped with a complete set of isochrones from the
 `BaSTI library <http://basti-iac.oa-abruzzo.inaf.it/index.html>`_. The input physics, science cases, and mixtures are
 described in the `Solar scaled paper <https://ui.adsabs.harvard.edu/abs/2018ApJ...856..125H/abstract>`_ and the
@@ -32,8 +30,8 @@ are the following:
         "DEC",
         "Teff",
         "Teff_err",
-        "MeH",
-        "MeH_err",
+        "FeH",
+        "FeH_err",
         "dnu",
         "dnu_err",
         "numax",
@@ -73,7 +71,7 @@ The first fit in this example will include spectroscopy, global asteroseismic pa
     # ==================================================================================
     # BLOCK 2: Fitting control
     # ==================================================================================
-    define_fit["fitparams"] = ("Teff", "MeH", "dnuSer", "numax", "parallax")
+    define_fit["fitparams"] = ("Teff", "FeH", "dnuSer", "numax", "parallax")
 
 Note that we fit the large frequency separation using the
 `Serenelli et al. 2017 correction <https://ui.adsabs.harvard.edu/abs/2017ApJS..233...23S/abstract>`_ (called ``dnuSer``)
@@ -111,12 +109,12 @@ isochrones library:
 The file ``${BASTADIR}/examples/xmlinput/create_inputfile_parallax.py`` has been modified accordingly to produce a fit
 to the isochrones. After running the associated ``input_parallax.xml`` the output should look as follows:
 
-.. figure:: ../examples/reference/parallax/11129442_kiel.pdf
+.. figure:: figures/parallax/11129442_kiel.png
    :alt: Kiel diagram of the fit to KIC 11129442 including asteroseismology and parallaxes.
 
    Kiel diagram of the fit to KIC 11129442 including asteroseismology and parallaxes.
 
-.. figure:: ../examples/reference/parallax/11129442_corner.pdf
+.. figure:: figures/parallax/11129442_corner.png
    :alt: Corner plot of the fit to KIC 11129442 including asteroseismology and parallaxes.
 
    Corner plot of the fit to KIC 11129442 including asteroseismology and parallaxes.
@@ -124,7 +122,7 @@ to the isochrones. After running the associated ``input_parallax.xml`` the outpu
 You may have noticed that there is one additional figure to this fit that has not appeared before. This is the corner
 plot of the parameters associated to the distance called ``11129442_distance_corner.pdf`` and looks like this:
 
-.. figure:: ../examples/reference/parallax/11129442_distance_corner.pdf
+.. figure:: figures/parallax/11129442_distance_corner.png
    :alt: Corner plot of distance properties of the fit to KIC 11129442 including asteroseismology and parallaxes.
 
    Corner plot of distance properties of the fit to KIC 11129442 including asteroseismology and parallaxes.
@@ -151,7 +149,7 @@ parallax must not be present in ``fitparams``:
     # ==================================================================================
     # BLOCK 2: Fitting control
     # ==================================================================================
-    define_fit["fitparams"] = ("Teff", "MeH", "dnuSer", "numax")
+    define_fit["fitparams"] = ("Teff", "FeH", "dnuSer", "numax")
 
 At least one magnitude must be specified (and of course its observed value and uncertainty should be included in the
 ``ascii`` file containing the inout data). In this case we will use only the 2MASS *J* magnitude:
@@ -170,11 +168,11 @@ And finally the parameter ``distance`` must be included in the output parameters
     # ==================================================================================
     # BLOCK 3: Output control
     # ==================================================================================
-    define_output["outparams"] = ("Teff", "MeH", "rho", "radPhot", "massfin", "age", "distance")
+    define_output["outparams"] = ("Teff", "FeH", "dnuSer", "numax", "radPhot", "massfin", "age", "distance")
 
 The distance corner plot including only the *J* magnitude looks as follows:
 
-.. figure:: ../examples/reference/distance/11129442_distance_corner.pdf
+.. figure:: figures/distance/11129442_distance_corner.png
    :alt: Distance corner plot for KIC 11129442 estimated using asteroseismology, spectroscopy, and one apparent magnitude.
 
    Distance corner plot for KIC 11129442 estimated using asteroseismology, spectroscopy, and one apparent magnitude.
@@ -208,12 +206,12 @@ procedure. The only difference (except for reading the pre-modified data file) i
 
 .. code-block:: python
 
-    define_fit["fitparams"] = ("Teff", "MeH", "dnuSer", "numax", "phase")
+    define_fit["fitparams"] = ("Teff", "FeH", "dnuSer", "numax", "phase")
 
 The resulting fit is displayed in the following Kiel diagram where it is clear that the red clump phase has been
 selected instead of the RGB
 
-.. figure:: ../examples/reference/phase/11129442_kiel.pdf
+.. figure:: figures/phase/11129442_kiel.png
    :alt: Kiel diagram of the fit to KIC 11129442 forcing the star into the red clump.
 
    Kiel diagram of the fit to KIC 11129442 forcing the star into the red clump.
diff --git a/docs/figures/distance/11129442_corner.png b/docs/figures/distance/11129442_corner.png
new file mode 100644
index 0000000..8898aef
Binary files /dev/null and b/docs/figures/distance/11129442_corner.png differ
diff --git a/docs/figures/distance/11129442_distance_corner.png b/docs/figures/distance/11129442_distance_corner.png
new file mode 100644
index 0000000..fe0d819
Binary files /dev/null and b/docs/figures/distance/11129442_distance_corner.png differ
diff --git a/docs/figures/distance/11129442_kiel.png b/docs/figures/distance/11129442_kiel.png
new file mode 100644
index 0000000..da61176
Binary files /dev/null and b/docs/figures/distance/11129442_kiel.png differ
diff --git a/docs/figures/epsilon/16CygA_corner.png b/docs/figures/epsilon/16CygA_corner.png
new file mode 100644
index 0000000..e732e89
Binary files /dev/null and b/docs/figures/epsilon/16CygA_corner.png differ
diff --git a/docs/figures/epsilon/16CygA_dupechelle.png b/docs/figures/epsilon/16CygA_dupechelle.png
new file mode 100644
index 0000000..b2ab8e3
Binary files /dev/null and b/docs/figures/epsilon/16CygA_dupechelle.png differ
diff --git a/docs/figures/epsilon/16CygA_dupechelle_uncorrected.png b/docs/figures/epsilon/16CygA_dupechelle_uncorrected.png
new file mode 100644
index 0000000..47ff5fd
Binary files /dev/null and b/docs/figures/epsilon/16CygA_dupechelle_uncorrected.png differ
diff --git a/docs/figures/epsilon/16CygA_echelle.png b/docs/figures/epsilon/16CygA_echelle.png
new file mode 100644
index 0000000..b585aad
Binary files /dev/null and b/docs/figures/epsilon/16CygA_echelle.png differ
diff --git a/docs/figures/epsilon/16CygA_echelle_uncorrected.png b/docs/figures/epsilon/16CygA_echelle_uncorrected.png
new file mode 100644
index 0000000..afa72fb
Binary files /dev/null and b/docs/figures/epsilon/16CygA_echelle_uncorrected.png differ
diff --git a/docs/figures/epsilon/16CygA_epsdiff_e012.png b/docs/figures/epsilon/16CygA_epsdiff_e012.png
new file mode 100644
index 0000000..dd238a3
Binary files /dev/null and b/docs/figures/epsilon/16CygA_epsdiff_e012.png differ
diff --git a/docs/figures/epsilon/16CygA_epsdiff_e012_cormap.png b/docs/figures/epsilon/16CygA_epsdiff_e012_cormap.png
new file mode 100644
index 0000000..b3e2744
Binary files /dev/null and b/docs/figures/epsilon/16CygA_epsdiff_e012_cormap.png differ
diff --git a/docs/figures/epsilon/16CygA_kiel.png b/docs/figures/epsilon/16CygA_kiel.png
new file mode 100644
index 0000000..a3b9c2e
Binary files /dev/null and b/docs/figures/epsilon/16CygA_kiel.png differ
diff --git a/docs/figures/epsilon/16CygA_pairechelle.png b/docs/figures/epsilon/16CygA_pairechelle.png
new file mode 100644
index 0000000..b4f3711
Binary files /dev/null and b/docs/figures/epsilon/16CygA_pairechelle.png differ
diff --git a/docs/figures/epsilon/16CygA_pairechelle_uncorrected.png b/docs/figures/epsilon/16CygA_pairechelle_uncorrected.png
new file mode 100644
index 0000000..07448c1
Binary files /dev/null and b/docs/figures/epsilon/16CygA_pairechelle_uncorrected.png differ
diff --git a/docs/figures/epsilon/16CygA_ratios_r012.png b/docs/figures/epsilon/16CygA_ratios_r012.png
new file mode 100644
index 0000000..00108b5
Binary files /dev/null and b/docs/figures/epsilon/16CygA_ratios_r012.png differ
diff --git a/docs/figures/epsilon/16CygA_ratios_r012_cormap.png b/docs/figures/epsilon/16CygA_ratios_r012_cormap.png
new file mode 100644
index 0000000..fca3a8f
Binary files /dev/null and b/docs/figures/epsilon/16CygA_ratios_r012_cormap.png differ
diff --git a/docs/figures/freqs/16CygA_corner.png b/docs/figures/freqs/16CygA_corner.png
new file mode 100644
index 0000000..1495703
Binary files /dev/null and b/docs/figures/freqs/16CygA_corner.png differ
diff --git a/docs/figures/freqs/16CygA_dupechelle.png b/docs/figures/freqs/16CygA_dupechelle.png
new file mode 100644
index 0000000..3f31f53
Binary files /dev/null and b/docs/figures/freqs/16CygA_dupechelle.png differ
diff --git a/docs/figures/freqs/16CygA_dupechelle_uncorrected.png b/docs/figures/freqs/16CygA_dupechelle_uncorrected.png
new file mode 100644
index 0000000..3f71bdc
Binary files /dev/null and b/docs/figures/freqs/16CygA_dupechelle_uncorrected.png differ
diff --git a/docs/figures/freqs/16CygA_echelle.png b/docs/figures/freqs/16CygA_echelle.png
new file mode 100644
index 0000000..b38fdf4
Binary files /dev/null and b/docs/figures/freqs/16CygA_echelle.png differ
diff --git a/docs/figures/freqs/16CygA_echelle_uncorrected.png b/docs/figures/freqs/16CygA_echelle_uncorrected.png
new file mode 100644
index 0000000..5b39019
Binary files /dev/null and b/docs/figures/freqs/16CygA_echelle_uncorrected.png differ
diff --git a/docs/figures/freqs/16CygA_kiel.png b/docs/figures/freqs/16CygA_kiel.png
new file mode 100644
index 0000000..6d0e6a8
Binary files /dev/null and b/docs/figures/freqs/16CygA_kiel.png differ
diff --git a/docs/figures/freqs/16CygA_pairechelle.png b/docs/figures/freqs/16CygA_pairechelle.png
new file mode 100644
index 0000000..60ae9bf
Binary files /dev/null and b/docs/figures/freqs/16CygA_pairechelle.png differ
diff --git a/docs/figures/freqs/16CygA_pairechelle_uncorrected.png b/docs/figures/freqs/16CygA_pairechelle_uncorrected.png
new file mode 100644
index 0000000..a8e3b80
Binary files /dev/null and b/docs/figures/freqs/16CygA_pairechelle_uncorrected.png differ
diff --git a/docs/figures/glitches/16CygA_corner.png b/docs/figures/glitches/16CygA_corner.png
new file mode 100644
index 0000000..e5ff3cc
Binary files /dev/null and b/docs/figures/glitches/16CygA_corner.png differ
diff --git a/docs/figures/glitches/16CygA_dupechelle.png b/docs/figures/glitches/16CygA_dupechelle.png
new file mode 100644
index 0000000..55b1296
Binary files /dev/null and b/docs/figures/glitches/16CygA_dupechelle.png differ
diff --git a/docs/figures/glitches/16CygA_dupechelle_uncorrected.png b/docs/figures/glitches/16CygA_dupechelle_uncorrected.png
new file mode 100644
index 0000000..5a9be15
Binary files /dev/null and b/docs/figures/glitches/16CygA_dupechelle_uncorrected.png differ
diff --git a/docs/figures/glitches/16CygA_echelle.png b/docs/figures/glitches/16CygA_echelle.png
new file mode 100644
index 0000000..f7599ac
Binary files /dev/null and b/docs/figures/glitches/16CygA_echelle.png differ
diff --git a/docs/figures/glitches/16CygA_echelle_uncorrected.png b/docs/figures/glitches/16CygA_echelle_uncorrected.png
new file mode 100644
index 0000000..56c05a0
Binary files /dev/null and b/docs/figures/glitches/16CygA_echelle_uncorrected.png differ
diff --git a/docs/figures/glitches/16CygA_kiel.png b/docs/figures/glitches/16CygA_kiel.png
new file mode 100644
index 0000000..034d46e
Binary files /dev/null and b/docs/figures/glitches/16CygA_kiel.png differ
diff --git a/docs/figures/glitches/16CygA_pairechelle.png b/docs/figures/glitches/16CygA_pairechelle.png
new file mode 100644
index 0000000..03d5e9b
Binary files /dev/null and b/docs/figures/glitches/16CygA_pairechelle.png differ
diff --git a/docs/figures/glitches/16CygA_pairechelle_uncorrected.png b/docs/figures/glitches/16CygA_pairechelle_uncorrected.png
new file mode 100644
index 0000000..924044a
Binary files /dev/null and b/docs/figures/glitches/16CygA_pairechelle_uncorrected.png differ
diff --git a/docs/figures/global/16CygA_corner.png b/docs/figures/global/16CygA_corner.png
new file mode 100644
index 0000000..0b56952
Binary files /dev/null and b/docs/figures/global/16CygA_corner.png differ
diff --git a/docs/figures/global/16CygA_kiel.png b/docs/figures/global/16CygA_kiel.png
new file mode 100644
index 0000000..5ab53a9
Binary files /dev/null and b/docs/figures/global/16CygA_kiel.png differ
diff --git a/docs/figures/interp_MS/16CygA_corner.png b/docs/figures/interp_MS/16CygA_corner.png
new file mode 100644
index 0000000..0480327
Binary files /dev/null and b/docs/figures/interp_MS/16CygA_corner.png differ
diff --git a/docs/figures/interp_MS/16CygA_dupechelle.png b/docs/figures/interp_MS/16CygA_dupechelle.png
new file mode 100644
index 0000000..7626139
Binary files /dev/null and b/docs/figures/interp_MS/16CygA_dupechelle.png differ
diff --git a/docs/figures/interp_MS/16CygA_dupechelle_uncorrected.png b/docs/figures/interp_MS/16CygA_dupechelle_uncorrected.png
new file mode 100644
index 0000000..0184c65
Binary files /dev/null and b/docs/figures/interp_MS/16CygA_dupechelle_uncorrected.png differ
diff --git a/docs/figures/interp_MS/16CygA_echelle.png b/docs/figures/interp_MS/16CygA_echelle.png
new file mode 100644
index 0000000..66e95b2
Binary files /dev/null and b/docs/figures/interp_MS/16CygA_echelle.png differ
diff --git a/docs/figures/interp_MS/16CygA_echelle_uncorrected.png b/docs/figures/interp_MS/16CygA_echelle_uncorrected.png
new file mode 100644
index 0000000..279c894
Binary files /dev/null and b/docs/figures/interp_MS/16CygA_echelle_uncorrected.png differ
diff --git a/docs/figures/interp_MS/16CygA_kiel.png b/docs/figures/interp_MS/16CygA_kiel.png
new file mode 100644
index 0000000..173ce13
Binary files /dev/null and b/docs/figures/interp_MS/16CygA_kiel.png differ
diff --git a/docs/figures/interp_MS/16CygA_pairechelle.png b/docs/figures/interp_MS/16CygA_pairechelle.png
new file mode 100644
index 0000000..11da85c
Binary files /dev/null and b/docs/figures/interp_MS/16CygA_pairechelle.png differ
diff --git a/docs/figures/interp_MS/16CygA_pairechelle_uncorrected.png b/docs/figures/interp_MS/16CygA_pairechelle_uncorrected.png
new file mode 100644
index 0000000..01e6cb1
Binary files /dev/null and b/docs/figures/interp_MS/16CygA_pairechelle_uncorrected.png differ
diff --git a/docs/figures/interp_MS/intpol_example_16CygA_intpolbase.png b/docs/figures/interp_MS/intpol_example_16CygA_intpolbase.png
new file mode 100644
index 0000000..e8eadcb
Binary files /dev/null and b/docs/figures/interp_MS/intpol_example_16CygA_intpolbase.png differ
diff --git a/docs/figures/json/16CygA_corner.png b/docs/figures/json/16CygA_corner.png
new file mode 100644
index 0000000..f24ac5e
Binary files /dev/null and b/docs/figures/json/16CygA_corner.png differ
diff --git a/docs/figures/json/16CygA_kiel.png b/docs/figures/json/16CygA_kiel.png
new file mode 100644
index 0000000..b270c79
Binary files /dev/null and b/docs/figures/json/16CygA_kiel.png differ
diff --git a/docs/figures/parallax/11129442_corner.png b/docs/figures/parallax/11129442_corner.png
new file mode 100644
index 0000000..ad3295c
Binary files /dev/null and b/docs/figures/parallax/11129442_corner.png differ
diff --git a/docs/figures/parallax/11129442_distance_corner.png b/docs/figures/parallax/11129442_distance_corner.png
new file mode 100644
index 0000000..b99c0ec
Binary files /dev/null and b/docs/figures/parallax/11129442_distance_corner.png differ
diff --git a/docs/figures/parallax/11129442_kiel.png b/docs/figures/parallax/11129442_kiel.png
new file mode 100644
index 0000000..9166538
Binary files /dev/null and b/docs/figures/parallax/11129442_kiel.png differ
diff --git a/docs/figures/phase/11129442_corner.png b/docs/figures/phase/11129442_corner.png
new file mode 100644
index 0000000..268acef
Binary files /dev/null and b/docs/figures/phase/11129442_corner.png differ
diff --git a/docs/figures/phase/11129442_kiel.png b/docs/figures/phase/11129442_kiel.png
new file mode 100644
index 0000000..50c8372
Binary files /dev/null and b/docs/figures/phase/11129442_kiel.png differ
diff --git a/docs/figures/preview_interp_MS/16CygA_interp_preview_across_resolution.png b/docs/figures/preview_interp_MS/16CygA_interp_preview_across_resolution.png
new file mode 100644
index 0000000..76f4596
Binary files /dev/null and b/docs/figures/preview_interp_MS/16CygA_interp_preview_across_resolution.png differ
diff --git a/docs/figures/preview_interp_MS/16CygA_interp_preview_along_resolution.png b/docs/figures/preview_interp_MS/16CygA_interp_preview_along_resolution.png
new file mode 100644
index 0000000..ab35656
Binary files /dev/null and b/docs/figures/preview_interp_MS/16CygA_interp_preview_along_resolution.png differ
diff --git a/docs/figures/ratios/16CygA_corner.png b/docs/figures/ratios/16CygA_corner.png
new file mode 100644
index 0000000..6f34db6
Binary files /dev/null and b/docs/figures/ratios/16CygA_corner.png differ
diff --git a/docs/figures/ratios/16CygA_kiel.png b/docs/figures/ratios/16CygA_kiel.png
new file mode 100644
index 0000000..f9c09a5
Binary files /dev/null and b/docs/figures/ratios/16CygA_kiel.png differ
diff --git a/docs/figures/ratios/16CygA_ratios_r012.png b/docs/figures/ratios/16CygA_ratios_r012.png
new file mode 100644
index 0000000..1c61f04
Binary files /dev/null and b/docs/figures/ratios/16CygA_ratios_r012.png differ
diff --git a/docs/figures/ratios/16CygA_ratios_r012_cormap.png b/docs/figures/ratios/16CygA_ratios_r012_cormap.png
new file mode 100644
index 0000000..9c1caf6
Binary files /dev/null and b/docs/figures/ratios/16CygA_ratios_r012_cormap.png differ
diff --git a/docs/figures/subgiant/Valid_245_corner.png b/docs/figures/subgiant/Valid_245_corner.png
new file mode 100644
index 0000000..828163e
Binary files /dev/null and b/docs/figures/subgiant/Valid_245_corner.png differ
diff --git a/docs/figures/subgiant/Valid_245_dupechelle.png b/docs/figures/subgiant/Valid_245_dupechelle.png
new file mode 100644
index 0000000..a89a33c
Binary files /dev/null and b/docs/figures/subgiant/Valid_245_dupechelle.png differ
diff --git a/docs/figures/subgiant/Valid_245_dupechelle_uncorrected.png b/docs/figures/subgiant/Valid_245_dupechelle_uncorrected.png
new file mode 100644
index 0000000..ccca79a
Binary files /dev/null and b/docs/figures/subgiant/Valid_245_dupechelle_uncorrected.png differ
diff --git a/docs/figures/subgiant/Valid_245_echelle.png b/docs/figures/subgiant/Valid_245_echelle.png
new file mode 100644
index 0000000..ea574a0
Binary files /dev/null and b/docs/figures/subgiant/Valid_245_echelle.png differ
diff --git a/docs/figures/subgiant/Valid_245_echelle_uncorrected.png b/docs/figures/subgiant/Valid_245_echelle_uncorrected.png
new file mode 100644
index 0000000..862523d
Binary files /dev/null and b/docs/figures/subgiant/Valid_245_echelle_uncorrected.png differ
diff --git a/docs/figures/subgiant/Valid_245_kiel.png b/docs/figures/subgiant/Valid_245_kiel.png
new file mode 100644
index 0000000..fe69bad
Binary files /dev/null and b/docs/figures/subgiant/Valid_245_kiel.png differ
diff --git a/docs/figures/subgiant/Valid_245_pairechelle.png b/docs/figures/subgiant/Valid_245_pairechelle.png
new file mode 100644
index 0000000..cca5398
Binary files /dev/null and b/docs/figures/subgiant/Valid_245_pairechelle.png differ
diff --git a/docs/figures/subgiant/Valid_245_pairechelle_uncorrected.png b/docs/figures/subgiant/Valid_245_pairechelle_uncorrected.png
new file mode 100644
index 0000000..1bc5ede
Binary files /dev/null and b/docs/figures/subgiant/Valid_245_pairechelle_uncorrected.png differ
diff --git a/docs/grids.rst b/docs/grids.rst
index 52cb7aa..ad38861 100644
--- a/docs/grids.rst
+++ b/docs/grids.rst
@@ -3,7 +3,7 @@
 Grids of models
 ===============
 
-BASTA runs over grids of stellar tracks or isochrones stored in hierarchical data format ``hdf5``. The list of
+BASTA uses grids of stellar tracks or isochrones stored in hierarchical data format ``hdf5``. The list of
 parameters included in our grids can be seen in :py:meth:`constants.parameters`, but one should keep in mind that
 not all parameters are included in all grids. As an example, isochrones do not contain individual frequencies of
 oscillations while stellar tracks do, and fits to e.g., frequency ratios can only be performed with grids of stellar
diff --git a/docs/index.rst b/docs/index.rst
index 0f9fcde..4df628b 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -3,13 +3,11 @@ BASTA: The BAyesian STellar Algorithm
 
 **Welcome to BASTA's documentation!**
 
-*We are currently having issues with the examples pages (images are not showing in Firefox and Chrome). We are working on a release to fix it. Until then, please try with Safari. We are sorry for the inconvenience!*
-
-BASTA is a python-based fitting tool designed to determine properties of stars using a pre-computed grid of stellar
+BASTA is a Python-based fitting tool designed to determine properties of stars using a pre-computed grid of stellar
 models. It calculates the probability density function of a given stellar property based on a set of observational
 constraints defined by the user.
 
-Current stable version: v1.2.0
+Current stable version: v1.2.1
 
 *Please follow the repository on* `GitHub <https://github.com/BASTAcode/BASTA>`_ *to get notifications on new releases: Click "Watch" then "Custom" and tick "Releases.*
 
@@ -55,7 +53,7 @@ The original author of the code is:
 
 The current core developing team are:
 
-* Jakob Lysgaard Rørsted
+* Jakob Lysgaard Rørsted (maintainer)
 * Mark Lykke Winther
 * Amalie Stokholm
 * Kuldeep Verma
diff --git a/docs/install.rst b/docs/install.rst
index ea20caf..b8cd416 100644
--- a/docs/install.rst
+++ b/docs/install.rst
@@ -34,10 +34,9 @@ Now, assuming you have downloaded the code, run the following to setup the virtu
     cd BASTA
     python3 -m venv venv
     source venv/bin/activate
-    pip install --upgrade pip
+    pip install --upgrade pip setuptools wheel
     deactivate
     source venv/bin/activate
-    pip install wheel
     pip install -r requirements.txt
     deactivate
 
@@ -66,6 +65,8 @@ If your installation of BASTA is in a folder different than the default ``~/BAST
 Consistent f2py
 ---------------
 
+The following should not be necessary if you use *pyenv* or similar to manage your Python environment.
+
 BASTA needs the executable ``f2py3`` to exist in the path to be able to compile external Fortran routines required for
 different parts of the code. In order to avoid issues, this tool *must* match the version of NumPy used to run BASTA.
 It is possible to ensure this by utilizing the version installed in the virtual environment (since the executable is
diff --git a/docs/prepare_doc_figures.py b/docs/prepare_doc_figures.py
new file mode 100644
index 0000000..4556968
--- /dev/null
+++ b/docs/prepare_doc_figures.py
@@ -0,0 +1,42 @@
+# Convert reference figures from the examples to png to be rendered by the documentation
+import os
+import subprocess
+from contextlib import contextmanager
+
+from tqdm import tqdm
+
+
+@contextmanager
+def cd(newdir):
+    """Change directory"""
+    prevdir = os.getcwd()
+    os.chdir(os.path.expanduser(newdir))
+    try:
+        yield
+    finally:
+        os.chdir(prevdir)
+
+
+# Settings for imagemagic (change name on Windows)
+imname = ["convert"]
+settings = ["-density", "96", "-quality", "85"]
+imgcmd = imname + settings
+
+# Paths
+indir = os.path.abspath("../examples/reference")
+outdir = os.path.abspath("figures")
+
+# Transverse all example cases and process all pdf -> png, and store in output dir
+with cd(indir):
+    cases = next(os.walk("."))[1]
+    for case in tqdm(cases):
+        os.makedirs(os.path.join(outdir, case), exist_ok=True)
+        with cd(os.path.join(indir, case)):
+            files = next(os.walk("."))[2]
+            for file in files:
+                if not file.endswith(".pdf"):
+                    continue
+                outname = file[:-4] + ".png"
+                outfile = os.path.join(outdir, case, outname)
+                imgcall = imgcmd + [f"{file}", f"{outfile}"]
+                subprocess.check_call(imgcall)
diff --git a/examples/create_inputfile.py b/examples/create_inputfile.py
index 4b68b03..424cb70 100644
--- a/examples/create_inputfile.py
+++ b/examples/create_inputfile.py
@@ -222,6 +222,15 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     #                    * "1/N-dof": Normalisation by number of frequencies minus the
     #                                 (user-specified) degrees of freedom. This option
     #                                 must be supplemented by a specification of "dof".
+    #
+    # - "threepoint": Set 'True' to change to the three-point small frequency separation
+    #                 formulation for the frequency ratios, instead of the five-point
+    #                 small frequency separation. The default, if not set or set as
+    #                 'False', is to use the five-point formulation.
+    #
+    # - "dnuprior": Set to 'False' to disable the automatic prior on dnu (mimicking the
+    #               range defined by dnufrac), which speeds up the computation of the
+    #               fit. Recommended to keep as 'True' (default).
 
     # Example of typical settings for a frequency fit (with default seismic weights):
     # define_fit["freqparams"] = {
diff --git a/examples/preview_interpolation.py b/examples/preview_interpolation.py
index 3bd7998..7579121 100644
--- a/examples/preview_interpolation.py
+++ b/examples/preview_interpolation.py
@@ -219,7 +219,7 @@ def test_across_interpolation(
         sobol = resolution["scale"]
 
     if len(baseparams) == 0:
-        baseparams = [par.decode("UTF-8") for par in grid["header/active_weights"]]
+        baseparams = [par.decode("UTF-8") for par in grid["header/pars_sampled"]]
 
     base = np.zeros((len(selectedmodels), len(baseparams)))
     for i, name in enumerate(selectedmodels):
@@ -227,7 +227,7 @@ def test_across_interpolation(
             parm = grid[os.path.join(bpath, name)][bpar][0]
             base[i, j] = parm
     tri = spatial.Delaunay(base)
-    newbase, _ = iac._calc_across_points(base, baseparams, tri, sobol, outname)
+    newbase, _, _ = iac._calc_sobol_points(base, baseparams, tri, sobol, outname)
     success = ip.base_corner(baseparams, base, newbase, tri, sobol, outname)
 
 
diff --git a/examples/reference/distance/11129442_corner.pdf b/examples/reference/distance/11129442_corner.pdf
index b5c2421..77d0cf6 100644
Binary files a/examples/reference/distance/11129442_corner.pdf and b/examples/reference/distance/11129442_corner.pdf differ
diff --git a/examples/reference/distance/11129442_distance_corner.pdf b/examples/reference/distance/11129442_distance_corner.pdf
index fe02c69..9cc4904 100644
Binary files a/examples/reference/distance/11129442_distance_corner.pdf and b/examples/reference/distance/11129442_distance_corner.pdf differ
diff --git a/examples/reference/distance/11129442_kiel.pdf b/examples/reference/distance/11129442_kiel.pdf
index c450043..4957358 100644
Binary files a/examples/reference/distance/11129442_kiel.pdf and b/examples/reference/distance/11129442_kiel.pdf differ
diff --git a/examples/reference/epsilon/16CygA_corner.pdf b/examples/reference/epsilon/16CygA_corner.pdf
index fccf013..707212a 100644
Binary files a/examples/reference/epsilon/16CygA_corner.pdf and b/examples/reference/epsilon/16CygA_corner.pdf differ
diff --git a/examples/reference/epsilon/16CygA_dupechelle.pdf b/examples/reference/epsilon/16CygA_dupechelle.pdf
index 3ea6bac..cd24835 100644
Binary files a/examples/reference/epsilon/16CygA_dupechelle.pdf and b/examples/reference/epsilon/16CygA_dupechelle.pdf differ
diff --git a/examples/reference/epsilon/16CygA_dupechelle_uncorrected.pdf b/examples/reference/epsilon/16CygA_dupechelle_uncorrected.pdf
index 1774078..4891bac 100644
Binary files a/examples/reference/epsilon/16CygA_dupechelle_uncorrected.pdf and b/examples/reference/epsilon/16CygA_dupechelle_uncorrected.pdf differ
diff --git a/examples/reference/epsilon/16CygA_echelle.pdf b/examples/reference/epsilon/16CygA_echelle.pdf
index c09332a..f5cc8b5 100644
Binary files a/examples/reference/epsilon/16CygA_echelle.pdf and b/examples/reference/epsilon/16CygA_echelle.pdf differ
diff --git a/examples/reference/epsilon/16CygA_echelle_uncorrected.pdf b/examples/reference/epsilon/16CygA_echelle_uncorrected.pdf
index 7a37209..0c759d8 100644
Binary files a/examples/reference/epsilon/16CygA_echelle_uncorrected.pdf and b/examples/reference/epsilon/16CygA_echelle_uncorrected.pdf differ
diff --git a/examples/reference/epsilon/16CygA_epsdiff_e012.pdf b/examples/reference/epsilon/16CygA_epsdiff_e012.pdf
index f71a042..93b4050 100644
Binary files a/examples/reference/epsilon/16CygA_epsdiff_e012.pdf and b/examples/reference/epsilon/16CygA_epsdiff_e012.pdf differ
diff --git a/examples/reference/epsilon/16CygA_epsdiff_e012_cormap.pdf b/examples/reference/epsilon/16CygA_epsdiff_e012_cormap.pdf
index 4138216..0d03a17 100644
Binary files a/examples/reference/epsilon/16CygA_epsdiff_e012_cormap.pdf and b/examples/reference/epsilon/16CygA_epsdiff_e012_cormap.pdf differ
diff --git a/examples/reference/epsilon/16CygA_kiel.pdf b/examples/reference/epsilon/16CygA_kiel.pdf
index 4705ea6..9420d6e 100644
Binary files a/examples/reference/epsilon/16CygA_kiel.pdf and b/examples/reference/epsilon/16CygA_kiel.pdf differ
diff --git a/examples/reference/epsilon/16CygA_pairechelle.pdf b/examples/reference/epsilon/16CygA_pairechelle.pdf
index 2624983..464787a 100644
Binary files a/examples/reference/epsilon/16CygA_pairechelle.pdf and b/examples/reference/epsilon/16CygA_pairechelle.pdf differ
diff --git a/examples/reference/epsilon/16CygA_pairechelle_uncorrected.pdf b/examples/reference/epsilon/16CygA_pairechelle_uncorrected.pdf
index de0f649..aa10ba0 100644
Binary files a/examples/reference/epsilon/16CygA_pairechelle_uncorrected.pdf and b/examples/reference/epsilon/16CygA_pairechelle_uncorrected.pdf differ
diff --git a/examples/reference/epsilon/16CygA_ratios_r012.pdf b/examples/reference/epsilon/16CygA_ratios_r012.pdf
index 6bb7346..1405ed7 100644
Binary files a/examples/reference/epsilon/16CygA_ratios_r012.pdf and b/examples/reference/epsilon/16CygA_ratios_r012.pdf differ
diff --git a/examples/reference/epsilon/16CygA_ratios_r012_cormap.pdf b/examples/reference/epsilon/16CygA_ratios_r012_cormap.pdf
index 085b8dc..67902f2 100644
Binary files a/examples/reference/epsilon/16CygA_ratios_r012_cormap.pdf and b/examples/reference/epsilon/16CygA_ratios_r012_cormap.pdf differ
diff --git a/examples/reference/freqs/16CygA_corner.pdf b/examples/reference/freqs/16CygA_corner.pdf
index 5061481..134e95b 100644
Binary files a/examples/reference/freqs/16CygA_corner.pdf and b/examples/reference/freqs/16CygA_corner.pdf differ
diff --git a/examples/reference/freqs/16CygA_dupechelle.pdf b/examples/reference/freqs/16CygA_dupechelle.pdf
index fc5ea54..0be9ac9 100644
Binary files a/examples/reference/freqs/16CygA_dupechelle.pdf and b/examples/reference/freqs/16CygA_dupechelle.pdf differ
diff --git a/examples/reference/freqs/16CygA_dupechelle_uncorrected.pdf b/examples/reference/freqs/16CygA_dupechelle_uncorrected.pdf
index f57ac1a..7b3e5ed 100644
Binary files a/examples/reference/freqs/16CygA_dupechelle_uncorrected.pdf and b/examples/reference/freqs/16CygA_dupechelle_uncorrected.pdf differ
diff --git a/examples/reference/freqs/16CygA_echelle.pdf b/examples/reference/freqs/16CygA_echelle.pdf
index 9882186..63ad056 100644
Binary files a/examples/reference/freqs/16CygA_echelle.pdf and b/examples/reference/freqs/16CygA_echelle.pdf differ
diff --git a/examples/reference/freqs/16CygA_echelle_uncorrected.pdf b/examples/reference/freqs/16CygA_echelle_uncorrected.pdf
index 933d59f..4662dbc 100644
Binary files a/examples/reference/freqs/16CygA_echelle_uncorrected.pdf and b/examples/reference/freqs/16CygA_echelle_uncorrected.pdf differ
diff --git a/examples/reference/freqs/16CygA_kiel.pdf b/examples/reference/freqs/16CygA_kiel.pdf
index 80a2d61..a92f08c 100644
Binary files a/examples/reference/freqs/16CygA_kiel.pdf and b/examples/reference/freqs/16CygA_kiel.pdf differ
diff --git a/examples/reference/freqs/16CygA_pairechelle.pdf b/examples/reference/freqs/16CygA_pairechelle.pdf
index aba2542..a858c47 100644
Binary files a/examples/reference/freqs/16CygA_pairechelle.pdf and b/examples/reference/freqs/16CygA_pairechelle.pdf differ
diff --git a/examples/reference/freqs/16CygA_pairechelle_uncorrected.pdf b/examples/reference/freqs/16CygA_pairechelle_uncorrected.pdf
index c653d47..ff0400d 100644
Binary files a/examples/reference/freqs/16CygA_pairechelle_uncorrected.pdf and b/examples/reference/freqs/16CygA_pairechelle_uncorrected.pdf differ
diff --git a/examples/reference/glitches/16CygA_corner.pdf b/examples/reference/glitches/16CygA_corner.pdf
index 1b4ec7b..fb685ae 100644
Binary files a/examples/reference/glitches/16CygA_corner.pdf and b/examples/reference/glitches/16CygA_corner.pdf differ
diff --git a/examples/reference/glitches/16CygA_dupechelle.pdf b/examples/reference/glitches/16CygA_dupechelle.pdf
index 4b85a7c..5d6c6a4 100644
Binary files a/examples/reference/glitches/16CygA_dupechelle.pdf and b/examples/reference/glitches/16CygA_dupechelle.pdf differ
diff --git a/examples/reference/glitches/16CygA_dupechelle_uncorrected.pdf b/examples/reference/glitches/16CygA_dupechelle_uncorrected.pdf
index 091f426..31c37d3 100644
Binary files a/examples/reference/glitches/16CygA_dupechelle_uncorrected.pdf and b/examples/reference/glitches/16CygA_dupechelle_uncorrected.pdf differ
diff --git a/examples/reference/glitches/16CygA_echelle.pdf b/examples/reference/glitches/16CygA_echelle.pdf
index f3e29e5..995eb8a 100644
Binary files a/examples/reference/glitches/16CygA_echelle.pdf and b/examples/reference/glitches/16CygA_echelle.pdf differ
diff --git a/examples/reference/glitches/16CygA_echelle_uncorrected.pdf b/examples/reference/glitches/16CygA_echelle_uncorrected.pdf
index 3f266dc..b505961 100644
Binary files a/examples/reference/glitches/16CygA_echelle_uncorrected.pdf and b/examples/reference/glitches/16CygA_echelle_uncorrected.pdf differ
diff --git a/examples/reference/glitches/16CygA_kiel.pdf b/examples/reference/glitches/16CygA_kiel.pdf
index 55c27e0..cd7fbc6 100644
Binary files a/examples/reference/glitches/16CygA_kiel.pdf and b/examples/reference/glitches/16CygA_kiel.pdf differ
diff --git a/examples/reference/glitches/16CygA_pairechelle.pdf b/examples/reference/glitches/16CygA_pairechelle.pdf
index 95a90b2..b6fa21a 100644
Binary files a/examples/reference/glitches/16CygA_pairechelle.pdf and b/examples/reference/glitches/16CygA_pairechelle.pdf differ
diff --git a/examples/reference/glitches/16CygA_pairechelle_uncorrected.pdf b/examples/reference/glitches/16CygA_pairechelle_uncorrected.pdf
index 4a418b0..c72ccd6 100644
Binary files a/examples/reference/glitches/16CygA_pairechelle_uncorrected.pdf and b/examples/reference/glitches/16CygA_pairechelle_uncorrected.pdf differ
diff --git a/examples/reference/global/16CygA_corner.pdf b/examples/reference/global/16CygA_corner.pdf
index ded2e30..2f7ec02 100644
Binary files a/examples/reference/global/16CygA_corner.pdf and b/examples/reference/global/16CygA_corner.pdf differ
diff --git a/examples/reference/global/16CygA_kiel.pdf b/examples/reference/global/16CygA_kiel.pdf
index 7be180e..c697a04 100644
Binary files a/examples/reference/global/16CygA_kiel.pdf and b/examples/reference/global/16CygA_kiel.pdf differ
diff --git a/examples/reference/interp_MS/16CygA_corner.pdf b/examples/reference/interp_MS/16CygA_corner.pdf
index e176fb3..a03fbfd 100644
Binary files a/examples/reference/interp_MS/16CygA_corner.pdf and b/examples/reference/interp_MS/16CygA_corner.pdf differ
diff --git a/examples/reference/interp_MS/16CygA_dupechelle.pdf b/examples/reference/interp_MS/16CygA_dupechelle.pdf
index 7020738..4b91678 100644
Binary files a/examples/reference/interp_MS/16CygA_dupechelle.pdf and b/examples/reference/interp_MS/16CygA_dupechelle.pdf differ
diff --git a/examples/reference/interp_MS/16CygA_dupechelle_uncorrected.pdf b/examples/reference/interp_MS/16CygA_dupechelle_uncorrected.pdf
index e09467f..5fff86d 100644
Binary files a/examples/reference/interp_MS/16CygA_dupechelle_uncorrected.pdf and b/examples/reference/interp_MS/16CygA_dupechelle_uncorrected.pdf differ
diff --git a/examples/reference/interp_MS/16CygA_echelle.pdf b/examples/reference/interp_MS/16CygA_echelle.pdf
index 0b9e17a..5c53124 100644
Binary files a/examples/reference/interp_MS/16CygA_echelle.pdf and b/examples/reference/interp_MS/16CygA_echelle.pdf differ
diff --git a/examples/reference/interp_MS/16CygA_echelle_uncorrected.pdf b/examples/reference/interp_MS/16CygA_echelle_uncorrected.pdf
index f01a288..47a1ad1 100644
Binary files a/examples/reference/interp_MS/16CygA_echelle_uncorrected.pdf and b/examples/reference/interp_MS/16CygA_echelle_uncorrected.pdf differ
diff --git a/examples/reference/interp_MS/16CygA_kiel.pdf b/examples/reference/interp_MS/16CygA_kiel.pdf
index 3fb8042..726abc6 100644
Binary files a/examples/reference/interp_MS/16CygA_kiel.pdf and b/examples/reference/interp_MS/16CygA_kiel.pdf differ
diff --git a/examples/reference/interp_MS/16CygA_pairechelle.pdf b/examples/reference/interp_MS/16CygA_pairechelle.pdf
index 7e23b84..32609cc 100644
Binary files a/examples/reference/interp_MS/16CygA_pairechelle.pdf and b/examples/reference/interp_MS/16CygA_pairechelle.pdf differ
diff --git a/examples/reference/interp_MS/16CygA_pairechelle_uncorrected.pdf b/examples/reference/interp_MS/16CygA_pairechelle_uncorrected.pdf
index 8439115..b6a4fcf 100644
Binary files a/examples/reference/interp_MS/16CygA_pairechelle_uncorrected.pdf and b/examples/reference/interp_MS/16CygA_pairechelle_uncorrected.pdf differ
diff --git a/examples/reference/interp_MS/intpol_example_16CygA_intpolbase.pdf b/examples/reference/interp_MS/intpol_example_16CygA_intpolbase.pdf
index e9cd925..06b3cf7 100644
Binary files a/examples/reference/interp_MS/intpol_example_16CygA_intpolbase.pdf and b/examples/reference/interp_MS/intpol_example_16CygA_intpolbase.pdf differ
diff --git a/examples/reference/json/16CygA_corner.pdf b/examples/reference/json/16CygA_corner.pdf
index 7874810..fa33f56 100644
Binary files a/examples/reference/json/16CygA_corner.pdf and b/examples/reference/json/16CygA_corner.pdf differ
diff --git a/examples/reference/json/16CygA_kiel.pdf b/examples/reference/json/16CygA_kiel.pdf
index 1430b55..f46752a 100644
Binary files a/examples/reference/json/16CygA_kiel.pdf and b/examples/reference/json/16CygA_kiel.pdf differ
diff --git a/examples/reference/parallax/11129442_corner.pdf b/examples/reference/parallax/11129442_corner.pdf
index bba55bb..3e4c57e 100644
Binary files a/examples/reference/parallax/11129442_corner.pdf and b/examples/reference/parallax/11129442_corner.pdf differ
diff --git a/examples/reference/parallax/11129442_distance_corner.pdf b/examples/reference/parallax/11129442_distance_corner.pdf
index b435e5d..b799cf4 100644
Binary files a/examples/reference/parallax/11129442_distance_corner.pdf and b/examples/reference/parallax/11129442_distance_corner.pdf differ
diff --git a/examples/reference/parallax/11129442_kiel.pdf b/examples/reference/parallax/11129442_kiel.pdf
index 2dc0da0..c476a7c 100644
Binary files a/examples/reference/parallax/11129442_kiel.pdf and b/examples/reference/parallax/11129442_kiel.pdf differ
diff --git a/examples/reference/phase/11129442_corner.pdf b/examples/reference/phase/11129442_corner.pdf
index 2784c99..c5eb40c 100644
Binary files a/examples/reference/phase/11129442_corner.pdf and b/examples/reference/phase/11129442_corner.pdf differ
diff --git a/examples/reference/phase/11129442_kiel.pdf b/examples/reference/phase/11129442_kiel.pdf
index b0b4e84..0d63caf 100644
Binary files a/examples/reference/phase/11129442_kiel.pdf and b/examples/reference/phase/11129442_kiel.pdf differ
diff --git a/examples/reference/preview_interp_MS/16CygA_interp_preview_across_resolution.pdf b/examples/reference/preview_interp_MS/16CygA_interp_preview_across_resolution.pdf
index e6180a1..9e2da4a 100644
Binary files a/examples/reference/preview_interp_MS/16CygA_interp_preview_across_resolution.pdf and b/examples/reference/preview_interp_MS/16CygA_interp_preview_across_resolution.pdf differ
diff --git a/examples/reference/preview_interp_MS/16CygA_interp_preview_along_resolution.pdf b/examples/reference/preview_interp_MS/16CygA_interp_preview_along_resolution.pdf
index 7db2889..97b803a 100644
Binary files a/examples/reference/preview_interp_MS/16CygA_interp_preview_along_resolution.pdf and b/examples/reference/preview_interp_MS/16CygA_interp_preview_along_resolution.pdf differ
diff --git a/examples/reference/ratios/16CygA_corner.pdf b/examples/reference/ratios/16CygA_corner.pdf
index cc259a0..a4ab530 100644
Binary files a/examples/reference/ratios/16CygA_corner.pdf and b/examples/reference/ratios/16CygA_corner.pdf differ
diff --git a/examples/reference/ratios/16CygA_kiel.pdf b/examples/reference/ratios/16CygA_kiel.pdf
index 29a0b3e..520fd87 100644
Binary files a/examples/reference/ratios/16CygA_kiel.pdf and b/examples/reference/ratios/16CygA_kiel.pdf differ
diff --git a/examples/reference/ratios/16CygA_ratios_r012.pdf b/examples/reference/ratios/16CygA_ratios_r012.pdf
index 7f39f18..be68dbe 100644
Binary files a/examples/reference/ratios/16CygA_ratios_r012.pdf and b/examples/reference/ratios/16CygA_ratios_r012.pdf differ
diff --git a/examples/reference/ratios/16CygA_ratios_r012_cormap.pdf b/examples/reference/ratios/16CygA_ratios_r012_cormap.pdf
index 4c7702a..f4a569b 100644
Binary files a/examples/reference/ratios/16CygA_ratios_r012_cormap.pdf and b/examples/reference/ratios/16CygA_ratios_r012_cormap.pdf differ
diff --git a/examples/reference/subgiant/Valid_245_corner.pdf b/examples/reference/subgiant/Valid_245_corner.pdf
index 2ffc6ac..a009544 100644
Binary files a/examples/reference/subgiant/Valid_245_corner.pdf and b/examples/reference/subgiant/Valid_245_corner.pdf differ
diff --git a/examples/reference/subgiant/Valid_245_dupechelle.pdf b/examples/reference/subgiant/Valid_245_dupechelle.pdf
index d611788..cdad34c 100644
Binary files a/examples/reference/subgiant/Valid_245_dupechelle.pdf and b/examples/reference/subgiant/Valid_245_dupechelle.pdf differ
diff --git a/examples/reference/subgiant/Valid_245_dupechelle_uncorrected.pdf b/examples/reference/subgiant/Valid_245_dupechelle_uncorrected.pdf
index b7d42e5..55f1f21 100644
Binary files a/examples/reference/subgiant/Valid_245_dupechelle_uncorrected.pdf and b/examples/reference/subgiant/Valid_245_dupechelle_uncorrected.pdf differ
diff --git a/examples/reference/subgiant/Valid_245_echelle.pdf b/examples/reference/subgiant/Valid_245_echelle.pdf
index 2f26d2b..7e885a4 100644
Binary files a/examples/reference/subgiant/Valid_245_echelle.pdf and b/examples/reference/subgiant/Valid_245_echelle.pdf differ
diff --git a/examples/reference/subgiant/Valid_245_echelle_uncorrected.pdf b/examples/reference/subgiant/Valid_245_echelle_uncorrected.pdf
index 981b04b..f32287c 100644
Binary files a/examples/reference/subgiant/Valid_245_echelle_uncorrected.pdf and b/examples/reference/subgiant/Valid_245_echelle_uncorrected.pdf differ
diff --git a/examples/reference/subgiant/Valid_245_kiel.pdf b/examples/reference/subgiant/Valid_245_kiel.pdf
index 923b77f..416d013 100644
Binary files a/examples/reference/subgiant/Valid_245_kiel.pdf and b/examples/reference/subgiant/Valid_245_kiel.pdf differ
diff --git a/examples/reference/subgiant/Valid_245_pairechelle.pdf b/examples/reference/subgiant/Valid_245_pairechelle.pdf
index 108b159..19092c4 100644
Binary files a/examples/reference/subgiant/Valid_245_pairechelle.pdf and b/examples/reference/subgiant/Valid_245_pairechelle.pdf differ
diff --git a/examples/reference/subgiant/Valid_245_pairechelle_uncorrected.pdf b/examples/reference/subgiant/Valid_245_pairechelle_uncorrected.pdf
index 238d764..a932cc7 100644
Binary files a/examples/reference/subgiant/Valid_245_pairechelle_uncorrected.pdf and b/examples/reference/subgiant/Valid_245_pairechelle_uncorrected.pdf differ
diff --git a/examples/xmlinput/create_inputfile_distance.py b/examples/xmlinput/create_inputfile_distance.py
index f4b996a..6fd82d1 100644
--- a/examples/xmlinput/create_inputfile_distance.py
+++ b/examples/xmlinput/create_inputfile_distance.py
@@ -231,14 +231,18 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     #                    * "1/N-dof": Normalisation by number of frequencies minus the
     #                                 (user-specified) degrees of freedom. This option
     #                                 must be supplemented by a specification of "dof".
+    #
     # - "threepoint": Set 'True' to change to the three-point small frequency separation
     #                 formulation for the frequency ratios, instead of the five-point
     #                 small frequency separation. The default, if not set or set as
     #                 'False', is to use the five-point formulation.
+    #
+    # - "dnuprior": Set to 'False' to disable the automatic prior on dnu (mimicking the
+    #               range defined by dnufrac), which speeds up the computation of the
+    #               fit. Recommended to keep as 'True' (default).
 
     # Example of typical settings for a frequency fit (with default seismic weights):
     # define_fit["freqparams"] = {
-    # define_fit["freqparams"] = {
     #     "freqpath": os.path.join(BASTADIR, "examples/data/freqs"),
     #     "fcor": "BG14",
     #     "correlations": False,
@@ -287,6 +291,8 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     define_output["outparams"] = (
         "Teff",
         "FeH",
+        "dnuSer",
+        "numax",
         "radPhot",
         "massfin",
         "age",
diff --git a/examples/xmlinput/create_inputfile_epsilondifference.py b/examples/xmlinput/create_inputfile_epsilondifference.py
index feee5ed..92f78d4 100644
--- a/examples/xmlinput/create_inputfile_epsilondifference.py
+++ b/examples/xmlinput/create_inputfile_epsilondifference.py
@@ -224,13 +224,17 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     #                    * "1/N-dof": Normalisation by number of frequencies minus the
     #                                 (user-specified) degrees of freedom. This option
     #                                 must be supplemented by a specification of "dof".
+    #
     # - "threepoint": Set 'True' to change to the three-point small frequency separation
     #                 formulation for the frequency ratios, instead of the five-point
     #                 small frequency separation. The default, if not set or set as
     #                 'False', is to use the five-point formulation.
+    #
+    # - "dnuprior": Set to 'False' to disable the automatic prior on dnu (mimicking the
+    #               range defined by dnufrac), which speeds up the computation of the
+    #               fit. Recommended to keep as 'True' (default).
 
     # Example of typical settings for a frequency fit (with default seismic weights):
-    # define_fit["freqparams"] = {
     define_fit["freqparams"] = {
         "freqpath": os.path.join(BASTADIR, "examples/data/freqs"),
         "fcor": "BG14",
@@ -277,7 +281,15 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     # --> The names must match entries in the parameter list (basta/constants.py)
     # --> A reasonable choice is to (as a minimum) output the parameters used in the fit
     # --> If you want to predict distance, add the special keyword "distance".
-    define_output["outparams"] = ("Teff", "FeH", "radPhot", "massfin", "age")
+    define_output["outparams"] = (
+        "Teff",
+        "FeH",
+        "dnufit",
+        "numax",
+        "radPhot",
+        "massfin",
+        "age",
+    )
 
     # Name of the output file containing the results of the fit in ascii format.
     # --> A version in xml-format will be automatically created
diff --git a/examples/xmlinput/create_inputfile_freqs.py b/examples/xmlinput/create_inputfile_freqs.py
index 102eaa5..395e6b8 100644
--- a/examples/xmlinput/create_inputfile_freqs.py
+++ b/examples/xmlinput/create_inputfile_freqs.py
@@ -223,13 +223,17 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     #                    * "1/N-dof": Normalisation by number of frequencies minus the
     #                                 (user-specified) degrees of freedom. This option
     #                                 must be supplemented by a specification of "dof".
+    #
     # - "threepoint": Set 'True' to change to the three-point small frequency separation
     #                 formulation for the frequency ratios, instead of the five-point
     #                 small frequency separation. The default, if not set or set as
     #                 'False', is to use the five-point formulation.
+    #
+    # - "dnuprior": Set to 'False' to disable the automatic prior on dnu (mimicking the
+    #               range defined by dnufrac), which speeds up the computation of the
+    #               fit. Recommended to keep as 'True' (default).
 
     # Example of typical settings for a frequency fit (with default seismic weights):
-    # define_fit["freqparams"] = {
     define_fit["freqparams"] = {
         "freqpath": os.path.join(BASTADIR, "examples/data/freqs"),
         "fcor": "BG14",
@@ -276,7 +280,15 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     # --> The names must match entries in the parameter list (basta/constants.py)
     # --> A reasonable choice is to (as a minimum) output the parameters used in the fit
     # --> If you want to predict distance, add the special keyword "distance".
-    define_output["outparams"] = ("Teff", "FeH", "radPhot", "massfin", "age")
+    define_output["outparams"] = (
+        "Teff",
+        "FeH",
+        "dnufit",
+        "numax",
+        "radPhot",
+        "massfin",
+        "age",
+    )
 
     # Name of the output file containing the results of the fit in ascii format.
     # --> A version in xml-format will be automatically created
diff --git a/examples/xmlinput/create_inputfile_glitches.py b/examples/xmlinput/create_inputfile_glitches.py
index 0534ba1..b958bf0 100644
--- a/examples/xmlinput/create_inputfile_glitches.py
+++ b/examples/xmlinput/create_inputfile_glitches.py
@@ -223,13 +223,17 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     #                    * "1/N-dof": Normalisation by number of frequencies minus the
     #                                 (user-specified) degrees of freedom. This option
     #                                 must be supplemented by a specification of "dof".
+    #
     # - "threepoint": Set 'True' to change to the three-point small frequency separation
     #                 formulation for the frequency ratios, instead of the five-point
     #                 small frequency separation. The default, if not set or set as
     #                 'False', is to use the five-point formulation.
+    #
+    # - "dnuprior": Set to 'False' to disable the automatic prior on dnu (mimicking the
+    #               range defined by dnufrac), which speeds up the computation of the
+    #               fit. Recommended to keep as 'True' (default).
 
     # Example of typical settings for a frequency fit (with default seismic weights):
-    # define_fit["freqparams"] = {
     define_fit["freqparams"] = {
         "freqpath": os.path.join(BASTADIR, "examples/data/freqs"),
         "fcor": "BG14",
@@ -276,7 +280,15 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     # --> The names must match entries in the parameter list (basta/constants.py)
     # --> A reasonable choice is to (as a minimum) output the parameters used in the fit
     # --> If you want to predict distance, add the special keyword "distance".
-    define_output["outparams"] = ("Teff", "FeH", "radPhot", "massfin", "age")
+    define_output["outparams"] = (
+        "Teff",
+        "FeH",
+        "dnufit",
+        "numax",
+        "radPhot",
+        "massfin",
+        "age",
+    )
 
     # Name of the output file containing the results of the fit in ascii format.
     # --> A version in xml-format will be automatically created
diff --git a/examples/xmlinput/create_inputfile_global.py b/examples/xmlinput/create_inputfile_global.py
index c612c03..3ed2181 100644
--- a/examples/xmlinput/create_inputfile_global.py
+++ b/examples/xmlinput/create_inputfile_global.py
@@ -224,10 +224,15 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     #                    * "1/N-dof": Normalisation by number of frequencies minus the
     #                                 (user-specified) degrees of freedom. This option
     #                                 must be supplemented by a specification of "dof".
+    #
     # - "threepoint": Set 'True' to change to the three-point small frequency separation
     #                 formulation for the frequency ratios, instead of the five-point
     #                 small frequency separation. The default, if not set or set as
     #                 'False', is to use the five-point formulation.
+    #
+    # - "dnuprior": Set to 'False' to disable the automatic prior on dnu (mimicking the
+    #               range defined by dnufrac), which speeds up the computation of the
+    #               fit. Recommended to keep as 'True' (default).
 
     # Example of typical settings for a frequency fit (with default seismic weights):
     # define_fit["freqparams"] = {
@@ -276,7 +281,15 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     # --> The names must match entries in the parameter list (basta/constants.py)
     # --> A reasonable choice is to (as a minimum) output the parameters used in the fit
     # --> If you want to predict distance, add the special keyword "distance".
-    define_output["outparams"] = ("Teff", "FeH", "radPhot", "massfin", "age")
+    define_output["outparams"] = (
+        "Teff",
+        "FeH",
+        "dnufit",
+        "numax",
+        "radPhot",
+        "massfin",
+        "age",
+    )
 
     # Name of the output file containing the results of the fit in ascii format.
     # --> A version in xml-format will be automatically created
diff --git a/examples/xmlinput/create_inputfile_interp_MS.py b/examples/xmlinput/create_inputfile_interp_MS.py
index 9bd3bed..0dd9ce7 100644
--- a/examples/xmlinput/create_inputfile_interp_MS.py
+++ b/examples/xmlinput/create_inputfile_interp_MS.py
@@ -29,7 +29,7 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     # The path to the grid to be used by BASTA for the fitting.
     # --> If using isochrones, remember to also specify physics settings in BLOCK 3c
     # --> If you need the location of BASTA, it is in BASTADIR
-    define_io["gridfile"] = os.path.join(BASTADIR, "grids", "Garstec_16CygA_v1.hdf5")
+    define_io["gridfile"] = os.path.join(BASTADIR, "grids", "Garstec_16CygA.hdf5")
 
     # Where to store the output of the BASTA run
     define_io["outputpath"] = os.path.join(BASTADIR, "examples", "output/interp_MS")
@@ -224,13 +224,17 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     #                    * "1/N-dof": Normalisation by number of frequencies minus the
     #                                 (user-specified) degrees of freedom. This option
     #                                 must be supplemented by a specification of "dof".
+    #
     # - "threepoint": Set 'True' to change to the three-point small frequency separation
     #                 formulation for the frequency ratios, instead of the five-point
     #                 small frequency separation. The default, if not set or set as
     #                 'False', is to use the five-point formulation.
+    #
+    # - "dnuprior": Set to 'False' to disable the automatic prior on dnu (mimicking the
+    #               range defined by dnufrac), which speeds up the computation of the
+    #               fit. Recommended to keep as 'True' (default).
 
     # Example of typical settings for a frequency fit (with default seismic weights):
-    # define_fit["freqparams"] = {
     define_fit["freqparams"] = {
         "freqpath": os.path.join(BASTADIR, "examples/data/freqs"),
         "fcor": "BG14",
@@ -277,7 +281,15 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     # --> The names must match entries in the parameter list (basta/constants.py)
     # --> A reasonable choice is to (as a minimum) output the parameters used in the fit
     # --> If you want to predict distance, add the special keyword "distance".
-    define_output["outparams"] = ("Teff", "FeH", "radPhot", "massfin", "age")
+    define_output["outparams"] = (
+        "Teff",
+        "FeH",
+        "dnufit",
+        "numax",
+        "radPhot",
+        "massfin",
+        "age",
+    )
 
     # Name of the output file containing the results of the fit in ascii format.
     # --> A version in xml-format will be automatically created
diff --git a/examples/xmlinput/create_inputfile_json.py b/examples/xmlinput/create_inputfile_json.py
index 85a412c..e4583b9 100644
--- a/examples/xmlinput/create_inputfile_json.py
+++ b/examples/xmlinput/create_inputfile_json.py
@@ -224,10 +224,15 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     #                    * "1/N-dof": Normalisation by number of frequencies minus the
     #                                 (user-specified) degrees of freedom. This option
     #                                 must be supplemented by a specification of "dof".
+    #
     # - "threepoint": Set 'True' to change to the three-point small frequency separation
     #                 formulation for the frequency ratios, instead of the five-point
     #                 small frequency separation. The default, if not set or set as
     #                 'False', is to use the five-point formulation.
+    #
+    # - "dnuprior": Set to 'False' to disable the automatic prior on dnu (mimicking the
+    #               range defined by dnufrac), which speeds up the computation of the
+    #               fit. Recommended to keep as 'True' (default).
 
     # Example of typical settings for a frequency fit (with default seismic weights):
     # define_fit["freqparams"] = {
@@ -276,7 +281,15 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     # --> The names must match entries in the parameter list (basta/constants.py)
     # --> A reasonable choice is to (as a minimum) output the parameters used in the fit
     # --> If you want to predict distance, add the special keyword "distance".
-    define_output["outparams"] = ("Teff", "FeH", "radPhot", "massfin", "age")
+    define_output["outparams"] = (
+        "Teff",
+        "FeH",
+        "dnufit",
+        "numax",
+        "radPhot",
+        "massfin",
+        "age",
+    )
 
     # Name of the output file containing the results of the fit in ascii format.
     # --> A version in xml-format will be automatically created
diff --git a/examples/xmlinput/create_inputfile_parallax.py b/examples/xmlinput/create_inputfile_parallax.py
index f1067bc..350eae5 100644
--- a/examples/xmlinput/create_inputfile_parallax.py
+++ b/examples/xmlinput/create_inputfile_parallax.py
@@ -231,14 +231,18 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     #                    * "1/N-dof": Normalisation by number of frequencies minus the
     #                                 (user-specified) degrees of freedom. This option
     #                                 must be supplemented by a specification of "dof".
+    #
     # - "threepoint": Set 'True' to change to the three-point small frequency separation
     #                 formulation for the frequency ratios, instead of the five-point
     #                 small frequency separation. The default, if not set or set as
     #                 'False', is to use the five-point formulation.
+    #
+    # - "dnuprior": Set to 'False' to disable the automatic prior on dnu (mimicking the
+    #               range defined by dnufrac), which speeds up the computation of the
+    #               fit. Recommended to keep as 'True' (default).
 
     # Example of typical settings for a frequency fit (with default seismic weights):
     # define_fit["freqparams"] = {
-    # define_fit["freqparams"] = {
     #     "freqpath": os.path.join(BASTADIR, "examples/data/freqs"),
     #     "fcor": "BG14",
     #     "correlations": False,
@@ -287,6 +291,8 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     define_output["outparams"] = (
         "Teff",
         "FeH",
+        "dnuSer",
+        "numax",
         "radPhot",
         "massfin",
         "age",
diff --git a/examples/xmlinput/create_inputfile_phase.py b/examples/xmlinput/create_inputfile_phase.py
index 2d4e0a2..473a582 100644
--- a/examples/xmlinput/create_inputfile_phase.py
+++ b/examples/xmlinput/create_inputfile_phase.py
@@ -233,14 +233,18 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     #                    * "1/N-dof": Normalisation by number of frequencies minus the
     #                                 (user-specified) degrees of freedom. This option
     #                                 must be supplemented by a specification of "dof".
+    #
     # - "threepoint": Set 'True' to change to the three-point small frequency separation
     #                 formulation for the frequency ratios, instead of the five-point
     #                 small frequency separation. The default, if not set or set as
     #                 'False', is to use the five-point formulation.
+    #
+    # - "dnuprior": Set to 'False' to disable the automatic prior on dnu (mimicking the
+    #               range defined by dnufrac), which speeds up the computation of the
+    #               fit. Recommended to keep as 'True' (default).
 
     # Example of typical settings for a frequency fit (with default seismic weights):
     # define_fit["freqparams"] = {
-    # define_fit["freqparams"] = {
     #     "freqpath": os.path.join(BASTADIR, "examples/data/freqs"),
     #     "fcor": "BG14",
     #     "correlations": False,
@@ -289,6 +293,8 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     define_output["outparams"] = (
         "Teff",
         "FeH",
+        "dnuSer",
+        "numax",
         "radPhot",
         "massfin",
         "age",
diff --git a/examples/xmlinput/create_inputfile_ratios.py b/examples/xmlinput/create_inputfile_ratios.py
index 51e95dc..a04f6ae 100644
--- a/examples/xmlinput/create_inputfile_ratios.py
+++ b/examples/xmlinput/create_inputfile_ratios.py
@@ -224,13 +224,17 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     #                    * "1/N-dof": Normalisation by number of frequencies minus the
     #                                 (user-specified) degrees of freedom. This option
     #                                 must be supplemented by a specification of "dof".
+    #
     # - "threepoint": Set 'True' to change to the three-point small frequency separation
     #                 formulation for the frequency ratios, instead of the five-point
     #                 small frequency separation. The default, if not set or set as
     #                 'False', is to use the five-point formulation.
+    #
+    # - "dnuprior": Set to 'False' to disable the automatic prior on dnu (mimicking the
+    #               range defined by dnufrac), which speeds up the computation of the
+    #               fit. Recommended to keep as 'True' (default).
 
     # Example of typical settings for a frequency fit (with default seismic weights):
-    # define_fit["freqparams"] = {
     define_fit["freqparams"] = {
         "freqpath": os.path.join(BASTADIR, "examples/data/freqs"),
         "fcor": "BG14",
@@ -277,7 +281,15 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     # --> The names must match entries in the parameter list (basta/constants.py)
     # --> A reasonable choice is to (as a minimum) output the parameters used in the fit
     # --> If you want to predict distance, add the special keyword "distance".
-    define_output["outparams"] = ("Teff", "FeH", "radPhot", "massfin", "age")
+    define_output["outparams"] = (
+        "Teff",
+        "FeH",
+        "dnufit",
+        "numax",
+        "radPhot",
+        "massfin",
+        "age",
+    )
 
     # Name of the output file containing the results of the fit in ascii format.
     # --> A version in xml-format will be automatically created
diff --git a/examples/xmlinput/create_inputfile_subgiant.py b/examples/xmlinput/create_inputfile_subgiant.py
index 08434fa..4802fb3 100644
--- a/examples/xmlinput/create_inputfile_subgiant.py
+++ b/examples/xmlinput/create_inputfile_subgiant.py
@@ -220,13 +220,17 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     #                    * "1/N-dof": Normalisation by number of frequencies minus the
     #                                 (user-specified) degrees of freedom. This option
     #                                 must be supplemented by a specification of "dof".
+    #
     # - "threepoint": Set 'True' to change to the three-point small frequency separation
     #                 formulation for the frequency ratios, instead of the five-point
     #                 small frequency separation. The default, if not set or set as
     #                 'False', is to use the five-point formulation.
+    #
+    # - "dnuprior": Set to 'False' to disable the automatic prior on dnu (mimicking the
+    #               range defined by dnufrac), which speeds up the computation of the
+    #               fit. Recommended to keep as 'True' (default).
 
     # Example of typical settings for a frequency fit (with default seismic weights):
-    # define_fit["freqparams"] = {
     define_fit["freqparams"] = {
         "freqpath": os.path.join(BASTADIR, "examples/data/freqs"),
         "fcor": "BG14",
@@ -273,7 +277,15 @@ def define_input(define_io, define_fit, define_output, define_plots, define_intp
     # --> The names must match entries in the parameter list (basta/constants.py)
     # --> A reasonable choice is to (as a minimum) output the parameters used in the fit
     # --> If you want to predict distance, add the special keyword "distance".
-    define_output["outparams"] = ("Teff", "FeH", "radPhot", "massfin", "age")
+    define_output["outparams"] = (
+        "Teff",
+        "FeH",
+        "dnufit",
+        "numax",
+        "radPhot",
+        "massfin",
+        "age",
+    )
 
     # Name of the output file containing the results of the fit in ascii format.
     # --> A version in xml-format will be automatically created
diff --git a/requirements.txt b/requirements.txt
index 1423149..3a6d094 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,12 +1,12 @@
-astropy==5.1.1
-black==22.10.0
-bottleneck==1.3.5
+astropy==5.2.1
+black==23.1.0
+bottleneck==1.3.6
 dustmaps==1.0.10
-h5py==3.7.0
-healpy==1.16.1
-matplotlib==3.6.2
-numpy==1.23.4
-pre-commit==2.20.0
-scikit-learn==1.1.3
-scipy==1.9.3
+h5py==3.8.0
+healpy==1.16.2
+matplotlib==3.7.0
+numpy==1.24.2
+pre-commit==3.1.1
+scikit-learn==1.2.1
+scipy==1.10.1
 tqdm==4.64.1