diff --git a/test/conversion_tests.py b/test/conversion_tests.py
index 3130af9534630b925550bc4ca3082ed9f63cfd11..f0e4409f2d360c93ba93a96f6ac33c7b4089dc88 100644
--- a/test/conversion_tests.py
+++ b/test/conversion_tests.py
@@ -8,14 +8,35 @@ import numpy as np
class TestBasicConversions(unittest.TestCase):
def setUp(self):
- 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.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.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
@@ -27,7 +48,8 @@ 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.assertTupleEqual((mass_1, mass_2), (self.mass_1, self.mass_2))
+ self.assertTrue(all([abs(mass_1 - self.mass_1) < 1e-5,
+ abs(mass_2 - self.mass_2) < 1e-5]))
def test_symmetric_mass_ratio_to_mass_ratio(self):
mass_ratio = tupak.gw.conversion.symmetric_mass_ratio_to_mass_ratio(self.symmetric_mass_ratio)
@@ -61,6 +83,20 @@ 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 feb7a0b6a483c7fd7193420ecd5f65cd2ec18219..561d35c241794b0fd346155db616da0e269c6eea 100644
--- a/tupak/gw/conversion.py
+++ b/tupak/gw/conversion.py
@@ -171,6 +171,71 @@ 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.
@@ -360,6 +425,85 @@ 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 a9bbc9f8c94118c9a5e41e7458768c98e5544f65..9d5acf7354896873517dd9ff180df3706cd4ee25 100644
--- a/tupak/gw/prior.py
+++ b/tupak/gw/prior.py
@@ -96,6 +96,45 @@ 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$',
@@ -118,7 +157,11 @@ Prior._default_latex_labels = {
'cos_iota': '$\cos\iota$',
'psi': '$\psi$',
'phase': '$\phi$',
- 'geocent_time': '$t_c$'}
+ 'geocent_time': '$t_c$',
+ 'lambda_1': '$\\Lambda_1$',
+ 'lambda_2': '$\\Lambda_2$',
+ 'lambda_tilde': '$\\tilde{\\Lambda}$',
+ 'delta_lambda': '$\\delta\\Lambda$'}
class CalibrationPriorSet(PriorSet):
diff --git a/tupak/gw/prior_files/binary_neutron_stars.prior b/tupak/gw/prior_files/binary_neutron_stars.prior
new file mode 100644
index 0000000000000000000000000000000000000000..f7b427c3c324c00ea8246d5150e6df9362f709c9
--- /dev/null
+++ b/tupak/gw/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 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 760ee473cc6da49c05e3ebcba37af1bf1e6db5db..770e53b885f1698a900c30dc78e496096810ea4d 100644
--- a/tupak/gw/source.py
+++ b/tupak/gw/source.py
@@ -7,6 +7,7 @@ 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.")
@@ -251,3 +252,96 @@ 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}