diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index b583fe99f65618ad741c23eba6dacc0ac2f6a6b2..e3843f1983639a2d060539b117f8ebc3c22756dd 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -79,6 +79,28 @@ python-3.7:
       - coverage_badge.svg
       - docs/_build/html/
 
+# Tests run at a fixed schedule rather than on push
+scheduled-python-3.7:
+  stage: test
+  image: bilbydev/bilby-test-suite-python37
+  only:
+    - schedules
+  before_script:
+    # Install the dependencies specified in the Pipfile
+    - pipenv install --three --python=/opt/conda/bin/python --system --deploy
+  script:
+    - python setup.py install
+
+    # Run pyflakes
+    - flake8 .
+
+    # Run tests
+    - pytest
+
+    # Run tests which are only done on schedule
+    - pytest test/example_test.py
+    - pytest test/gw_example_test.py
+
 pages:
   stage: deploy
   dependencies:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5d30ca96f4a6f2a22b87bf040cf9a6717c942016..e895e5f50f7d1d39311fccf115a6768fac30e4be 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -14,6 +14,7 @@
 ## [0.4.0] 2019-02-15
 
 ### Changed
+- Changed the logic around redundancy tests in the `PriorDict` classes
 - Fixed an accidental addition of astropy as a first-class dependency and added a check for missing dependencies to the C.I.
 - Fixed a bug in the "create-your-own-time-domain-model" example
 - Added citation guide to the readme
@@ -33,7 +34,7 @@
 - Improve the load_data_from_cache_file method
 
 ### Removed
--
+- Removed deprecated `PriorSet` classes. Use `PriorDict` instead.
 
 ## [0.3.5] 2019-01-25
 
diff --git a/bilby/core/prior.py b/bilby/core/prior.py
index 228148d6c11dae4671d6b424c09e7959a0fe469b..88ca370832769482cdab9881c0c0bbb862f2b708 100644
--- a/bilby/core/prior.py
+++ b/bilby/core/prior.py
@@ -245,10 +245,28 @@ class PriorDict(OrderedDict):
         """
         return [self[key].rescale(sample) for key, sample in zip(keys, theta)]
 
-    def test_redundancy(self, key):
+    def test_redundancy(self, key, disable_logging=False):
         """Empty redundancy test, should be overwritten in subclasses"""
         return False
 
+    def test_has_redundant_keys(self):
+        """
+        Test whether there are redundant keys in self.
+
+        Return
+        ------
+        bool: Whether there are redundancies or not
+        """
+        redundant = False
+        for key in self:
+            temp = self.copy()
+            del temp[key]
+            if temp.test_redundancy(key, disable_logging=True):
+                logger.warning('{} is a redundant key in this {}.'
+                               .format(key, self.__class__.__name__))
+                redundant = True
+        return redundant
+
     def copy(self):
         """
         We have to overwrite the copy method as it fails due to the presence of
diff --git a/bilby/core/result.py b/bilby/core/result.py
index 5c95f1e834fc7eb0b43f2c28d7a5b38335d1d725..3813e664eb753fad4516463720aa1c95beeb17db 100644
--- a/bilby/core/result.py
+++ b/bilby/core/result.py
@@ -8,6 +8,7 @@ import numpy as np
 import deepdish
 import pandas as pd
 import corner
+import json
 import scipy.stats
 import matplotlib
 import matplotlib.pyplot as plt
@@ -19,7 +20,7 @@ from .utils import (logger, infer_parameters_from_function,
 from .prior import Prior, PriorDict, DeltaFunction
 
 
-def result_file_name(outdir, label):
+def result_file_name(outdir, label, extension='json'):
     """ Returns the standard filename used for a result file
 
     Parameters
@@ -28,17 +29,27 @@ def result_file_name(outdir, label):
         Name of the output directory
     label: str
         Naming scheme of the output file
+    extension: str, optional
+        Whether to save as `hdf5` or `json`
 
     Returns
     -------
     str: File name of the output file
     """
-    return '{}/{}_result.h5'.format(outdir, label)
+    if extension == 'hdf5':
+        return '{}/{}_result.h5'.format(outdir, label)
+    else:
+        return '{}/{}_result.json'.format(outdir, label)
 
 
 def read_in_result(filename=None, outdir=None, label=None):
-    """ Wrapper to bilby.core.result.Result.from_hdf5 """
-    return Result.from_hdf5(filename=filename, outdir=outdir, label=label)
+    """ Wrapper to bilby.core.result.Result.from_hdf5
+        or bilby.core.result.Result.from_json """
+    try:
+        result = Result.from_json(filename=filename, outdir=outdir, label=label)
+    except (IOError, ValueError):
+        result = Result.from_hdf5(filename=filename, outdir=outdir, label=label)
+    return result
 
 
 class Result(object):
@@ -155,7 +166,7 @@ class Result(object):
             if (outdir is None) and (label is None):
                 raise ValueError("No information given to load file")
             else:
-                filename = result_file_name(outdir, label)
+                filename = result_file_name(outdir, label, extension='hdf5')
         if os.path.isfile(filename):
             dictionary = deepdish.io.load(filename)
             # Some versions of deepdish/pytables return the dictionanary as
@@ -169,6 +180,50 @@ class Result(object):
         else:
             raise IOError("No result '{}' found".format(filename))
 
+    @classmethod
+    def from_json(cls, filename=None, outdir=None, label=None):
+        """ Read in a saved .json data file
+
+        Parameters
+        ----------
+        filename: str
+            If given, try to load from this filename
+        outdir, label: str
+            If given, use the default naming convention for saved results file
+
+        Returns
+        -------
+        result: bilby.core.result.Result
+
+        Raises
+        -------
+        ValueError: If no filename is given and either outdir or label is None
+                    If no bilby.core.result.Result is found in the path
+
+        """
+        if filename is None:
+            if (outdir is None) and (label is None):
+                raise ValueError("No information given to load file")
+            else:
+                filename = result_file_name(outdir, label)
+        if os.path.isfile(filename):
+            dictionary = json.load(open(filename, 'r'))
+            for key in dictionary.keys():
+                # Convert some dictionaries back to DataFrames
+                if key in ['posterior', 'nested_samples']:
+                    dictionary[key] = pd.DataFrame.from_dict(dictionary[key])
+                # Convert the loaded priors to bilby prior type
+                if key == 'priors':
+                    for param in dictionary[key].keys():
+                        dictionary[key][param] = str(dictionary[key][param])
+                    dictionary[key] = PriorDict(dictionary[key])
+            try:
+                return cls(**dictionary)
+            except TypeError as e:
+                raise IOError("Unable to load dictionary, error={}".format(e))
+        else:
+            raise IOError("No result '{}' found".format(filename))
+
     def __str__(self):
         """Print a summary """
         if getattr(self, 'posterior', None) is not None:
@@ -303,9 +358,9 @@ class Result(object):
                 pass
         return dictionary
 
-    def save_to_file(self, overwrite=False, outdir=None):
+    def save_to_file(self, overwrite=False, outdir=None, extension='json'):
         """
-        Writes the Result to a deepdish h5 file
+        Writes the Result to a json or deepdish h5 file
 
         Parameters
         ----------
@@ -314,9 +369,11 @@ class Result(object):
             default=False
         outdir: str, optional
             Path to the outdir. Default is the one stored in the result object.
+        extension: str, optional
+            Whether to save as hdf5 instead of json
         """
         outdir = self._safe_outdir_creation(outdir, self.save_to_file)
-        file_name = result_file_name(outdir, self.label)
+        file_name = result_file_name(outdir, self.label, extension)
 
         if os.path.isfile(file_name):
             if overwrite:
@@ -341,8 +398,19 @@ class Result(object):
                 if hasattr(dictionary['sampler_kwargs'][key], '__call__'):
                     dictionary['sampler_kwargs'][key] = str(dictionary['sampler_kwargs'])
 
+        # Convert to json saveable format
+        if extension != 'hdf5':
+            for key in dictionary.keys():
+                if isinstance(dictionary[key], pd.core.frame.DataFrame):
+                    dictionary[key] = dictionary[key].to_dict()
+                elif isinstance(dictionary[key], np.ndarray):
+                    dictionary[key] = dictionary[key].tolist()
+
         try:
-            deepdish.io.save(file_name, dictionary)
+            if extension == 'hdf5':
+                deepdish.io.save(file_name, dictionary)
+            else:
+                json.dump(dictionary, open(file_name, 'w'), indent=2)
         except Exception as e:
             logger.error("\n\n Saving the data has failed with the "
                          "following message:\n {} \n\n".format(e))
diff --git a/bilby/core/sampler/__init__.py b/bilby/core/sampler/__init__.py
index 7a15b410cb1fb9ef4d425610a16bbd3c9d4ccc3b..147d95c5d43001dfb618dd77a8ced8a6e7c39258 100644
--- a/bilby/core/sampler/__init__.py
+++ b/bilby/core/sampler/__init__.py
@@ -85,6 +85,7 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
         overwritten.
     save: bool
         If true, save the priors and results to disk.
+        If hdf5, save as an hdf5 file instead of json.
     result_class: bilby.core.result.Result, or child of
         The result class to use. By default, `bilby.core.result.Result` is used,
         but objects which inherit from this class can be given providing
@@ -183,7 +184,10 @@ def run_sampler(likelihood, priors=None, label='label', outdir='outdir',
 
     result.samples_to_posterior(likelihood=likelihood, priors=priors,
                                 conversion_function=conversion_function)
-    if save:
+    if save == 'hdf5':
+        result.save_to_file(extension='hdf5')
+        logger.info("Results saved to {}/".format(outdir))
+    elif save:
         result.save_to_file()
         logger.info("Results saved to {}/".format(outdir))
     if plot:
diff --git a/bilby/core/sampler/base_sampler.py b/bilby/core/sampler/base_sampler.py
index 20f4366668dd29cc7b99a0ebc5c924cbe3ccd835..c07f208aec77fa1d27c41501d3c3833776f3d556 100644
--- a/bilby/core/sampler/base_sampler.py
+++ b/bilby/core/sampler/base_sampler.py
@@ -241,6 +241,10 @@ class Sampler(object):
             Likelihood can't be evaluated.
 
         """
+
+        if self.priors.test_has_redundant_keys():
+            raise IllegalSamplingSetError("Your sampling set contains redundant parameters.")
+
         self._check_if_priors_can_be_sampled()
         try:
             t1 = datetime.datetime.now()
@@ -502,3 +506,7 @@ class SamplerError(Error):
 
 class SamplerNotInstalledError(SamplerError):
     """ Base class for Error raised by not installed samplers """
+
+
+class IllegalSamplingSetError(Error):
+    """ Class for illegal sets of sampling parameters """
diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py
index ea5853d55eeac0440eae8cae539898fe9030ac10..eb55e3c9957a5b9649c07fa17f58e4b60fc1da47 100644
--- a/bilby/gw/conversion.py
+++ b/bilby/gw/conversion.py
@@ -215,7 +215,7 @@ def convert_to_lal_binary_black_hole_parameters(parameters):
             converted_parameters['phi_jl'] = 0.0
             converted_parameters['phi_12'] = 0.0
 
-    for angle in ['tilt_1', 'tilt_2', 'iota']:
+    for angle in ['tilt_1', 'tilt_2', 'theta_jn']:
         cos_angle = str('cos_' + angle)
         if cos_angle in converted_parameters.keys():
             converted_parameters[angle] =\
@@ -827,7 +827,7 @@ def generate_component_spins(sample):
     Parameters
     ----------
     sample: A dictionary with the necessary spin conversion parameters:
-    'iota', 'phi_jl', 'tilt_1', 'tilt_2', 'phi_12', 'a_1', 'a_2', 'mass_1',
+    'theta_jn', 'phi_jl', 'tilt_1', 'tilt_2', 'phi_12', 'a_1', 'a_2', 'mass_1',
     'mass_2', 'reference_frequency', 'phase'
 
     Returns
@@ -837,7 +837,7 @@ def generate_component_spins(sample):
     """
     output_sample = sample.copy()
     spin_conversion_parameters =\
-        ['iota', 'phi_jl', 'tilt_1', 'tilt_2', 'phi_12', 'a_1', 'a_2', 'mass_1',
+        ['theta_jn', 'phi_jl', 'tilt_1', 'tilt_2', 'phi_12', 'a_1', 'a_2', 'mass_1',
          'mass_2', 'reference_frequency', 'phase']
     if all(key in output_sample.keys() for key in spin_conversion_parameters):
         output_sample['iota'], output_sample['spin_1x'],\
@@ -845,7 +845,7 @@ def generate_component_spins(sample):
             output_sample['spin_2x'], output_sample['spin_2y'],\
             output_sample['spin_2z'] =\
             transform_precessing_spins(
-                output_sample['iota'], output_sample['phi_jl'],
+                output_sample['theta_jn'], output_sample['phi_jl'],
                 output_sample['tilt_1'], output_sample['tilt_2'],
                 output_sample['phi_12'], output_sample['a_1'],
                 output_sample['a_2'],
diff --git a/bilby/gw/detector.py b/bilby/gw/detector.py
index 9b01722dc445c8c0c1cea58b12ead3a864045cde..2d68559921f0c696e386c92be7a5b29af2ccf7d3 100644
--- a/bilby/gw/detector.py
+++ b/bilby/gw/detector.py
@@ -438,7 +438,9 @@ class InterferometerStrainData(object):
         strain = strain.filter(bp, filtfilt=True)
         self._time_domain_strain = strain.value
 
-    def create_power_spectral_density(self, fft_length, name='unknown', outdir=None):
+    def create_power_spectral_density(
+            self, fft_length, overlap=0, name='unknown', outdir=None,
+            analysis_segment_start_time=None):
         """ Use the time domain strain to generate a power spectral density
 
         This create a Tukey-windowed power spectral density and writes it to a
@@ -448,12 +450,17 @@ class InterferometerStrainData(object):
         ----------
         fft_length: float
             Duration of the analysis segment.
+        overlap: float
+            Number of seconds of overlap between FFTs.
         name: str
             The name of the detector, used in storing the PSD. Defaults to
             "unknown".
         outdir: str
             The output directory to write the PSD file too. If not given,
             the PSD will not be written to file.
+        analysis_segment_start_time: float
+            The start time of the analysis segment, if given, this data will
+            be removed before creating the PSD.
 
         Returns
         -------
@@ -461,11 +468,28 @@ class InterferometerStrainData(object):
             The frequencies and power spectral density array
 
         """
-        strain = gwpy.timeseries.TimeSeries(self.time_domain_strain, sample_rate=self.sampling_frequency)
+
+        data = self.time_domain_strain
+
+        if analysis_segment_start_time is not None:
+            analysis_segment_end_time = analysis_segment_start_time + fft_length
+            inside = (analysis_segment_start_time > self.time_array[0] +
+                      analysis_segment_end_time < self.time_array[-1])
+            if inside:
+                logger.info("Removing analysis segment data from the PSD data")
+                idxs = (
+                    (self.time_array < analysis_segment_start_time) +
+                    (self.time_array > analysis_segment_end_time))
+                data = data[idxs]
+
+        # WARNING this line can cause issues if the data is non-contiguous
+        strain = gwpy.timeseries.TimeSeries(data=data, sample_rate=self.sampling_frequency)
         psd_alpha = 2 * self.roll_off / fft_length
-        logger.info("Creating PSD with non-overlapping tukey window, alpha={}, roll off={}".format(
-            psd_alpha, self.roll_off))
-        psd = strain.psd(fftlength=fft_length, overlap=0, window=('tukey', psd_alpha))
+        logger.info(
+            "Tukey window PSD data with alpha={}, roll off={}".format(
+                psd_alpha, self.roll_off))
+        psd = strain.psd(
+            fftlength=fft_length, overlap=overlap, window=('tukey', psd_alpha))
 
         if outdir:
             psd_file = '{}/{}_PSD_{}_{}.txt'.format(outdir, name, self.start_time, self.duration)
@@ -1801,7 +1825,8 @@ class PowerSpectralDensity(object):
     @staticmethod
     def from_frame_file(frame_file, psd_start_time, psd_duration,
                         fft_length=4, sampling_frequency=4096, roll_off=0.2,
-                        channel=None):
+                        overlap=0, channel=None, name=None, outdir=None,
+                        analysis_segment_start_time=None):
         """ Generate power spectral density from a frame file
 
         Parameters
@@ -1819,15 +1844,25 @@ class PowerSpectralDensity(object):
             This is twice the maximum frequency.
         roll_off: float, optional
             Rise time in seconds of tukey window.
+        overlap: float,
+            Number of seconds of overlap between FFTs.
         channel: str, optional
             Name of channel to use to generate PSD.
+        name, outdir: str, optional
+            Name (and outdir) of the detector for which a PSD is to be
+            generated.
+        analysis_segment_start_time: float, optional
+            The start time of the analysis segment, if given, this data will
+            be removed before creating the PSD.
 
         """
         strain = InterferometerStrainData(roll_off=roll_off)
         strain.set_from_frame_file(
             frame_file, start_time=psd_start_time, duration=psd_duration,
             channel=channel, sampling_frequency=sampling_frequency)
-        frequency_array, psd_array = strain.create_power_spectral_density(fft_length=fft_length)
+        frequency_array, psd_array = strain.create_power_spectral_density(
+            fft_length=fft_length, name=name, outdir=outdir, overlap=overlap,
+            analysis_segment_start_time=analysis_segment_start_time)
         return PowerSpectralDensity(frequency_array=frequency_array, psd_array=psd_array)
 
     @staticmethod
@@ -2256,23 +2291,30 @@ def get_event_data(
 
 
 def load_data_from_cache_file(
-        cache_file, start_time, segment_duration, psd_duration,
-        channel_name=None, sampling_frequency=4096):
+        cache_file, start_time, segment_duration, psd_duration, psd_start_time,
+        channel_name=None, sampling_frequency=4096, roll_off=0.2,
+        overlap=0, outdir=None):
     """ Helper routine to generate an interferometer from a cache file
 
     Parameters
     ----------
     cache_file: str
         Path to the location of the cache file
-    start_time: float
-        GPS start time of the segment
+    start_time, psd_start_time: float
+        GPS start time of the segment and data stretch used for the PSD
     segment_duration, psd_duration: float
         Segment duration and duration of data to use to generate the PSD (in
         seconds).
+    roll_off: float, optional
+        Rise time in seconds of tukey window.
+    overlap: float,
+        Number of seconds of overlap between FFTs.
     channel_name: str
         Channel name
     sampling_frequency: int
         Sampling frequency
+    outdir: str, optional
+        The output directory in which the data is saved
 
     Returns
     -------
@@ -2283,7 +2325,6 @@ def load_data_from_cache_file(
 
     data_set = False
     psd_set = False
-    psd_start_time = start_time - psd_duration
 
     with open(cache_file, 'r') as ff:
         for line in ff:
@@ -2307,10 +2348,17 @@ def load_data_from_cache_file(
             if not psd_set & psd_in_cache:
                 ifo.power_spectral_density = \
                     PowerSpectralDensity.from_frame_file(
-                        cache.path, psd_start_time=psd_start_time,
+                        cache.path,
+                        psd_start_time=psd_start_time,
                         psd_duration=psd_duration,
+                        fft_length=segment_duration,
+                        sampling_frequency=sampling_frequency,
+                        roll_off=roll_off,
+                        overlap=overlap,
                         channel=channel_name,
-                        sampling_frequency=sampling_frequency)
+                        name=cache.observatory,
+                        outdir=outdir,
+                        analysis_segment_start_time=start_time)
                 psd_set = True
     if data_set and psd_set:
         return ifo
diff --git a/bilby/gw/likelihood.py b/bilby/gw/likelihood.py
index 1758efbd3978f907367e823878f3d09a320910ab..43d7fcf7a83cdcd94158dbf0f246cea8a6fbf5ec 100644
--- a/bilby/gw/likelihood.py
+++ b/bilby/gw/likelihood.py
@@ -412,10 +412,13 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
 
     """
     def __init__(self, interferometers, waveform_generator,
-                 linear_matrix, quadratic_matrix, priors):
+                 linear_matrix, quadratic_matrix, priors,
+                 distance_marginalization=False, phase_marginalization=False):
         GravitationalWaveTransient.__init__(
             self, interferometers=interferometers,
-            waveform_generator=waveform_generator, priors=priors)
+            waveform_generator=waveform_generator, priors=priors,
+            distance_marginalization=distance_marginalization,
+            phase_marginalization=phase_marginalization)
 
         if isinstance(linear_matrix, str):
             logger.info("Loading linear matrix from {}".format(linear_matrix))
@@ -434,7 +437,7 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
 
     def log_likelihood_ratio(self):
         optimal_snr_squared = 0.
-        matched_filter_snr_squared = 0.
+        d_inner_h = 0.
 
         indices, in_bounds = self._closest_time_indices(
             self.parameters['geocent_time'] - self.interferometers.start_time)
@@ -470,19 +473,27 @@ class ROQGravitationalWaveTransient(GravitationalWaveTransient):
             if not in_bounds:
                 return np.nan_to_num(-np.inf)
 
-            matched_filter_snr_squared_array = np.einsum(
+            d_inner_h_tc_array = np.einsum(
                 'i,ji->j', np.conjugate(h_plus_linear + h_cross_linear),
                 self.weights[ifo.name + '_linear'][indices])
 
-            matched_filter_snr_squared += interp1d(
+            d_inner_h += interp1d(
                 self.time_samples[indices],
-                matched_filter_snr_squared_array, kind='quadratic')(ifo_time)
+                d_inner_h_tc_array, kind='quadratic')(ifo_time)
 
             optimal_snr_squared += \
                 np.vdot(np.abs(h_plus_quadratic + h_cross_quadratic)**2,
                         self.weights[ifo.name + '_quadratic'])
 
-        log_l = matched_filter_snr_squared - optimal_snr_squared / 2
+        if self.distance_marginalization:
+            rho_mf_ref, rho_opt_ref = self._setup_rho(d_inner_h, optimal_snr_squared)
+            if self.phase_marginalization:
+                rho_mf_ref = abs(rho_mf_ref)
+            log_l = self._interp_dist_margd_loglikelihood(rho_mf_ref.real, rho_opt_ref)[0]
+        else:
+            if self.phase_marginalization:
+                d_inner_h = self._bessel_function_interped(abs(d_inner_h))
+            log_l = d_inner_h - optimal_snr_squared / 2
 
         return log_l.real
 
diff --git a/bilby/gw/prior.py b/bilby/gw/prior.py
index 68f8b34508985dbf2c3cc8a531c7514e5c8cdb42..99c353b44a2dd4ebc4b73669a483b2a2438b4e88 100644
--- a/bilby/gw/prior.py
+++ b/bilby/gw/prior.py
@@ -207,63 +207,48 @@ class BBHPriorDict(PriorDict):
                 filename = os.path.join(os.path.dirname(__file__), 'prior_files', filename)
         PriorDict.__init__(self, dictionary=dictionary, filename=filename)
 
-    def test_redundancy(self, key):
+    def test_redundancy(self, key, disable_logging=False):
         """
         Test whether adding the key would add be redundant.
+        Already existing keys return True.
 
         Parameters
         ----------
         key: str
             The key to test.
+        disable_logging: bool, optional
+            Disable logging in this function call. Default is False.
 
         Return
         ------
         redundant: bool
             Whether the key is redundant or not
         """
-        redundant = False
         if key in self:
             logger.debug('{} already in prior'.format(key))
-            return redundant
+            return True
+
         mass_parameters = {'mass_1', 'mass_2', 'chirp_mass', 'total_mass', 'mass_ratio', 'symmetric_mass_ratio'}
-        spin_magnitude_parameters = {'a_1', 'a_2'}
         spin_tilt_1_parameters = {'tilt_1', 'cos_tilt_1'}
         spin_tilt_2_parameters = {'tilt_2', 'cos_tilt_2'}
         spin_azimuth_parameters = {'phi_1', 'phi_2', 'phi_12', 'phi_jl'}
-        inclination_parameters = {'iota', 'cos_iota'}
+        inclination_parameters = {'theta_jn', 'cos_theta_jn'}
         distance_parameters = {'luminosity_distance', 'comoving_distance', 'redshift'}
 
-        for parameter_set in [mass_parameters, spin_magnitude_parameters, spin_azimuth_parameters]:
-            if key in parameter_set:
-                if len(parameter_set.intersection(self)) > 2:
-                    redundant = True
-                    logger.warning('{} in prior. This may lead to unexpected behaviour.'.format(
-                        parameter_set.intersection(self)))
-                    break
-            elif len(parameter_set.intersection(self)) == 2:
-                redundant = True
-                break
-        for parameter_set in [inclination_parameters, distance_parameters, spin_tilt_1_parameters,
-                              spin_tilt_2_parameters]:
+        for independent_parameters, parameter_set in \
+                zip([2, 2, 1, 1, 1, 1],
+                    [mass_parameters, spin_azimuth_parameters,
+                     spin_tilt_1_parameters, spin_tilt_2_parameters,
+                     inclination_parameters, distance_parameters]):
             if key in parameter_set:
-                if len(parameter_set.intersection(self)) > 1:
-                    redundant = True
-                    logger.warning('{} in prior. This may lead to unexpected behaviour.'.format(
-                        parameter_set.intersection(self)))
-                    break
-                elif len(parameter_set.intersection(self)) == 1:
-                    redundant = True
-                    break
-
-        return redundant
-
-
-class BBHPriorSet(BBHPriorDict):
-
-    def __init__(self, dictionary=None, filename=None):
-        """ DEPRECATED: USE BBHPriorDict INSTEAD"""
-        logger.warning("The name 'BBHPriorSet' is deprecated use 'BBHPriorDict' instead")
-        super(BBHPriorSet, self).__init__(dictionary, filename)
+                if len(parameter_set.intersection(self)) >= independent_parameters:
+                    logger.disabled = disable_logging
+                    logger.warning('{} already in prior. '
+                                   'This may lead to unexpected behaviour.'
+                                   .format(parameter_set.intersection(self)))
+                    logger.disabled = False
+                    return True
+        return False
 
 
 class BNSPriorDict(PriorDict):
@@ -286,34 +271,32 @@ class BNSPriorDict(PriorDict):
                 filename = os.path.join(os.path.dirname(__file__), 'prior_files', filename)
         PriorDict.__init__(self, dictionary=dictionary, filename=filename)
 
-    def test_redundancy(self, key):
-        logger.info("Performing redundancy check using BBHPriorDict().test_redundancy")
-        bbh_redundancy = BBHPriorDict().test_redundancy(key)
+    def test_redundancy(self, key, disable_logging=False):
+        logger.disabled = disable_logging
+        logger.info("Performing redundancy check using BBHPriorDict(self).test_redundancy")
+        logger.disabled = False
+        bbh_redundancy = BBHPriorDict(self).test_redundancy(key)
+
         if bbh_redundancy:
             return True
         redundant = False
 
-        tidal_parameters =\
+        tidal_parameters = \
             {'lambda_1', 'lambda_2', 'lambda_tilde', 'delta_lambda'}
 
         if key in tidal_parameters:
             if len(tidal_parameters.intersection(self)) > 2:
                 redundant = True
-                logger.warning('{} in prior. This may lead to unexpected behaviour.'.format(
-                    tidal_parameters.intersection(self)))
+                logger.disabled = disable_logging
+                logger.warning('{} already in prior. '
+                               'This may lead to unexpected behaviour.'
+                               .format(tidal_parameters.intersection(self)))
+                logger.disabled = False
             elif len(tidal_parameters.intersection(self)) == 2:
                 redundant = True
         return redundant
 
 
-class BNSPriorSet(BNSPriorDict):
-
-    def __init__(self, dictionary=None, filename=None):
-        """ DEPRECATED: USE BNSPriorDict INSTEAD"""
-        super(BNSPriorSet, self).__init__(dictionary, filename)
-        logger.warning("The name 'BNSPriorSet' is deprecated use 'BNSPriorDict' instead")
-
-
 Prior._default_latex_labels = {
     'mass_1': '$m_1$',
     'mass_2': '$m_2$',
@@ -334,6 +317,8 @@ Prior._default_latex_labels = {
     'ra': '$\mathrm{RA}$',
     'iota': '$\iota$',
     'cos_iota': '$\cos\iota$',
+    'theta_jn': '$\\theta_{JN}$',
+    'cos_theta_jn': '$\cos\\theta_{JN}$',
     'psi': '$\psi$',
     'phase': '$\phi$',
     'geocent_time': '$t_c$',
@@ -418,13 +403,13 @@ class CalibrationPriorDict(PriorDict):
         nodes = np.logspace(np.log10(minimum_frequency),
                             np.log10(maximum_frequency), n_nodes)
 
-        amplitude_mean_nodes =\
+        amplitude_mean_nodes = \
             UnivariateSpline(frequency_array, amplitude_median)(nodes)
-        amplitude_sigma_nodes =\
+        amplitude_sigma_nodes = \
             UnivariateSpline(frequency_array, amplitude_sigma)(nodes)
-        phase_mean_nodes =\
+        phase_mean_nodes = \
             UnivariateSpline(frequency_array, phase_median)(nodes)
-        phase_sigma_nodes =\
+        phase_sigma_nodes = \
             UnivariateSpline(frequency_array, phase_sigma)(nodes)
 
         prior = CalibrationPriorDict()
@@ -506,11 +491,3 @@ class CalibrationPriorDict(PriorDict):
                                         latex_label=latex_label)
 
         return prior
-
-
-class CalibrationPriorSet(CalibrationPriorDict):
-
-    def __init__(self, dictionary=None, filename=None):
-        """ DEPRECATED: USE BNSPriorDict INSTEAD"""
-        super(CalibrationPriorSet, self).__init__(dictionary, filename)
-        logger.warning("The name 'CalibrationPriorSet' is deprecated use 'CalibrationPriorDict' instead")
diff --git a/bilby/gw/prior_files/GW150914.prior b/bilby/gw/prior_files/GW150914.prior
index 35b12524c498e5a21127e91c74472c9dffe48430..0eef13a9734f8f0e37dc6c578e3d3a0721e28dad 100644
--- a/bilby/gw/prior_files/GW150914.prior
+++ b/bilby/gw/prior_files/GW150914.prior
@@ -10,7 +10,7 @@ phi_jl =  Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi)
 luminosity_distance =  bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=1e3, unit='Mpc')
 dec =  Cosine(name='dec')
 ra =  Uniform(name='ra', minimum=0, maximum=2 * np.pi)
-iota =  Sine(name='iota')
+theta_jn =  Sine(name='theta_jn')
 psi =  Uniform(name='psi', minimum=0, maximum=np.pi)
 phase =  Uniform(name='phase', minimum=0, maximum=2 * np.pi)
 geocent_time = Uniform(1126259462.322, 1126259462.522, name='geocent_time', unit='$s$')
diff --git a/bilby/gw/prior_files/binary_black_holes.prior b/bilby/gw/prior_files/binary_black_holes.prior
index 5314635bd779f09098e3c8c2002f3367f046547d..e41cc73e2f765fbddd821693604ac15e302586f2 100644
--- a/bilby/gw/prior_files/binary_black_holes.prior
+++ b/bilby/gw/prior_files/binary_black_holes.prior
@@ -19,7 +19,7 @@ phi_jl =  Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi)
 luminosity_distance =  bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=5e3, unit='Mpc')
 dec =  Cosine(name='dec')
 ra =  Uniform(name='ra', minimum=0, maximum=2 * np.pi)
-iota =  Sine(name='iota')
-# cos_iota =  Uniform(name='cos_iota', minimum=-1, maximum=1)
+theta_jn =  Sine(name='theta_jn')
+# cos_theta_jn =  Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
 psi =  Uniform(name='psi', minimum=0, maximum=np.pi)
 phase =  Uniform(name='phase', minimum=0, maximum=2 * np.pi)
diff --git a/bilby/gw/prior_files/binary_neutron_stars.prior b/bilby/gw/prior_files/binary_neutron_stars.prior
index 1bc4d485c1cfdf89ddee3ab4b7b33d6cabd80654..08ec88f2ffb8bd6536563cb20723aabb589c7c13 100644
--- a/bilby/gw/prior_files/binary_neutron_stars.prior
+++ b/bilby/gw/prior_files/binary_neutron_stars.prior
@@ -13,8 +13,8 @@ chi_2 =  bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1
 luminosity_distance =  bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=10, maximum=500, unit='Mpc')
 dec =  Cosine(name='dec')
 ra =  Uniform(name='ra', minimum=0, maximum=2 * np.pi)
-iota =  Sine(name='iota')
-# cos_iota =  Uniform(name='cos_iota', minimum=-1, maximum=1)
+theta_jn =  Sine(name='theta_jn')
+# cos_theta_jn =  Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
 psi =  Uniform(name='psi', minimum=0, maximum=np.pi)
 phase =  Uniform(name='phase', minimum=0, maximum=2 * np.pi)
 lambda_1 = Uniform(name='lambda_1', minimum=0, maximum=3000 )
diff --git a/bilby/gw/source.py b/bilby/gw/source.py
index a3c68a544449ae96108faa5df2a968066062b131..cfc50597c8ea56fdda7c41acc6118af0774b16f0 100644
--- a/bilby/gw/source.py
+++ b/bilby/gw/source.py
@@ -23,7 +23,7 @@ except ImportError:
 
 def lal_binary_black_hole(
         frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
-        phi_12, a_2, tilt_2, phi_jl, iota, phase, **kwargs):
+        phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **kwargs):
     """ A Binary Black Hole waveform model using lalsimulation
 
     Parameters
@@ -41,15 +41,16 @@ def lal_binary_black_hole(
     tilt_1: float
         Primary tilt angle
     phi_12: float
-
+        Azimuthal angle between the two component spins
     a_2: float
         Dimensionless secondary spin magnitude
     tilt_2: float
         Secondary tilt angle
     phi_jl: float
-
-    iota: float
-        Orbital inclination
+        Azimuthal angle between the total binary angular momentum and the
+        orbital angular momentum
+    theta_jn: float
+        Angle between the total binary angular momentum and the line of sight
     phase: float
         The phase at coalescence
     kwargs: dict
@@ -65,14 +66,14 @@ def lal_binary_black_hole(
     waveform_kwargs.update(kwargs)
     return _base_lal_cbc_fd_waveform(
         frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
-        luminosity_distance=luminosity_distance, iota=iota, phase=phase,
+        luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
         a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_12=phi_12,
         phi_jl=phi_jl, **waveform_kwargs)
 
 
 def lal_binary_neutron_star(
         frequency_array, mass_1, mass_2, luminosity_distance, chi_1, chi_2,
-        iota, phase, lambda_1, lambda_2, **kwargs):
+        theta_jn, phase, lambda_1, lambda_2, **kwargs):
     """ A Binary Neutron Star waveform model using lalsimulation
 
     Parameters
@@ -89,7 +90,7 @@ def lal_binary_neutron_star(
         Dimensionless aligned spin
     chi_2: float
         Dimensionless aligned spin
-    iota: float
+    theta_jn: float
         Orbital inclination
     phase: float
         The phase at coalescence
@@ -123,16 +124,15 @@ def lal_binary_neutron_star(
     waveform_kwargs.update(kwargs)
     return _base_lal_cbc_fd_waveform(
         frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
-        luminosity_distance=luminosity_distance, iota=iota, phase=phase,
+        luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
         a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, lambda_1=lambda_1,
         lambda_2=lambda_2, **waveform_kwargs)
 
 
 def lal_eccentric_binary_black_hole_no_spins(
         frequency_array, mass_1, mass_2, eccentricity, luminosity_distance,
-        iota, phase, **kwargs):
-    """
-    Eccentric binary black hole waveform model using lalsimulation (EccentricFD)
+        theta_jn, phase, **kwargs):
+    """ Eccentric binary black hole waveform model using lalsimulation (EccentricFD)
 
     Parameters
     ----------
@@ -146,7 +146,7 @@ def lal_eccentric_binary_black_hole_no_spins(
         The orbital eccentricity of the system
     luminosity_distance: float
         The luminosity distance in megaparsec
-    iota: float
+    theta_jn: float
         Orbital inclination
     phase: float
         The phase at coalescence
@@ -163,12 +163,12 @@ def lal_eccentric_binary_black_hole_no_spins(
     waveform_kwargs.update(kwargs)
     return _base_lal_cbc_fd_waveform(
         frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
-        luminosity_distance=luminosity_distance, iota=iota, phase=phase,
+        luminosity_distance=luminosity_distance, theta_jn=theta_jn, phase=phase,
         eccentricity=eccentricity, **waveform_kwargs)
 
 
 def _base_lal_cbc_fd_waveform(
-        frequency_array, mass_1, mass_2, luminosity_distance, iota, phase,
+        frequency_array, mass_1, mass_2, luminosity_distance, theta_jn, phase,
         a_1=0.0, a_2=0.0, tilt_1=0.0, tilt_2=0.0, phi_12=0.0, phi_jl=0.0,
         lambda_1=0.0, lambda_2=0.0, eccentricity=0.0, **waveform_kwargs):
     """ Generate a cbc waveform model using lalsimulation
@@ -239,7 +239,7 @@ def _base_lal_cbc_fd_waveform(
     else:
         iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = (
             transform_precessing_spins(
-                iota, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1,
+                theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1,
                 mass_2, reference_frequency, phase))
 
     longitude_ascending_nodes = 0.0
@@ -328,7 +328,7 @@ def supernova_pca_model(
 
 
 def roq(frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
-        phi_12, a_2, tilt_2, phi_jl, iota, phase, **waveform_arguments):
+        phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **waveform_arguments):
     """
     See https://git.ligo.org/lscsoft/lalsuite/blob/master/lalsimulation/src/LALSimInspiral.c#L1460
 
@@ -354,7 +354,7 @@ def roq(frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
         Secondary tilt angle
     phi_jl: float
 
-    iota: float
+    theta_jn: float
         Orbital inclination
     phase: float
         The phase at coalescence
@@ -399,11 +399,12 @@ def roq(frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
         spin_2x = 0
         spin_2y = 0
         spin_2z = a_2
+        iota = theta_jn
     else:
         iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = \
-            transform_precessing_spins(
-                iota, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1, mass_2,
-                reference_frequency, phase)
+            lalsim_SimInspiralTransformPrecessingNewInitialConditions(
+                theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1,
+                mass_2, reference_frequency, phase)
 
     chi_1_l, chi_2_l, chi_p, theta_jn, alpha, phase_aligned, zeta =\
         lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame(
diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py
index 94c3cafb2c0cdb63a73e72f6cef8db7c9757974a..82d6c57b58fd29ef27e800909caf4e67b8e6f448 100644
--- a/bilby/gw/utils.py
+++ b/bilby/gw/utils.py
@@ -130,6 +130,7 @@ def get_polarization_tensor(ra, dec, time, psi, mode):
     elif mode.lower() == 'breathing':
         return np.einsum('i,j->ij', m, m) + np.einsum('i,j->ij', n, n)
 
+    # Calculating omega here to avoid calculation when model in [plus, cross, breathing]
     omega = np.cross(m, n)
     if mode.lower() == 'longitudinal':
         return np.sqrt(2) * np.einsum('i,j->ij', omega, omega)
@@ -138,8 +139,7 @@ def get_polarization_tensor(ra, dec, time, psi, mode):
     elif mode.lower() == 'y':
         return np.einsum('i,j->ij', n, omega) + np.einsum('i,j->ij', omega, n)
     else:
-        logger.warning("{} not a polarization mode!".format(mode))
-        return None
+        raise ValueError("{} not a polarization mode!".format(mode))
 
 
 def get_vertex_position_geocentric(latitude, longitude, elevation):
@@ -380,7 +380,7 @@ def get_open_strain_data(
     return strain
 
 
-def read_frame_file(file_name, start_time, end_time, channel=None, buffer_time=1, **kwargs):
+def read_frame_file(file_name, start_time, end_time, channel=None, buffer_time=0, **kwargs):
     """ A function which accesses the open strain data
 
     This uses `gwpy` to download the open data and then saves a cached copy for
@@ -416,23 +416,22 @@ def read_frame_file(file_name, start_time, end_time, channel=None, buffer_time=1
         except RuntimeError:
             logger.warning('Channel {} not found. Trying preset channel names'.format(channel))
 
-    while not loaded:
-        ligo_channel_types = ['GDS-CALIB_STRAIN', 'DCS-CALIB_STRAIN_C01', 'DCS-CALIB_STRAIN_C02',
-                              'DCH-CLEAN_STRAIN_C02']
-        virgo_channel_types = ['Hrec_hoft_V1O2Repro2A_16384Hz', 'FAKE_h_16384Hz_4R']
-        channel_types = dict(H1=ligo_channel_types, L1=ligo_channel_types, V1=virgo_channel_types)
-        for detector in channel_types.keys():
-            for channel_type in channel_types[detector]:
-                if loaded:
-                    break
-                channel = '{}:{}'.format(detector, channel_type)
-                try:
-                    strain = TimeSeries.read(source=file_name, channel=channel, start=start_time, end=end_time,
-                                             **kwargs)
-                    loaded = True
-                    logger.info('Successfully read strain data for channel {}.'.format(channel))
-                except RuntimeError:
-                    pass
+    ligo_channel_types = ['GDS-CALIB_STRAIN', 'DCS-CALIB_STRAIN_C01', 'DCS-CALIB_STRAIN_C02',
+                          'DCH-CLEAN_STRAIN_C02']
+    virgo_channel_types = ['Hrec_hoft_V1O2Repro2A_16384Hz', 'FAKE_h_16384Hz_4R']
+    channel_types = dict(H1=ligo_channel_types, L1=ligo_channel_types, V1=virgo_channel_types)
+    for detector in channel_types.keys():
+        for channel_type in channel_types[detector]:
+            if loaded:
+                break
+            channel = '{}:{}'.format(detector, channel_type)
+            try:
+                strain = TimeSeries.read(source=file_name, channel=channel, start=start_time, end=end_time,
+                                         **kwargs)
+                loaded = True
+                logger.info('Successfully read strain data for channel {}.'.format(channel))
+            except RuntimeError:
+                pass
 
     if loaded:
         return strain
diff --git a/examples/injection_examples/australian_detector.py b/examples/injection_examples/australian_detector.py
index d770af182ce106e4b0e760a78413cb624bdf9aa7..3cd8c4aa94dd039e77644f9e69680892c3b6a147 100644
--- a/examples/injection_examples/australian_detector.py
+++ b/examples/injection_examples/australian_detector.py
@@ -59,7 +59,7 @@ interferometers.append(AusIFO)
 # signal at 4 Gpc
 injection_parameters = dict(
     mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=0.2108)
 
 
diff --git a/examples/injection_examples/binary_neutron_star_example.py b/examples/injection_examples/binary_neutron_star_example.py
index f1ad6d2d372ec18049630a2b0d92f199c9fa7ad1..cd540368845817ff1b7162d0cf63f5ba427fddf1 100644
--- a/examples/injection_examples/binary_neutron_star_example.py
+++ b/examples/injection_examples/binary_neutron_star_example.py
@@ -28,7 +28,7 @@ np.random.seed(88170235)
 # aligned spins of both black holes (chi_1, chi_2), etc.
 injection_parameters = dict(
     mass_1=1.5, mass_2=1.3, chi_1=0.02, chi_2=0.02, luminosity_distance=50.,
-    iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
+    theta_jn=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
     ra=1.375, dec=-1.2108, lambda_1=400, lambda_2=450)
 
 # Set the duration and sampling frequency of the data segment that we're going
@@ -66,7 +66,7 @@ interferometers.inject_signal(parameters=injection_parameters,
 # delta_lambda rather than mass_1, mass_2, lambda_1, and lambda_2.
 priors = bilby.gw.prior.BNSPriorDict()
 for key in ['psi', 'geocent_time', 'ra', 'dec', 'chi_1', 'chi_2',
-            'iota', 'luminosity_distance', 'phase']:
+            'theta_jn', 'luminosity_distance', 'phase']:
     priors[key] = injection_parameters[key]
 priors.pop('mass_1')
 priors.pop('mass_2')
diff --git a/examples/injection_examples/calibration_example.py b/examples/injection_examples/calibration_example.py
index ec3f3491e9ad61fda38512557d91a5345cf68f2d..ca58994696f4233ceb60d17bf500e144671e6c72 100644
--- a/examples/injection_examples/calibration_example.py
+++ b/examples/injection_examples/calibration_example.py
@@ -28,7 +28,7 @@ np.random.seed(88170235)
 # spins of both black holes (a, tilt, phi), etc.
 injection_parameters = dict(
     mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
 
 # Fixed arguments passed into the source model
diff --git a/examples/injection_examples/change_sampled_parameters.py b/examples/injection_examples/change_sampled_parameters.py
index b2c848576bcde36985796e28c206c887411114d3..9f432ef52910fd125b1930ea90420e7eae0b6245 100644
--- a/examples/injection_examples/change_sampled_parameters.py
+++ b/examples/injection_examples/change_sampled_parameters.py
@@ -22,7 +22,7 @@ np.random.seed(151226)
 
 injection_parameters = dict(
     total_mass=66., mass_ratio=0.9, a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=2000, iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=2000, theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
 
 waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
@@ -63,8 +63,8 @@ priors['redshift'] = bilby.prior.Uniform(
 for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'psi',
             'ra', 'dec', 'geocent_time', 'phase']:
     priors[key] = injection_parameters[key]
-priors.pop('iota')
-priors['cos_iota'] = np.cos(injection_parameters['iota'])
+priors.pop('theta_jn')
+priors['cos_theta_jn'] = np.cos(injection_parameters['theta_jn'])
 print(priors)
 
 # Initialise GravitationalWaveTransient
diff --git a/examples/injection_examples/eccentric_inspiral.py b/examples/injection_examples/eccentric_inspiral.py
index ce4635c7b982cc4a44e0f7d976c1b300d0e1bc19..8d3053f353fb769384f25a17c3b1d406860fd88e 100644
--- a/examples/injection_examples/eccentric_inspiral.py
+++ b/examples/injection_examples/eccentric_inspiral.py
@@ -27,7 +27,7 @@ np.random.seed(150914)
 
 injection_parameters = dict(
     mass_1=35., mass_2=30., eccentricity=0.1, luminosity_distance=440.,
-    iota=0.4, psi=0.1, phase=1.2, geocent_time=1180002601.0, ra=45, dec=5.73)
+    theta_jn=0.4, psi=0.1, phase=1.2, geocent_time=1180002601.0, ra=45, dec=5.73)
 
 waveform_arguments = dict(waveform_approximant='EccentricFD',
                           reference_frequency=10., minimum_frequency=10.)
@@ -70,7 +70,7 @@ priors["luminosity_distance"] = bilby.gw.prior.UniformComovingVolume(
 priors["dec"] = bilby.core.prior.Cosine(name='dec')
 priors["ra"] = bilby.core.prior.Uniform(
     name='ra', minimum=0, maximum=2 * np.pi)
-priors["iota"] = bilby.core.prior.Sine(name='iota')
+priors["theta_jn"] = bilby.core.prior.Sine(name='theta_jn')
 priors["psi"] = bilby.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi)
 priors["phase"] = bilby.core.prior.Uniform(
     name='phase', minimum=0, maximum=2 * np.pi)
diff --git a/examples/injection_examples/fast_tutorial.py b/examples/injection_examples/fast_tutorial.py
index c472ba44b62af83b0858ae407bc4653cbb2341f2..6bfd0a0ed8bab100ee78d0fd03476a40b83b04fb 100644
--- a/examples/injection_examples/fast_tutorial.py
+++ b/examples/injection_examples/fast_tutorial.py
@@ -31,7 +31,7 @@ np.random.seed(88170235)
 # spins of both black holes (a, tilt, phi), etc.
 injection_parameters = dict(
     mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
 
 # Fixed arguments passed into the source model
@@ -61,7 +61,7 @@ ifos.inject_signal(waveform_generator=waveform_generator,
 # prior is a delta function at the true, injected value.  In reality, the
 # sampler implementation is smart enough to not sample any parameter that has
 # a delta-function prior.
-# The above list does *not* include mass_1, mass_2, iota and luminosity
+# The above list does *not* include mass_1, mass_2, theta_jn and luminosity
 # distance, which means those are the parameters that will be included in the
 # sampler.  If we do nothing, then the default priors get used.
 priors = bilby.gw.prior.BBHPriorDict()
diff --git a/examples/injection_examples/how_to_specify_the_prior.py b/examples/injection_examples/how_to_specify_the_prior.py
index 885ab328b955bc5794cf63bc194dc428972ac1fd..68081b4d2eddabec080e48005c6fc5536bb546e1 100644
--- a/examples/injection_examples/how_to_specify_the_prior.py
+++ b/examples/injection_examples/how_to_specify_the_prior.py
@@ -17,7 +17,7 @@ np.random.seed(151012)
 
 injection_parameters = dict(
     mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
 
 waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
@@ -41,7 +41,7 @@ ifos.inject_signal(waveform_generator=waveform_generator,
 # This loads in a predefined set of priors for BBHs.
 priors = bilby.gw.prior.BBHPriorDict()
 # These parameters will not be sampled
-for key in ['tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'phase', 'iota', 'ra',
+for key in ['tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'phase', 'theta_jn', 'ra',
             'dec', 'geocent_time', 'psi']:
     priors[key] = injection_parameters[key]
 # We can make uniform distributions.
diff --git a/examples/injection_examples/marginalized_likelihood.py b/examples/injection_examples/marginalized_likelihood.py
index 06ef88b6731622146e984d5a750861ac4967bd56..6d84b1893c6b1f08419172c5e3ff77a9f5eba915 100644
--- a/examples/injection_examples/marginalized_likelihood.py
+++ b/examples/injection_examples/marginalized_likelihood.py
@@ -17,7 +17,7 @@ np.random.seed(170608)
 
 injection_parameters = dict(
     mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
 
 waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
@@ -40,7 +40,7 @@ ifos.inject_signal(waveform_generator=waveform_generator,
 # Set up prior
 priors = bilby.gw.prior.BBHPriorDict()
 # These parameters will not be sampled
-for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'iota', 'ra',
+for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'theta_jn', 'ra',
             'dec']:
     priors[key] = injection_parameters[key]
 
diff --git a/examples/injection_examples/plot_skymap.py b/examples/injection_examples/plot_skymap.py
index 6453669ba0bc52396e892f1312351d69562e201a..9ebee7725fb0b913a2329c56622367f51033e129 100644
--- a/examples/injection_examples/plot_skymap.py
+++ b/examples/injection_examples/plot_skymap.py
@@ -12,7 +12,7 @@ outdir = 'outdir'
 label = 'plot_skymap'
 injection_parameters = dict(
     mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-0.2108)
 
 waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
@@ -33,7 +33,7 @@ ifos.inject_signal(waveform_generator=waveform_generator,
 priors = bilby.gw.prior.BBHPriorDict()
 for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'phi_12', 'phi_jl', 'psi',
             'mass_1', 'mass_2', 'phase', 'geocent_time', 'luminosity_distance',
-            'iota']:
+            'theta_jn']:
     priors[key] = injection_parameters[key]
 
 likelihood = bilby.gw.GravitationalWaveTransient(
diff --git a/examples/injection_examples/plot_time_domain_data.py b/examples/injection_examples/plot_time_domain_data.py
index f3d056b7a9a594a328a200c84ac334ec50775093..eee3c8557bd447714bd0fe9025f838928f9cd542 100644
--- a/examples/injection_examples/plot_time_domain_data.py
+++ b/examples/injection_examples/plot_time_domain_data.py
@@ -15,7 +15,7 @@ label = 'example'
 
 injection_parameters = dict(
     mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=1000., iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=1000., theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
 
 waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
diff --git a/examples/injection_examples/roq_example.py b/examples/injection_examples/roq_example.py
index cd8475b6c81bc84b4623f37db70a6653026246da..2ad8c99d4d020289e1e148f63829c1dbcc0c6560 100644
--- a/examples/injection_examples/roq_example.py
+++ b/examples/injection_examples/roq_example.py
@@ -33,7 +33,7 @@ sampling_frequency = 2048
 
 injection_parameters = dict(
     chirp_mass=36., mass_ratio=0.9, a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=1000., iota=0.4, psi=0.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=1000., theta_jn=0.4, psi=0.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
 
 waveform_arguments = dict(waveform_approximant='IMRPhenomPv2',
@@ -63,7 +63,7 @@ search_waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
     parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters)
 
 priors = bilby.gw.prior.BBHPriorDict()
-for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'iota', 'phase', 'psi', 'ra',
+for key in ['a_1', 'a_2', 'tilt_1', 'tilt_2', 'theta_jn', 'phase', 'psi', 'ra',
             'dec', 'phi_12', 'phi_jl', 'luminosity_distance']:
     priors[key] = injection_parameters[key]
 priors.pop('mass_1')
diff --git a/examples/injection_examples/standard_15d_cbc_tutorial.py b/examples/injection_examples/standard_15d_cbc_tutorial.py
index aa17219acd2da42e20ec878c9a2b13ee03feb831..822c11b1002f1bdb1021e5d0d7daef57dc16a259 100644
--- a/examples/injection_examples/standard_15d_cbc_tutorial.py
+++ b/examples/injection_examples/standard_15d_cbc_tutorial.py
@@ -27,7 +27,7 @@ np.random.seed(88170235)
 # spins of both black holes (a, tilt, phi), etc.
 injection_parameters = dict(
     mass_1=36., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.5, tilt_2=1.0,
-    phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., iota=0.4, psi=2.659,
+    phi_12=1.7, phi_jl=0.3, luminosity_distance=2000., theta_jn=0.4, psi=2.659,
     phase=1.3, geocent_time=1126259642.413, ra=1.375, dec=-1.2108)
 
 # Fixed arguments passed into the source model
diff --git a/examples/injection_examples/using_gwin.py b/examples/injection_examples/using_gwin.py
index 1539195e60cf0e68f1b3d565f35161e91fbce84b..4d9c583ac6fa35f81b038232c566d2fbada5ad8a 100644
--- a/examples/injection_examples/using_gwin.py
+++ b/examples/injection_examples/using_gwin.py
@@ -29,7 +29,7 @@ label = 'using_gwin'
 # Search priors
 priors = dict()
 priors['distance'] = bilby.core.prior.Uniform(500, 2000, 'distance')
-priors['polarization'] = bilby.core.prior.Uniform(0, np.pi, 'iota')
+priors['polarization'] = bilby.core.prior.Uniform(0, np.pi, 'theta_jn')
 
 # Data variables
 seglen = 4
diff --git a/setup.cfg b/setup.cfg
index 2c867fe83d66a384df0dea28d12d225f8a60d90a..5a877a55c7907bbe1e011ebe111504508aa90b76 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -6,6 +6,8 @@ ignore = E129 W504 W605
 [tool:pytest]
 addopts =
     --ignore test/other_test.py
+    --ignore test/gw_example_test.py
+    --ignore test/example_test.py
 
 [metadata]
 license_file = LICENSE.md
diff --git a/test/detector_test.py b/test/detector_test.py
index edaa139fc74c12f694254986771decf16512ceee..07d107a074b9f534f1ab4bb45d747585519f4f7d 100644
--- a/test/detector_test.py
+++ b/test/detector_test.py
@@ -898,6 +898,13 @@ class TestInterferometerList(unittest.TestCase):
         with self.assertRaises(ValueError):
             self.ifo_list.inject_signal(injection_polarizations=None, waveform_generator=None)
 
+    def test_meta_data(self):
+        ifos_list = [self.ifo1, self.ifo2]
+        ifos = bilby.gw.detector.InterferometerList(ifos_list)
+        self.assertTrue(isinstance(ifos.meta_data, dict))
+        meta_data = {ifo.name: ifo.meta_data for ifo in ifos_list}
+        self.assertEqual(ifos.meta_data, meta_data)
+
     @patch.object(bilby.gw.waveform_generator.WaveformGenerator, 'frequency_domain_strain')
     def test_inject_signal_pol_none_calls_frequency_domain_strain(self, m):
         waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
diff --git a/test/gw_likelihood_test.py b/test/gw_likelihood_test.py
index b2d9a8128a0c2dcfa0089747dbd2d8a727a9422d..632aadac94536dcd425e7b2914ef1f63bcf17a94 100644
--- a/test/gw_likelihood_test.py
+++ b/test/gw_likelihood_test.py
@@ -10,7 +10,7 @@ class TestBasicGWTransient(unittest.TestCase):
         np.random.seed(500)
         self.parameters = dict(
             mass_1=31., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
-            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4,
+            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4,
             psi=2.659, phase=1.3, geocent_time=1126259642.413, ra=1.375,
             dec=-1.2108)
         self.interferometers = bilby.gw.detector.InterferometerList(['H1'])
@@ -73,7 +73,7 @@ class TestGWTransient(unittest.TestCase):
         self.sampling_frequency = 2048
         self.parameters = dict(
             mass_1=31., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
-            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4,
+            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4,
             psi=2.659, phase=1.3, geocent_time=1126259642.413, ra=1.375,
             dec=-1.2108)
         self.interferometers = bilby.gw.detector.InterferometerList(['H1'])
@@ -143,7 +143,7 @@ class TestTimeMarginalization(unittest.TestCase):
         self.sampling_frequency = 2048
         self.parameters = dict(
             mass_1=31., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
-            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4,
+            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4,
             psi=2.659, phase=1.3, geocent_time=1126259640, ra=1.375,
             dec=-1.2108)
 
@@ -242,7 +242,7 @@ class TestMarginalizedLikelihood(unittest.TestCase):
         self.sampling_frequency = 2048
         self.parameters = dict(
             mass_1=31., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
-            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4,
+            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4,
             psi=2.659, phase=1.3, geocent_time=1126259642.413, ra=1.375,
             dec=-1.2108)
 
@@ -306,7 +306,7 @@ class TestPhaseMarginalization(unittest.TestCase):
         self.sampling_frequency = 2048
         self.parameters = dict(
             mass_1=31., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
-            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4,
+            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4,
             psi=2.659, phase=1.3, geocent_time=1126259642.413, ra=1.375,
             dec=-1.2108)
 
@@ -369,7 +369,7 @@ class TestTimePhaseMarginalization(unittest.TestCase):
         self.sampling_frequency = 2048
         self.parameters = dict(
             mass_1=31., mass_2=29., a_1=0.4, a_2=0.3, tilt_1=0.0, tilt_2=0.0,
-            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., iota=0.4,
+            phi_12=1.7, phi_jl=0.3, luminosity_distance=4000., theta_jn=0.4,
             psi=2.659, phase=1.3, geocent_time=1126259642.413, ra=1.375,
             dec=-1.2108)
 
@@ -472,9 +472,9 @@ class TestROQLikelihood(unittest.TestCase):
         fnodes_quadratic = np.load(fnodes_quadratic_file).T
 
         self.test_parameters = dict(
-            mass_1=36.0, mass_2=36.0, a_1=0.0, a_2=0.0, tilt_1=0.0,
-            tilt_2=0.0, phi_12=1.7, phi_jl=0.3, luminosity_distance=5000.,
-            iota=0.4, psi=0.659, phase=1.3, geocent_time=1.2, ra=1.3, dec=-1.2)
+            mass_1=36.0, mass_2=36.0, a_1=0.0, a_2=0.0, tilt_1=0.0, tilt_2=0.0,
+            phi_12=1.7, phi_jl=0.3, luminosity_distance=5000., theta_jn=0.4,
+            psi=0.659, phase=1.3, geocent_time=1.2, ra=1.3, dec=-1.2)
 
         ifos = bilby.gw.detector.InterferometerList(['H1'])
         ifos.set_strain_data_from_power_spectral_densities(
@@ -509,7 +509,12 @@ class TestROQLikelihood(unittest.TestCase):
             interferometers=ifos, waveform_generator=roq_wfg,
             linear_matrix=linear_matrix_file,
             quadratic_matrix=quadratic_matrix_file, priors=self.priors)
-        pass
+
+        self.roq_phase_like = bilby.gw.likelihood.ROQGravitationalWaveTransient(
+            interferometers=ifos, waveform_generator=roq_wfg,
+            linear_matrix=linear_matrix_file,
+            quadratic_matrix=quadratic_matrix_file,
+            phase_marginalization=True, priors=self.priors.copy())
 
     def tearDown(self):
         pass
@@ -527,6 +532,20 @@ class TestROQLikelihood(unittest.TestCase):
         self.assertEqual(
             self.roq_likelihood.log_likelihood_ratio(), np.nan_to_num(-np.inf))
 
+    def test_phase_marginalisation_roq(self):
+        """Test phase marginalised likelihood matches brute force version"""
+        like = []
+        self.roq_likelihood.parameters.update(self.test_parameters)
+        phases = np.linspace(0, 2 * np.pi, 1000)
+        for phase in phases:
+            self.roq_likelihood.parameters['phase'] = phase
+            like.append(np.exp(self.roq_likelihood.log_likelihood_ratio()))
+
+        marg_like = np.log(np.trapz(like, phases) / (2 * np.pi))
+        self.roq_phase_like.parameters = self.test_parameters.copy()
+        self.assertAlmostEqual(
+            marg_like, self.roq_phase_like.log_likelihood_ratio(), delta=0.5)
+
 
 class TestBBHLikelihoodSetUp(unittest.TestCase):
 
diff --git a/test/gw_prior_test.py b/test/gw_prior_test.py
index f7c7454e9472f496b048045e1719b23f08321cfa..1eac2e76ca9311323edc6033668611cee5780a31 100644
--- a/test/gw_prior_test.py
+++ b/test/gw_prior_test.py
@@ -16,42 +16,174 @@ class TestBBHPriorDict(unittest.TestCase):
         self.base_directory =\
             '/'.join(os.path.dirname(
                 os.path.abspath(sys.argv[0])).split('/')[:-1])
-        self.filename = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'prior_files/binary_black_holes.prior')
-        self.default_prior = bilby.gw.prior.BBHPriorDict(
-            filename=self.filename)
+        self.filename = os.path.join(os.path.dirname(os.path.realpath(__file__)),
+                                     'prior_files/binary_black_holes.prior')
+        self.bbh_prior_dict = bilby.gw.prior.BBHPriorDict(filename=self.filename)
+        for key, value in self.bbh_prior_dict.items():
+            self.prior_dict[key] = value
 
     def tearDown(self):
         del self.prior_dict
         del self.filename
+        del self.bbh_prior_dict
+        del self.base_directory
 
     def test_create_default_prior(self):
         default = bilby.gw.prior.BBHPriorDict()
-        minima = all([self.default_prior[key].minimum == default[key].minimum
+        minima = all([self.bbh_prior_dict[key].minimum == default[key].minimum
                       for key in default.keys()])
-        maxima = all([self.default_prior[key].maximum == default[key].maximum
+        maxima = all([self.bbh_prior_dict[key].maximum == default[key].maximum
                       for key in default.keys()])
-        names = all([self.default_prior[key].name == default[key].name
+        names = all([self.bbh_prior_dict[key].name == default[key].name
                      for key in default.keys()])
 
         self.assertTrue(all([minima, maxima, names]))
 
     def test_create_from_dict(self):
-        bilby.gw.prior.BBHPriorDict(dictionary=self.prior_dict)
+        new_dict = bilby.gw.prior.BBHPriorDict(dictionary=self.prior_dict)
+        for key in self.bbh_prior_dict:
+            self.assertEqual(self.bbh_prior_dict[key], new_dict[key])
+
+    def test_redundant_priors_not_in_dict_before(self):
+        for prior in ['chirp_mass', 'total_mass', 'mass_ratio', 'symmetric_mass_ratio',
+                      'cos_tilt_1', 'cos_tilt_2', 'phi_1', 'phi_2', 'cos_theta_jn',
+                      'comoving_distance', 'redshift']:
+            self.assertTrue(self.bbh_prior_dict.test_redundancy(prior))
+
+    def test_redundant_priors_already_in_dict(self):
+        for prior in ['mass_1', 'mass_2', 'tilt_1', 'tilt_2',
+                      'phi_1', 'phi_2', 'theta_jn', 'luminosity_distance']:
+            self.assertTrue(self.bbh_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_masses(self):
+        del self.bbh_prior_dict['mass_2']
+        for prior in ['mass_2', 'chirp_mass', 'total_mass', 'mass_ratio',  'symmetric_mass_ratio']:
+            self.assertFalse(self.bbh_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_spin_magnitudes(self):
+        del self.bbh_prior_dict['a_2']
+        self.assertFalse(self.bbh_prior_dict.test_redundancy('a_2'))
+
+    def test_correct_not_redundant_priors_spin_tilt_1(self):
+        del self.bbh_prior_dict['tilt_1']
+        for prior in ['tilt_1', 'cos_tilt_1']:
+            self.assertFalse(self.bbh_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_spin_tilt_2(self):
+        del self.bbh_prior_dict['tilt_2']
+        for prior in ['tilt_2', 'cos_tilt_2']:
+            self.assertFalse(self.bbh_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_spin_azimuth(self):
+        del self.bbh_prior_dict['phi_12']
+        for prior in ['phi_1', 'phi_2', 'phi_12']:
+            self.assertFalse(self.bbh_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_inclination(self):
+        del self.bbh_prior_dict['theta_jn']
+        for prior in ['theta_jn', 'cos_theta_jn']:
+            self.assertFalse(self.bbh_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_distance(self):
+        del self.bbh_prior_dict['luminosity_distance']
+        for prior in ['luminosity_distance', 'comoving_distance',
+                      'redshift']:
+            self.assertFalse(self.bbh_prior_dict.test_redundancy(prior))
+
+    def test_add_unrelated_prior(self):
+        self.assertFalse(self.bbh_prior_dict.test_redundancy('abc'))
+
+    def test_test_has_redundant_priors(self):
+        self.assertFalse(self.bbh_prior_dict.test_has_redundant_keys())
+        for prior in ['chirp_mass', 'total_mass', 'mass_ratio', 'symmetric_mass_ratio',
+                      'cos_tilt_1', 'cos_tilt_2', 'phi_1', 'phi_2', 'cos_theta_jn',
+                      'comoving_distance', 'redshift']:
+            self.bbh_prior_dict[prior] = 0
+            self.assertTrue(self.bbh_prior_dict.test_has_redundant_keys())
+            del self.bbh_prior_dict[prior]
+
+
+class TestBNSPriorDict(unittest.TestCase):
 
-    def test_create_from_filename(self):
-        bilby.gw.prior.BBHPriorDict(filename=self.filename)
+    def setUp(self):
+        self.prior_dict = dict()
+        self.base_directory =\
+            '/'.join(os.path.dirname(
+                os.path.abspath(sys.argv[0])).split('/')[:-1])
+        self.filename = os.path.join(os.path.dirname(os.path.realpath(__file__)),
+                                     'prior_files/binary_neutron_stars.prior')
+        self.bns_prior_dict = bilby.gw.prior.BNSPriorDict(filename=self.filename)
+        for key, value in self.bns_prior_dict.items():
+            self.prior_dict[key] = value
 
-    def test_key_in_prior_not_redundant(self):
-        test = self.default_prior.test_redundancy('mass_1')
-        self.assertFalse(test)
+    def tearDown(self):
+        del self.prior_dict
+        del self.filename
+        del self.bns_prior_dict
+        del self.base_directory
 
-    def test_chirp_mass_redundant(self):
-        test = self.default_prior.test_redundancy('chirp_mass')
-        self.assertTrue(test)
+    def test_create_default_prior(self):
+        default = bilby.gw.prior.BNSPriorDict()
+        minima = all([self.bns_prior_dict[key].minimum == default[key].minimum
+                      for key in default.keys()])
+        maxima = all([self.bns_prior_dict[key].maximum == default[key].maximum
+                      for key in default.keys()])
+        names = all([self.bns_prior_dict[key].name == default[key].name
+                     for key in default.keys()])
 
-    def test_comoving_distance_redundant(self):
-        test = self.default_prior.test_redundancy('comoving_distance')
-        self.assertTrue(test)
+        self.assertTrue(all([minima, maxima, names]))
+
+    def test_create_from_dict(self):
+        new_dict = bilby.gw.prior.BNSPriorDict(dictionary=self.prior_dict)
+        self.assertDictEqual(self.bns_prior_dict, new_dict)
+
+    def test_redundant_priors_not_in_dict_before(self):
+        for prior in ['chirp_mass', 'total_mass', 'mass_ratio',
+                      'symmetric_mass_ratio', 'cos_theta_jn', 'comoving_distance',
+                      'redshift', 'lambda_tilde', 'delta_lambda']:
+            self.assertTrue(self.bns_prior_dict.test_redundancy(prior))
+
+    def test_redundant_priors_already_in_dict(self):
+        for prior in ['mass_1', 'mass_2', 'chi_1', 'chi_2',
+                      'theta_jn', 'luminosity_distance',
+                      'lambda_1', 'lambda_2']:
+            self.assertTrue(self.bns_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_masses(self):
+        del self.bns_prior_dict['mass_2']
+        for prior in ['mass_2', 'chirp_mass', 'total_mass', 'mass_ratio',  'symmetric_mass_ratio']:
+            self.assertFalse(self.bns_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_spin_magnitudes(self):
+        del self.bns_prior_dict['chi_2']
+        self.assertFalse(self.bns_prior_dict.test_redundancy('chi_2'))
+
+    def test_correct_not_redundant_priors_inclination(self):
+        del self.bns_prior_dict['theta_jn']
+        for prior in ['theta_jn', 'cos_theta_jn']:
+            self.assertFalse(self.bns_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_distance(self):
+        del self.bns_prior_dict['luminosity_distance']
+        for prior in ['luminosity_distance', 'comoving_distance',
+                      'redshift']:
+            self.assertFalse(self.bns_prior_dict.test_redundancy(prior))
+
+    def test_correct_not_redundant_priors_tidal(self):
+        del self.bns_prior_dict['lambda_1']
+        for prior in['lambda_1', 'lambda_tilde', 'delta_lambda']:
+            self.assertFalse(self.bns_prior_dict.test_redundancy(prior))
+
+    def test_add_unrelated_prior(self):
+        self.assertFalse(self.bns_prior_dict.test_redundancy('abc'))
+
+    def test_test_has_redundant_priors(self):
+        self.assertFalse(self.bns_prior_dict.test_has_redundant_keys())
+        for prior in ['chirp_mass', 'total_mass', 'mass_ratio', 'symmetric_mass_ratio',
+                      'cos_theta_jn', 'comoving_distance', 'redshift']:
+            self.bns_prior_dict[prior] = 0
+            self.assertTrue(self.bns_prior_dict.test_has_redundant_keys())
+            del self.bns_prior_dict[prior]
 
 
 class TestCalibrationPrior(unittest.TestCase):
diff --git a/test/gw_utils_test.py b/test/gw_utils_test.py
new file mode 100644
index 0000000000000000000000000000000000000000..bce6f97de0978fff982b1ee65d59e0559a8cfeab
--- /dev/null
+++ b/test/gw_utils_test.py
@@ -0,0 +1,188 @@
+from __future__ import absolute_import, division
+
+import unittest
+import os
+from shutil import rmtree
+
+import numpy as np
+import gwpy
+import lal
+import lalsimulation as lalsim
+
+import bilby
+from bilby.gw import utils as gwutils
+
+
+class TestGWUtils(unittest.TestCase):
+
+    def setUp(self):
+        self.outdir = 'outdir'
+        bilby.core.utils.check_directory_exists_and_if_not_mkdir(self.outdir)
+
+    def tearDown(self):
+        try:
+            rmtree(self.outdir)
+        except FileNotFoundError:
+            pass
+
+    def test_asd_from_freq_series(self):
+        freq_data = np.array([1, 2, 3])
+        df = 0.1
+        asd = gwutils.asd_from_freq_series(freq_data, df)
+        self.assertTrue(np.all(asd == freq_data * 2 * df ** 0.5))
+
+    def test_psd_from_freq_series(self):
+        freq_data = np.array([1, 2, 3])
+        df = 0.1
+        psd = gwutils.psd_from_freq_series(freq_data, df)
+        self.assertTrue(np.all(psd == (freq_data * 2 * df ** 0.5)**2))
+
+    def test_time_delay_from_geocenter(self):
+        det1 = np.array([0.1, 0.2, 0.3])
+        det2 = np.array([0.1, 0.2, 0.5])
+        ra = 0.5
+        dec = 0.2
+        time = 10
+        self.assertEqual(
+            gwutils.time_delay_geocentric(det1, det1, ra, dec, time), 0)
+        self.assertEqual(
+            gwutils.time_delay_geocentric(det1, det2, ra, dec, time),
+            1.3253791114055397e-10)
+
+    def test_get_polarization_tensor(self):
+        ra = 1
+        dec = 2.0
+        time = 10
+        psi = 0.1
+        for mode in ['plus', 'cross', 'breathing', 'longitudinal', 'x', 'y']:
+            p = gwutils.get_polarization_tensor(ra, dec, time, psi, mode)
+            self.assertEqual(p.shape, (3, 3))
+        with self.assertRaises(ValueError):
+            gwutils.get_polarization_tensor(ra, dec, time, psi, 'not-a-mode')
+
+    def test_inner_product(self):
+        aa = np.array([1, 2, 3])
+        bb = np.array([5, 6, 7])
+        frequency = np.array([0.2, 0.4, 0.6])
+        PSD = bilby.gw.detector.PowerSpectralDensity.from_aligo()
+        ip = gwutils.inner_product(aa, bb, frequency, PSD)
+        self.assertEqual(ip, 0)
+
+    def test_noise_weighted_inner_product(self):
+        aa = np.array([1e-23, 2e-23, 3e-23])
+        bb = np.array([5e-23, 6e-23, 7e-23])
+        frequency = np.array([100, 101, 102])
+        PSD = bilby.gw.detector.PowerSpectralDensity.from_aligo()
+        psd = PSD.power_spectral_density_interpolated(frequency)
+        duration = 4
+        nwip = gwutils.noise_weighted_inner_product(aa, bb, psd, duration)
+        self.assertEqual(nwip, 239.87768033598326)
+
+        self.assertEqual(
+            gwutils.optimal_snr_squared(aa, psd, duration),
+            gwutils.noise_weighted_inner_product(aa, aa, psd, duration))
+
+    def test_matched_filter_snr(self):
+        signal = np.array([1e-23, 2e-23, 3e-23])
+        frequency_domain_strain = np.array([5e-23, 6e-23, 7e-23])
+        frequency = np.array([100, 101, 102])
+        PSD = bilby.gw.detector.PowerSpectralDensity.from_aligo()
+        psd = PSD.power_spectral_density_interpolated(frequency)
+        duration = 4
+
+        mfsnr = gwutils.matched_filter_snr(
+            signal, frequency_domain_strain, psd, duration)
+        self.assertEqual(mfsnr, 25.510869054168282)
+
+    def test_get_event_time(self):
+        events = ['GW150914', 'GW151012', 'GW151226', 'GW170104', 'GW170608',
+                  'GW170729', 'GW170809', 'GW170814', 'GW170817', 'GW170818',
+                  'GW170823']
+        for event in events:
+            self.assertTrue(isinstance(gwutils.get_event_time(event), float))
+
+        self.assertTrue(gwutils.get_event_time('GW010290') is None)
+
+    def test_read_frame_file(self):
+        start_time = 0
+        end_time = 10
+        channel = 'H1:GDS-CALIB_STRAIN'
+        N = 100
+        times = np.linspace(start_time, end_time, N)
+        data = np.random.normal(0, 1, N)
+        ts = gwpy.timeseries.TimeSeries(data=data, times=times, t0=0)
+        ts.channel = gwpy.detector.Channel(channel)
+        ts.name = channel
+        filename = os.path.join(self.outdir, 'test.gwf')
+        ts.write(filename, format='gwf')
+
+        # Check reading without time limits
+        strain = gwutils.read_frame_file(
+            filename, start_time=None, end_time=None, channel=channel)
+        self.assertEqual(strain.channel.name, channel)
+        self.assertTrue(np.all(strain.value==data))
+
+        # Check reading with time limits
+        start_cut = 2
+        end_cut = 8
+        strain = gwutils.read_frame_file(
+            filename, start_time=start_cut, end_time=end_cut, channel=channel)
+        idxs = (times > start_cut) & (times < end_cut)
+        # Dropping the last element - for some reason gwpy drops the last element when reading in data
+        self.assertTrue(np.all(strain.value==data[idxs][:-1]))
+
+        # Check reading with unknown channels
+        strain = gwutils.read_frame_file(
+            filename, start_time=None, end_time=None)
+        self.assertTrue(np.all(strain.value==data))
+
+        # Check reading with incorrect channel
+        strain = gwutils.read_frame_file(
+            filename, start_time=None, end_time=None, channel='WRONG')
+        self.assertTrue(np.all(strain.value==data))
+
+        ts = gwpy.timeseries.TimeSeries(data=data, times=times, t0=0)
+        ts.name = 'NOT-A-KNOWN-CHANNEL'
+        ts.write(filename, format='gwf')
+        strain = gwutils.read_frame_file(
+            filename, start_time=None, end_time=None)
+        self.assertEqual(strain, None)
+
+    def test_convert_args_list_to_float(self):
+        self.assertEqual(
+            gwutils.convert_args_list_to_float(1, '2', 3.0), [1.0, 2.0, 3.0])
+        with self.assertRaises(ValueError):
+            gwutils.convert_args_list_to_float(1, '2', 'ten')
+
+    def test_lalsim_SimInspiralTransformPrecessingNewInitialConditions(self):
+        a = gwutils.lalsim_SimInspiralTransformPrecessingNewInitialConditions( 
+            0.1, 0, 0.6, 0.5, 0.6, 0.1, 0.8, 30.6, 23.2, 50, 0)
+        self.assertTrue(len(a) == 7)
+
+    def test_get_approximant(self):
+        with self.assertRaises(ValueError):
+            gwutils.lalsim_GetApproximantFromString(10)
+        self.assertEqual(gwutils.lalsim_GetApproximantFromString("IMRPhenomPV2"), 69)
+
+    def test_lalsim_SimInspiralChooseFDWaveform(self):
+        a = gwutils.lalsim_SimInspiralChooseFDWaveform(
+            35.2, 20.4, 0.1, 0.2, 0.2, 0.2, 0.2, 0.1, 1000, 2, 2.3,
+            45., 0.1, 10, 0.01, 10, 1000, 20, None, 69)
+        self.assertEqual(len(a), 2)
+        self.assertEqual(type(a[0]), lal.COMPLEX16FrequencySeries)
+        self.assertEqual(type(a[1]), lal.COMPLEX16FrequencySeries)
+
+        with self.assertRaises(ValueError):
+            a = gwutils.lalsim_SimInspiralChooseFDWaveform(
+                35.2, 20.4, 0.1, 0.2, 0.2, 0.2, 0.2, 0.1, 1000, 2, 2.3,
+                45., 0.1, 10, 0.01, 10, 1000, 20, None, 'IMRPhenomPV2')
+
+    def test_lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame(self):
+        version = lalsim.IMRPhenomPv2_V
+        a = gwutils.lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame(
+            25.6, 12.2, 50, 0.2, 0, 0, 0, 0, 0, 0, 0, version)
+        self.assertEqual(len(a), 7)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/test/make_standard_data.py b/test/make_standard_data.py
index fcab591cd11b14313fa87aa7d74483f58ac9823f..458ce69b4652ba7cd90fd90a3c746d055cc0b60f 100644
--- a/test/make_standard_data.py
+++ b/test/make_standard_data.py
@@ -22,7 +22,7 @@ simulation_parameters = dict(
     phi_12=0.,
     phi_jl=0.,
     luminosity_distance=100.,
-    iota=0.4,
+    theta_jn=0.4,
     phase=1.3,
     waveform_approximant='IMRPhenomPv2',
     reference_frequency=50.,
diff --git a/test/prior_files/binary_black_holes.prior b/test/prior_files/binary_black_holes.prior
index 5314635bd779f09098e3c8c2002f3367f046547d..e41cc73e2f765fbddd821693604ac15e302586f2 100644
--- a/test/prior_files/binary_black_holes.prior
+++ b/test/prior_files/binary_black_holes.prior
@@ -19,7 +19,7 @@ phi_jl =  Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi)
 luminosity_distance =  bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=5e3, unit='Mpc')
 dec =  Cosine(name='dec')
 ra =  Uniform(name='ra', minimum=0, maximum=2 * np.pi)
-iota =  Sine(name='iota')
-# cos_iota =  Uniform(name='cos_iota', minimum=-1, maximum=1)
+theta_jn =  Sine(name='theta_jn')
+# cos_theta_jn =  Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
 psi =  Uniform(name='psi', minimum=0, maximum=np.pi)
 phase =  Uniform(name='phase', minimum=0, maximum=2 * np.pi)
diff --git a/test/prior_files/binary_neutron_stars.prior b/test/prior_files/binary_neutron_stars.prior
new file mode 100644
index 0000000000000000000000000000000000000000..08ec88f2ffb8bd6536563cb20723aabb589c7c13
--- /dev/null
+++ b/test/prior_files/binary_neutron_stars.prior
@@ -0,0 +1,23 @@
+# These are the default priors we use for BNS systems.
+# Note that you may wish to use more specific mass and distance parameters.
+# These commands are all known to bilby.gw.prior.
+# Lines beginning "#" are ignored.
+mass_1 = Uniform(name='mass_1', minimum=1, maximum=2, unit='$M_{\\odot}$')
+mass_2 = Uniform(name='mass_2', minimum=1, maximum=2, unit='$M_{\\odot}$')
+# chirp_mass = Uniform(name='chirp_mass', minimum=0.87, maximum=1.74, unit='$M_{\\odot}$')
+# total_mass =  Uniform(name='total_mass', minimum=2, maximum=4, unit='$M_{\\odot}$')
+# mass_ratio =  Uniform(name='mass_ratio', minimum=0.5, maximum=1)
+# symmetric_mass_ratio =  Uniform(name='symmetric_mass_ratio', minimum=0.22, maximum=0.25)
+chi_1 =  bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1, 1), name='chi_1', latex_label='$\\chi_1$')
+chi_2 =  bilby.gw.prior.AlignedSpin(a_prior=Uniform(0, 0.05), z_prior=Uniform(-1, 1), name='chi_2', latex_label='$\\chi_2$')
+luminosity_distance =  bilby.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=10, maximum=500, unit='Mpc')
+dec =  Cosine(name='dec')
+ra =  Uniform(name='ra', minimum=0, maximum=2 * np.pi)
+theta_jn =  Sine(name='theta_jn')
+# cos_theta_jn =  Uniform(name='cos_theta_jn', minimum=-1, maximum=1)
+psi =  Uniform(name='psi', minimum=0, maximum=np.pi)
+phase =  Uniform(name='phase', minimum=0, maximum=2 * np.pi)
+lambda_1 = Uniform(name='lambda_1', minimum=0, maximum=3000 )
+lambda_2 = Uniform(name='lambda_2', minimum=0, maximum=3000 )
+# lambda_tilde = Uniform(name='lambda_tilde', minimum=0, maximum=5000)
+# delta_lambda = Uniform(name='delta_lambda', minimum=-5000, maximum=5000)
diff --git a/test/prior_test.py b/test/prior_test.py
index 4b85a66c625b43f1194fa79213278e807b65f399..d73d8d5b9001debc72c7cc7379c864a4619f90d1 100644
--- a/test/prior_test.py
+++ b/test/prior_test.py
@@ -380,7 +380,7 @@ class TestPriorDict(unittest.TestCase):
             dec=bilby.core.prior.Cosine(name='dec'),
             ra=bilby.core.prior.Uniform(
                 name='ra', minimum=0, maximum=2 * np.pi),
-            iota=bilby.core.prior.Sine(name='iota'),
+            theta_jn=bilby.core.prior.Sine(name='theta_jn'),
             psi=bilby.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi),
             phase=bilby.core.prior.Uniform(
                 name='phase', minimum=0, maximum=2 * np.pi)
@@ -436,7 +436,7 @@ class TestPriorDict(unittest.TestCase):
             dec=bilby.core.prior.Cosine(name='dec'),
             ra=bilby.core.prior.Uniform(
                 name='ra', minimum=0, maximum=2 * np.pi),
-            iota=bilby.core.prior.Sine(name='iota'),
+            theta_jn=bilby.core.prior.Sine(name='theta_jn'),
             psi=bilby.core.prior.Uniform(name='psi', minimum=0, maximum=np.pi),
             phase=bilby.core.prior.Uniform(
                 name='phase', minimum=0, maximum=2 * np.pi)
diff --git a/test/result_test.py b/test/result_test.py
index 6dbc5a2342391a57a9286b9d304894ec32b8746d..2aea035714b75ef31fa73ed65ce583365cdff2e5 100644
--- a/test/result_test.py
+++ b/test/result_test.py
@@ -45,10 +45,16 @@ class TestResult(unittest.TestCase):
         del self.result
         pass
 
-    def test_result_file_name(self):
+    def test_result_file_name_default(self):
         outdir = 'outdir'
         label = 'label'
         self.assertEqual(bilby.core.result.result_file_name(outdir, label),
+                         '{}/{}_result.json'.format(outdir, label))
+
+    def test_result_file_name_hdf5(self):
+        outdir = 'outdir'
+        label = 'label'
+        self.assertEqual(bilby.core.result.result_file_name(outdir, label, extension='hdf5'),
                          '{}/{}_result.h5'.format(outdir, label))
 
     def test_fail_save_and_load(self):
@@ -104,8 +110,8 @@ class TestResult(unittest.TestCase):
         with self.assertRaises(ValueError):
             _ = self.result.posterior
 
-    def test_save_and_load(self):
-        self.result.save_to_file()
+    def test_save_and_load_hdf5(self):
+        self.result.save_to_file(extension='hdf5')
         loaded_result = bilby.core.result.read_in_result(
             outdir=self.result.outdir, label=self.result.label)
         self.assertTrue(pd.DataFrame.equals
@@ -123,23 +129,61 @@ class TestResult(unittest.TestCase):
         self.assertEqual(self.result.priors['c'], loaded_result.priors['c'])
         self.assertEqual(self.result.priors['d'], loaded_result.priors['d'])
 
-    def test_save_and_dont_overwrite(self):
+    def test_save_and_load_default(self):
+        self.result.save_to_file()
+        loaded_result = bilby.core.result.read_in_result(
+            outdir=self.result.outdir, label=self.result.label)
+        self.assertTrue(np.array_equal
+                        (self.result.posterior.sort_values(by=['x']),
+                            loaded_result.posterior.sort_values(by=['x'])))
+        self.assertTrue(self.result.fixed_parameter_keys == loaded_result.fixed_parameter_keys)
+        self.assertTrue(self.result.search_parameter_keys == loaded_result.search_parameter_keys)
+        self.assertEqual(self.result.meta_data, loaded_result.meta_data)
+        self.assertEqual(self.result.injection_parameters, loaded_result.injection_parameters)
+        self.assertEqual(self.result.log_evidence, loaded_result.log_evidence)
+        self.assertEqual(self.result.log_noise_evidence, loaded_result.log_noise_evidence)
+        self.assertEqual(self.result.log_evidence_err, loaded_result.log_evidence_err)
+        self.assertEqual(self.result.log_bayes_factor, loaded_result.log_bayes_factor)
+        self.assertEqual(self.result.priors['x'], loaded_result.priors['x'])
+        self.assertEqual(self.result.priors['y'], loaded_result.priors['y'])
+        self.assertEqual(self.result.priors['c'], loaded_result.priors['c'])
+        self.assertEqual(self.result.priors['d'], loaded_result.priors['d'])
+
+    def test_save_and_dont_overwrite_default(self):
         shutil.rmtree(
-            '{}/{}_result.h5.old'.format(self.result.outdir, self.result.label),
+            '{}/{}_result.json.old'.format(self.result.outdir, self.result.label),
             ignore_errors=True)
         self.result.save_to_file(overwrite=False)
         self.result.save_to_file(overwrite=False)
+        self.assertTrue(os.path.isfile(
+            '{}/{}_result.json.old'.format(self.result.outdir, self.result.label)))
+
+    def test_save_and_dont_overwrite_hdf5(self):
+        shutil.rmtree(
+            '{}/{}_result.h5.old'.format(self.result.outdir, self.result.label),
+            ignore_errors=True)
+        self.result.save_to_file(overwrite=False, extension='hdf5')
+        self.result.save_to_file(overwrite=False, extension='hdf5')
         self.assertTrue(os.path.isfile(
             '{}/{}_result.h5.old'.format(self.result.outdir, self.result.label)))
 
-    def test_save_and_overwrite(self):
+    def test_save_and_overwrite_hdf5(self):
         shutil.rmtree(
             '{}/{}_result.h5.old'.format(self.result.outdir, self.result.label),
             ignore_errors=True)
+        self.result.save_to_file(overwrite=True, extension='hdf5')
+        self.result.save_to_file(overwrite=True, extension='hdf5')
+        self.assertFalse(os.path.isfile(
+            '{}/{}_result.h5.old'.format(self.result.outdir, self.result.label)))
+
+    def test_save_and_overwrite_default(self):
+        shutil.rmtree(
+            '{}/{}_result.json.old'.format(self.result.outdir, self.result.label),
+            ignore_errors=True)
         self.result.save_to_file(overwrite=True)
         self.result.save_to_file(overwrite=True)
         self.assertFalse(os.path.isfile(
-            '{}/{}_result.h5.old'.format(self.result.outdir, self.result.label)))
+            '{}/{}_result.json.old'.format(self.result.outdir, self.result.label)))
 
     def test_save_samples(self):
         self.result.save_posterior_samples()