diff --git a/examples/injection_examples/eccentric_GW150914_inspiral.py b/examples/injection_examples/eccentric_GW150914_inspiral.py
new file mode 100644
index 0000000000000000000000000000000000000000..b003c6b20999151030fc078ad57651e8e6bb0b29
--- /dev/null
+++ b/examples/injection_examples/eccentric_GW150914_inspiral.py
@@ -0,0 +1,113 @@
+#!/bin/python
+"""
+Tutorial to demonstrate running parameter estimation on a reduced parameter space for an injected eccentric binary
+black hole signal with masses & distnace similar to GW150914.
+
+This uses the same binary parameters that were used to make Figures 1, 2 & 5 in Lower et al. (2018) -> arXiv:1806.05350.
+
+For a more comprehensive look at what goes on in each step, refer to the "basic_tutorial.py" example.
+"""
+from __future__ import division, print_function
+
+import numpy as np
+
+import tupak
+
+import matplotlib.pyplot as plt
+
+duration = 64.
+sampling_frequency = 2048.
+
+outdir = 'outdir'
+label = 'eccentric_GW140914'
+tupak.core.utils.setup_logger(outdir=outdir, label=label)
+
+# Set up a random seed for result reproducibility.
+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)
+
+waveform_arguments = dict(waveform_approximant='EccentricFD', reference_frequency=10., minimum_frequency=10.)
+
+# Create the waveform_generator using the LAL eccentric black hole no spins source function
+waveform_generator = tupak.gw.WaveformGenerator(
+ duration=duration, sampling_frequency=sampling_frequency,
+ frequency_domain_source_model=tupak.gw.source.lal_eccentric_binary_black_hole_no_spins,
+ parameters=injection_parameters, waveform_arguments=waveform_arguments)
+
+hf_signal = waveform_generator.frequency_domain_strain()
+
+# Setting up three interferometers (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)) at their design sensitivities.
+# The maximum frequency is set just prior to the point at which the waveform model terminates. This is to avoid any biases
+# introduced from using a sharply terminating waveform model.
+minimum_frequency = 10.0
+maximum_frequency = 133.0
+
+def get_interferometer(name, injection_polarizations, injection_parameters, duration, sampling_frequency,
+ minimum_frequency, maximum_frequency, outdir):
+ """
+ Sets up the interferometers & injects the signal into them
+ """
+ start_time = injection_parameters['geocent_time'] + 2 - duration
+
+ ifo = tupak.gw.detector.get_empty_interferometer(name)
+ if name == 'V1':
+ ifo.power_spectral_density.set_from_power_spectral_density_file('AdV_psd.txt')
+ else:
+ ifo.power_spectral_density.set_from_power_spectral_density_file('aLIGO_ZERO_DET_high_P_psd.txt')
+
+ ifo.set_strain_data_from_power_spectral_density(sampling_frequency=sampling_frequency,
+ duration=duration, start_time=start_time)
+
+ injection_polarizations = ifo.inject_signal(parameters=injection_parameters,
+ injection_polarizations=injection_polarizations,
+ waveform_generator=waveform_generator)
+
+ signal = ifo.get_detector_response(injection_polarizations, injection_parameters)
+
+ ifo.minimum_frequency = minimum_frequency
+ ifo.maximum_frequency = maximum_frequency
+
+ ifo.plot_data(signal=signal, outdir=outdir, label=label)
+ ifo.save_data(outdir, label=label)
+
+ return ifo
+
+# IFOs = [tupak.gw.detector.get_interferometer_with_fake_noise_and_injection(name, injection_polarizations=hf_signal,
+# injection_parameters=injection_parameters, duration=duration,
+# sampling_frequency=sampling_frequency, outdir=outdir) for name in ['H1', 'L1', 'V1']]
+
+name = ['H1', 'L1', 'V1']
+IFOs = []
+
+for i in range(0,3):
+ IFOs.append(get_interferometer(name[i], injection_polarizations=hf_signal, injection_parameters=injection_parameters,
+ duration=duration, sampling_frequency=sampling_frequency,
+ minimum_frequency=minimum_frequency, maximum_frequency=maximum_frequency,
+ outdir=outdir))
+
+# Now we set up the priors on each of the binary parameters.
+priors = dict()
+priors["mass_1"] = tupak.core.prior.Uniform(name='mass_1', minimum=5, maximum=60)
+priors["mass_2"] = tupak.core.prior.Uniform(name='mass_2', minimum=5, maximum=60)
+priors["eccentricity"] = tupak.core.prior.PowerLaw(name='eccentricity', latex_label='$e$', alpha=-1, minimum=1e-4, maximum=0.4)
+priors["luminosity_distance"] = tupak.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=1e2, maximum=2e3)
+priors["dec"] = tupak.core.prior.Cosine(name='dec')
+priors["ra"] = tupak.core.prior.Uniform(name='ra', minimum=0, maximum=2 * np.pi)
+priors["iota"] = tupak.core.prior.Sine(name='iota')
+priors["psi"] = tupak.core.prior.Uniform(name='psi', minimum=0, maximum=2 * np.pi)
+priors["phase"] = tupak.core.prior.Uniform(name='phase', minimum=0, maximum=2 * np.pi)
+priors["geocent_time"] = tupak.core.prior.Uniform(1180002600.9, 1180002601.1, name='geocent_time')
+
+# Initialising the likelihood function.
+likelihood = tupak.gw.likelihood.GravitationalWaveTransient(interferometers=IFOs, waveform_generator=waveform_generator,
+ time_marginalization=False, phase_marginalization=False,
+ distance_marginalization=False, prior=priors)
+
+# Now we run sampler (PyMultiNest in our case).
+result = tupak.run_sampler(likelihood=likelihood, priors=priors, sampler='pymultinest', npoints=1000,
+ injection_parameters=injection_parameters, outdir=outdir, label=label)
+
+# And finally we make some plots of the output posteriors.
+result.plot_corner()
diff --git a/tupak/.version b/tupak/.version
new file mode 100644
index 0000000000000000000000000000000000000000..894d2803057a904be49f6af6079e0069bdd87b78
--- /dev/null
+++ b/tupak/.version
@@ -0,0 +1 @@
+0.2.1: (UNCLEAN) a85bcc8 2018-07-31 23:36:23 -0500
diff --git a/tupak/gw/source.py b/tupak/gw/source.py
index 86b26a157ef81deece15b313a99efd33be60bbb2..e7bdd2733f0ed3acce1003634f27c26aeafec141 100644
--- a/tupak/gw/source.py
+++ b/tupak/gw/source.py
@@ -108,6 +108,86 @@ def lal_binary_black_hole(
return {'plus': h_plus, 'cross': h_cross}
+def lal_eccentric_binary_black_hole_no_spins(
+ frequency_array, mass_1, mass_2, eccentricity, luminosity_distance, iota, phase, ra, dec,
+ geocent_time, psi, **kwargs):
+ """ Eccentric binary black hole waveform model using lalsimulation (EccentricFD)
+
+ Parameters
+ ----------
+ frequency_array: array_like
+ The frequencies at which we want to calculate the strain
+ mass_1: float
+ The mass of the heavier object in solar masses
+ mass_2: float
+ The mass of the lighter object in solar masses
+ eccentricity: float
+ The orbital eccentricity of the system
+ luminosity_distance: float
+ The luminosity distance in megaparsec
+ iota: float
+ Orbital inclination
+ phase: float
+ The phase at coalescence
+ ra: float
+ The right ascension of the binary
+ dec: float
+ The declination of the object
+ geocent_time: float
+ The time at coalescence
+ psi: float
+ Orbital polarisation
+ kwargs: dict
+ Optional keyword arguments
+
+ Returns
+ -------
+ dict: A dictionary with the plus and cross polarisation strain modes
+ """
+
+ waveform_kwargs = dict(waveform_approximant='EccentricFD', reference_frequency=10.0,
+ minimum_frequency=10.0)
+ waveform_kwargs.update(kwargs)
+ waveform_approximant = waveform_kwargs['waveform_approximant']
+ reference_frequency = waveform_kwargs['reference_frequency']
+ minimum_frequency = waveform_kwargs['minimum_frequency']
+
+ if mass_2 > mass_1:
+ return None
+
+ luminosity_distance = luminosity_distance * 1e6 * utils.parsec
+ mass_1 = mass_1 * utils.solar_mass
+ mass_2 = mass_2 * utils.solar_mass
+
+ spin_1x = 0.0
+ spin_1y = 0.0
+ spin_1z = 0.0
+ spin_2x = 0.0
+ spin_2y = 0.0
+ spin_2z = 0.0
+
+ longitude_ascending_nodes = 0.0
+ mean_per_ano = 0.0
+
+ waveform_dictionary = None
+
+ approximant = lalsim.GetApproximantFromString(waveform_approximant)
+
+ maximum_frequency = frequency_array[-1]
+ delta_frequency = frequency_array[1] - frequency_array[0]
+
+ hplus, hcross = lalsim.SimInspiralChooseFDWaveform(
+ mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y,
+ spin_2z, luminosity_distance, iota, phase,
+ longitude_ascending_nodes, eccentricity, mean_per_ano, delta_frequency,
+ minimum_frequency, maximum_frequency, reference_frequency,
+ waveform_dictionary, approximant)
+
+ h_plus = hplus.data.data
+ h_cross = hcross.data.data
+
+ return {'plus': h_plus, 'cross': h_cross}
+
def sinegaussian(frequency_array, hrss, Q, frequency, ra, dec, geocent_time, psi):
tau = Q / (np.sqrt(2.0)*np.pi*frequency)