diff --git a/examples/injection_examples/binary_neutron_star_example.py b/examples/injection_examples/binary_neutron_star_example.py
deleted file mode 100644
index 18263b1f7ceada3ac8700c5a14209196c2d38e9f..0000000000000000000000000000000000000000
--- a/examples/injection_examples/binary_neutron_star_example.py
+++ /dev/null
@@ -1,81 +0,0 @@
-#!/bin/python
-"""
-Tutorial to demonstrate running parameter estimation on a binary neutron star
-system taking into account tidal deformabilities.
-
-This example estimates the masses using a uniform prior in both component masses
-and also estimates the tidal deformabilities using a uniform prior in both
-tidal deformabilities
-"""
-
-from __future__ import division, print_function
-
-import numpy as np
-
-import tupak
-
-# Specify the output directory and the name of the simulation.
-outdir = 'outdir'
-label = 'bns_example'
-tupak.core.utils.setup_logger(outdir=outdir, label=label)
-
-# Set up a random seed for result reproducibility. This is optional!
-np.random.seed(88170235)
-
-# We are going to inject a binary neutron star waveform. We first establish a
-# dictionary of parameters that includes all of the different waveform
-# parameters, including masses of the two black holes (mass_1, mass_2),
-# spins of both black holes (a_1,a_2) , etc.
-injection_parameters = dict(
- mass_1=1.5, mass_2=1.3, a_1=0.0, a_2=0.0, luminosity_distance=50.,
- iota=0.4, psi=2.659, phase=1.3, geocent_time=1126259642.413,
- ra=1.375, dec=-1.2108, lambda1=400, lambda2=450)
-
-# Set the duration and sampling frequency of the data segment that we're going
-# to inject the signal into. For the
-# TaylorF2 waveform, we cut the signal close to the isco frequency
-duration = 8
-sampling_frequency = 2*1570.
-start_time = injection_parameters['geocent_time'] + 2 - duration
-
-# Fixed arguments passed into the source model. The analysis starts at 40 Hz.
-waveform_arguments = dict(waveform_approximant='TaylorF2',
- reference_frequency=50., minimum_frequency=40.0)
-
-# Create the waveform_generator using a LAL Binary Neutron Star source function
-waveform_generator = tupak.gw.WaveformGenerator(
- duration=duration, sampling_frequency=sampling_frequency,
- frequency_domain_source_model=tupak.gw.source.lal_binary_neutron_star,
- waveform_arguments=waveform_arguments)
-
-# Set up interferometers. In this case we'll use three interferometers
-# (LIGO-Hanford (H1), LIGO-Livingston (L1), and Virgo (V1)).
-# These default to their design sensitivity and start at 40 Hz.
-interferometers = tupak.gw.detector.InterferometerList(['H1', 'L1', 'V1'])
-for interferometer in interferometers:
- interferometer.minimum_frequency = 40
-interferometers.set_strain_data_from_power_spectral_densities(
- sampling_frequency=sampling_frequency, duration=duration,
- start_time=start_time)
-interferometers.inject_signal(parameters=injection_parameters,
- waveform_generator=waveform_generator)
-
-priors = tupak.gw.prior.BNSPriorSet()
-for key in ['a_1', 'a_2', 'psi', 'geocent_time', 'ra', 'dec',
- 'iota', 'luminosity_distance', 'phase']:
- priors[key] = injection_parameters[key]
-
-# Initialise the likelihood by passing in the interferometer data (IFOs)
-# and the waveoform generator
-likelihood = tupak.gw.GravitationalWaveTransient(
- interferometers=interferometers, waveform_generator=waveform_generator,
- time_marginalization=False, phase_marginalization=False,
- distance_marginalization=False, prior=priors)
-
-# Run sampler. In this case we're going to use the `nestle` sampler
-result = tupak.run_sampler(
- likelihood=likelihood, priors=priors, sampler='nestle', npoints=1000,
- injection_parameters=injection_parameters, outdir=outdir, label=label)
-
-result.plot_corner()
-
diff --git a/test/conversion_tests.py b/test/conversion_tests.py
index f0e4409f2d360c93ba93a96f6ac33c7b4089dc88..3130af9534630b925550bc4ca3082ed9f63cfd11 100644
--- a/test/conversion_tests.py
+++ b/test/conversion_tests.py
@@ -8,35 +8,14 @@ import numpy as np
class TestBasicConversions(unittest.TestCase):
def setUp(self):
- self.mass_1 = 1.4
- self.mass_2 = 1.3
- self.mass_ratio = 13/14
- self.total_mass = 2.7
- self.chirp_mass = (1.4 * 1.3)**0.6 / 2.7**0.2
- self.symmetric_mass_ratio = (1.4 * 1.3) / 2.7**2
+ self.mass_1 = 20
+ self.mass_2 = 10
+ self.mass_ratio = 0.5
+ self.total_mass = 30
+ self.chirp_mass = 200**0.6 / 30**0.2
+ self.symmetric_mass_ratio = 2/9
self.cos_angle = -1
self.angle = np.pi
- self.lambda_1 = 300
- self.lambda_2 = 300 * (14 / 13)**5
- self.lambda_tilde = 8 / 13 * (
- (1 + 7 * self.symmetric_mass_ratio
- - 31 * self.symmetric_mass_ratio**2)
- * (self.lambda_1 + self.lambda_2)
- + (1 - 4 * self.symmetric_mass_ratio)**0.5
- * (1 + 9 * self.symmetric_mass_ratio
- - 11 * self.symmetric_mass_ratio**2)
- * (self.lambda_1 - self.lambda_2)
- )
- self.delta_lambda = 1 / 2 * (
- (1 - 4 * self.symmetric_mass_ratio)**0.5
- * (1 - 13272 / 1319 * self.symmetric_mass_ratio
- + 8944 / 1319 * self.symmetric_mass_ratio**2)
- * (self.lambda_1 + self.lambda_2)
- + (1 - 15910 / 1319 * self.symmetric_mass_ratio
- + 32850 / 1319 * self.symmetric_mass_ratio**2
- + 3380 / 1319 * self.symmetric_mass_ratio**3)
- * (self.lambda_1 - self.lambda_2)
- )
def tearDown(self):
del self.mass_1
@@ -48,8 +27,7 @@ class TestBasicConversions(unittest.TestCase):
def test_total_mass_and_mass_ratio_to_component_masses(self):
mass_1, mass_2 = tupak.gw.conversion.total_mass_and_mass_ratio_to_component_masses(self.mass_ratio, self.total_mass)
- self.assertTrue(all([abs(mass_1 - self.mass_1) < 1e-5,
- abs(mass_2 - self.mass_2) < 1e-5]))
+ self.assertTupleEqual((mass_1, mass_2), (self.mass_1, self.mass_2))
def test_symmetric_mass_ratio_to_mass_ratio(self):
mass_ratio = tupak.gw.conversion.symmetric_mass_ratio_to_mass_ratio(self.symmetric_mass_ratio)
@@ -83,20 +61,6 @@ class TestBasicConversions(unittest.TestCase):
mass_ratio = tupak.gw.conversion.mass_1_and_chirp_mass_to_mass_ratio(self.mass_1, self.chirp_mass)
self.assertAlmostEqual(self.mass_ratio, mass_ratio)
- def test_lambda_tilde_to_lambda_1_lambda_2(self):
- lambda_1, lambda_2 =\
- tupak.gw.conversion.lambda_tilde_to_lambda_1_lambda_2(
- self.lambda_tilde, self.mass_1, self.mass_2)
- self.assertTrue(all([abs(self.lambda_1 - lambda_1) < 1e-5,
- abs(self.lambda_2 - lambda_2) < 1e-5]))
-
- def test_lambda_tilde_delta_lambda_to_lambda_1_lambda_2(self):
- lambda_1, lambda_2 =\
- tupak.gw.conversion.lambda_tilde_delta_lambda_to_lambda_1_lambda_2(
- self.lambda_tilde, self.delta_lambda, self.mass_1, self.mass_2)
- self.assertTrue(all([abs(self.lambda_1 - lambda_1) < 1e-5,
- abs(self.lambda_2 - lambda_2) < 1e-5]))
-
class TestConvertToLALBBHParams(unittest.TestCase):
diff --git a/tupak/gw/conversion.py b/tupak/gw/conversion.py
index 561d35c241794b0fd346155db616da0e269c6eea..feb7a0b6a483c7fd7193420ecd5f65cd2ec18219 100644
--- a/tupak/gw/conversion.py
+++ b/tupak/gw/conversion.py
@@ -171,71 +171,6 @@ def convert_to_lal_binary_black_hole_parameters(parameters, search_keys, remove=
return converted_parameters, added_keys
-def convert_to_lal_binary_neutron_star_parameters(parameters, search_keys, remove=True):
- """
- Convert parameters we have into parameters we need.
-
- This is defined by the parameters of tupak.source.lal_binary_black_hole()
-
-
- Mass: mass_1, mass_2
- Spin: a_1, a_2, tilt_1, tilt_2, phi_12, phi_jl
- Extrinsic: luminosity_distance, theta_jn, phase, ra, dec, geocent_time, psi
-
- This involves popping a lot of things from parameters.
- The keys in added_keys should be popped after evaluating the waveform.
-
- Parameters
- ----------
- parameters: dict
- dictionary of parameter values to convert into the required parameters
- search_keys: list
- parameters which are needed for the waveform generation
- remove: bool, optional
- Whether or not to remove the extra key, necessary for sampling, default=True.
-
- Return
- ------
- converted_parameters: dict
- dict of the required parameters
- added_keys: list
- keys which are added to parameters during function call
- """
- converted_parameters = parameters.copy()
- converted_parameters, added_keys = convert_to_lal_binary_black_hole_parameters(
- converted_parameters, search_keys, remove=remove)
-
- if 'lambda_1' not in search_keys and 'lambda_2' not in search_keys:
- if 'delta_lambda' in converted_parameters.keys():
- converted_parameters['lambda_1'], converted_parameters['lambda_2'] =\
- lambda_tilde_delta_lambda_to_lambda_1_lambda_2(
- converted_parameters['lambda_tilde'], parameters['delta_lambda'],
- converted_parameters['mass_1'], converted_parameters['mass_2'])
- added_keys.append('lambda_1')
- added_keys.append('lambda_2')
- elif 'lambda_tilde' in converted_parameters.keys():
- converted_parameters['lambda_1'], converted_parameters['lambda_2'] =\
- lambda_tilde_to_lambda_1_lambda_2(
- converted_parameters['lambda_tilde'],
- converted_parameters['mass_1'], converted_parameters['mass_2'])
- added_keys.append('lambda_1')
- added_keys.append('lambda_2')
- if 'lambda_2' not in converted_parameters.keys():
- converted_parameters['lambda_2'] =\
- converted_parameters['lambda_1']\
- * converted_parameters['mass_1']**5\
- / converted_parameters['mass_2']**5
- added_keys.append('lambda_2')
- elif converted_parameters['lambda_2'] is None:
- converted_parameters['lambda_2'] =\
- converted_parameters['lambda_1']\
- * converted_parameters['mass_1']**5\
- / converted_parameters['mass_2']**5
- added_keys.append('lambda_2')
-
- return converted_parameters, added_keys
-
-
def total_mass_and_mass_ratio_to_component_masses(mass_ratio, total_mass):
"""
Convert total mass and mass ratio of a binary to its component masses.
@@ -425,85 +360,6 @@ def mass_1_and_chirp_mass_to_mass_ratio(mass_1, chirp_mass):
return mass_ratio
-def lambda_tilde_delta_lambda_to_lambda_1_lambda_2(
- lambda_tilde, delta_lambda, mass_1, mass_2):
- """
- Convert from dominant tidal terms to individual tidal parameters.
-
- See, e.g., Wade et al., https://arxiv.org/pdf/1402.5156.pdf.
-
- Parameters
- ----------
- lambda_tilde: float
- Dominant tidal term.
- delta_lambda: float
- Secondary tidal term.
- mass_1: float
- Mass of more massive neutron star.
- mass_2: float
- Mass of less massive neutron star.
-
- Return
- ------
- lambda_1: float
- Tidal parameter of more massive neutron star.
- lambda_2: float
- Tidal parameter of less massive neutron star.
- """
- eta = component_masses_to_symmetric_mass_ratio(mass_1, mass_2)
- coefficient_1 = (1 + 7 * eta - 31 * eta**2)
- coefficient_2 = (1 - 4 * eta)**0.5 * (1 + 9 * eta - 11 * eta**2)
- coefficient_3 = (1 - 4 * eta)**0.5\
- * (1 - 13272 / 1319 * eta + 8944 / 1319 * eta**2)
- coefficient_4 = (1 - 15910 / 1319 * eta + 32850 / 1319 * eta**2
- + 3380 / 1319 * eta**3)
- lambda_1 =\
- (13 * lambda_tilde / 8 * (coefficient_3 - coefficient_4)
- - 2 * delta_lambda * (coefficient_1 - coefficient_2))\
- / ((coefficient_1 + coefficient_2) * (coefficient_3 - coefficient_4)
- - (coefficient_1 - coefficient_2) * (coefficient_3 + coefficient_4))
- lambda_2 =\
- (13 * lambda_tilde / 8 * (coefficient_3 + coefficient_4)
- - 2 * delta_lambda * (coefficient_1 + coefficient_2)) \
- / ((coefficient_1 - coefficient_2) * (coefficient_3 + coefficient_4)
- - (coefficient_1 + coefficient_2) * (coefficient_3 - coefficient_4))
- return lambda_1, lambda_2
-
-
-def lambda_tilde_to_lambda_1_lambda_2(
- lambda_tilde, mass_1, mass_2):
- """
- Convert from dominant tidal term to individual tidal parameters
- assuming lambda_1 * mass_1**5 = lambda_2 * mass_2**5.
-
- See, e.g., Wade et al., https://arxiv.org/pdf/1402.5156.pdf.
-
- Parameters
- ----------
- lambda_tilde: float
- Dominant tidal term.
- mass_1: float
- Mass of more massive neutron star.
- mass_2: float
- Mass of less massive neutron star.
-
- Return
- ------
- lambda_1: float
- Tidal parameter of more massive neutron star.
- lambda_2: float
- Tidal parameter of less massive neutron star.
- """
- eta = component_masses_to_symmetric_mass_ratio(mass_1, mass_2)
- q = mass_2 / mass_1
- lambda_1 = 13 / 8 * lambda_tilde / (
- (1 + 7 * eta - 31 * eta**2) * (1 + q**-5)
- + (1 - 4 * eta)**0.5 * (1 + 9 * eta - 11 * eta**2) * (1 - q**-5)
- )
- lambda_2 = lambda_1 / q**5
- return lambda_1, lambda_2
-
-
def generate_all_bbh_parameters(sample, likelihood=None, priors=None):
"""
From either a single sample or a set of samples fill in all missing BBH parameters, in place.
diff --git a/tupak/gw/prior.py b/tupak/gw/prior.py
index 9d5acf7354896873517dd9ff180df3706cd4ee25..a9bbc9f8c94118c9a5e41e7458768c98e5544f65 100644
--- a/tupak/gw/prior.py
+++ b/tupak/gw/prior.py
@@ -96,45 +96,6 @@ class BBHPriorSet(PriorSet):
return redundant
-class BNSPriorSet(PriorSet):
-
- def __init__(self, dictionary=None, filename=None):
- """ Initialises a Prior set for Binary Neutron Stars
-
- Parameters
- ----------
- dictionary: dict, optional
- See superclass
- filename: str, optional
- See superclass
- """
- if dictionary is None and filename is None:
- filename = os.path.join(os.path.dirname(__file__), 'prior_files', 'binary_neutron_stars.prior')
- logger.info('No prior given, using default BNS priors in {}.'.format(filename))
- elif filename is not None:
- if not os.path.isfile(filename):
- filename = os.path.join(os.path.dirname(__file__), 'prior_files', filename)
- PriorSet.__init__(self, dictionary=dictionary, filename=filename)
-
- def test_redundancy(self, key):
- bbh_redundancy = BBHPriorSet().test_redundancy(key)
- if bbh_redundancy:
- return True
- redundant = False
-
- 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)))
- elif len(tidal_parameters.intersection(self)) == 2:
- redundant = True
- return redundant
-
-
Prior._default_latex_labels = {
'mass_1': '$m_1$',
'mass_2': '$m_2$',
@@ -157,11 +118,7 @@ Prior._default_latex_labels = {
'cos_iota': '$\cos\iota$',
'psi': '$\psi$',
'phase': '$\phi$',
- 'geocent_time': '$t_c$',
- 'lambda_1': '$\\Lambda_1$',
- 'lambda_2': '$\\Lambda_2$',
- 'lambda_tilde': '$\\tilde{\\Lambda}$',
- 'delta_lambda': '$\\delta\\Lambda$'}
+ 'geocent_time': '$t_c$'}
class CalibrationPriorSet(PriorSet):
diff --git a/tupak/gw/prior_files/binary_neutron_stars.prior b/tupak/gw/prior_files/binary_neutron_stars.prior
deleted file mode 100644
index f7b427c3c324c00ea8246d5150e6df9362f709c9..0000000000000000000000000000000000000000
--- a/tupak/gw/prior_files/binary_neutron_stars.prior
+++ /dev/null
@@ -1,23 +0,0 @@
-# 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 tupak.gw.prior.
-# Lines beginning "#" are ignored.
-mass_1 = Uniform(name='mass_1', minimum=1, maximum=2)
-mass_2 = Uniform(name='mass_2', minimum=1, maximum=2)
-# chirp_mass = Uniform(name='chirp_mass', minimum=0.87, maximum=1.74)
-# total_mass = Uniform(name='total_mass', minimum=2, maximum=4)
-# 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)
-a_1 = Uniform(name='a_1', minimum= -0.05, maximum= 0.05)
-a_2 = Uniform(name='a_2', minimum= -0.05, maximum= 0.05)
-luminosity_distance = tupak.gw.prior.UniformComovingVolume(name='luminosity_distance', minimum=10, maximum=500)
-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)
-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/tupak/gw/source.py b/tupak/gw/source.py
index 770e53b885f1698a900c30dc78e496096810ea4d..760ee473cc6da49c05e3ebcba37af1bf1e6db5db 100644
--- a/tupak/gw/source.py
+++ b/tupak/gw/source.py
@@ -7,7 +7,6 @@ from tupak.core import utils
try:
import lalsimulation as lalsim
- import lal
except ImportError:
logger.warning("You do not have lalsuite installed currently. You will "
" not be able to use some of the prebuilt functions.")
@@ -252,96 +251,3 @@ def supernova_pca_model(
pc_coeff4 * pc4 + pc_coeff5 * pc5)
return {'plus': h_plus, 'cross': h_cross}
-
-def lal_binary_neutron_star(
- frequency_array, mass_1, mass_2, luminosity_distance, a_1, a_2,
- iota, phase, lambda_1, lambda_2, ra, dec, geocent_time, psi, **kwargs):
- """ A Binary Neutron Star waveform model using lalsimulation
-
- 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
- luminosity_distance: float
- The luminosity distance in megaparsec
- a_1: float
- Dimensionless spin magnitude
- a_2: float
- Dimensionless secondary spin magnitude.
- 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
- lambda_1: float
- Dimensionless tidal deformability of mass_1
- lambda_2: float
- Dimensionless tidal deformability of mass_2
-
- kwargs: dict
- Optional keyword arguments
-
- Returns
- -------
- dict: A dictionary with the plus and cross polarisation strain modes
- """
-
- waveform_kwargs = dict(waveform_approximant='TaylorF2', reference_frequency=50.0,
- minimum_frequency=20.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
- spin_1y = 0
- spin_1z = a_1
- spin_2x = 0
- spin_2y = 0
- spin_2z = a_2
-
- longitude_ascending_nodes = 0.0
- eccentricity = 0.0
- mean_per_ano = 0.0
-
- waveform_dictionary = lal.CreateDict()
- lalsim.SimInspiralWaveformParamsInsertTidalLambda1(waveform_dictionary, lambda_1)
- lalsim.SimInspiralWaveformParamsInsertTidalLambda2(waveform_dictionary, lambda_2)
-
- 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
-
- h_plus = h_plus[:len(frequency_array)]
- h_cross = h_cross[:len(frequency_array)]
-
- return {'plus': h_plus, 'cross': h_cross}