Commit b2f00873 authored by Gregory Ashton's avatar Gregory Ashton

Merge branch 'precessing-bns' into 'master'

allow precessing BNS

See merge request !491
parents 76f4e5b0 0776452d
Pipeline #64717 passed with stages
in 16 minutes and 28 seconds
......@@ -254,7 +254,7 @@ def convert_to_lal_binary_neutron_star_parameters(parameters):
Mass: mass_1, mass_2
Spin: chi_1, chi_2, tilt_1, tilt_2, phi_12, phi_jl
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.
......@@ -277,19 +277,6 @@ def convert_to_lal_binary_neutron_star_parameters(parameters):
converted_parameters, added_keys =\
for idx in ['1', '2']:
mag = 'a_{}'.format(idx)
if mag in original_keys:
tilt = 'tilt_{}'.format(idx)
if tilt in original_keys:
converted_parameters['chi_{}'.format(idx)] = (
converted_parameters[mag] *
converted_parameters['chi_{}'.format(idx)] = (
# catch if tidal parameters aren't present
if not any([key in converted_parameters for key in
['lambda_1', 'lambda_2', 'lambda_tilde', 'delta_lambda']]):
converted_parameters['lambda_1'] = 0
......@@ -336,10 +336,12 @@ class BNSPriorDict(PriorDict):
By default this generates many additional parameters, see
if not aligned_spin:
logger.warning('Non-aligned spins not yet supported for BNS.')
if aligned_spin:
default_file = 'binary_neutron_stars.prior'
default_file = 'precessing_binary_neutron_stars.prior'
if dictionary is None and filename is None:
filename = os.path.join(os.path.dirname(__file__), 'prior_files', 'binary_neutron_stars.prior')
filename = os.path.join(os.path.dirname(__file__), 'prior_files', default_file)'No prior given, using default BNS priors in {}.'.format(filename))
elif filename is not None:
if not os.path.isfile(filename):
# 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
# Lines beginning "#" are ignored.
mass_1 = Uniform(name='mass_1', minimum=1, maximum=2, unit='$M_{\\odot}$', boundary=None)
mass_2 = Uniform(name='mass_2', minimum=1, maximum=2, unit='$M_{\\odot}$', boundary=None)
mass_ratio = Constraint(name='mass_ratio', minimum=0.125, maximum=1)
# chirp_mass = Uniform(name='chirp_mass', minimum=0.87, maximum=1.74, unit='$M_{\\odot}$', boundary=None)
# total_mass = Uniform(name='total_mass', minimum=2, maximum=4, unit='$M_{\\odot}$', boundary=None)
# mass_ratio = Uniform(name='mass_ratio', minimum=0.5, maximum=1, boundary=None)
# symmetric_mass_ratio = Uniform(name='symmetric_mass_ratio', minimum=0.22, maximum=0.25, boundary=None)
a_1 = Uniform(name='a_1', minimum=0, maximum=0.05, boundary='reflective')
a_2 = Uniform(name='a_2', minimum=0, maximum=0.05, boundary='reflective')
tilt_1 = Sine(name='tilt_1', boundary='reflective')
tilt_2 = Sine(name='tilt_2', boundary='reflective')
# cos_tilt_1 = Uniform(name='cos_tilt_1', minimum=-1, maximum=1, boundary=None)
# cos_tilt_2 = Uniform(name='cos_tilt_2', minimum=-1, maximum=1, boundary=None)
phi_12 = Uniform(name='phi_12', minimum=0, maximum=2 * np.pi, boundary='periodic')
phi_jl = Uniform(name='phi_jl', minimum=0, maximum=2 * np.pi, boundary='periodic')
luminosity_distance ='luminosity_distance', minimum=10, maximum=500, unit='Mpc', boundary=None)
dec = Cosine(name='dec', boundary='reflective')
ra = Uniform(name='ra', minimum=0, maximum=2 * np.pi, boundary='periodic')
theta_jn = Sine(name='theta_jn', boundary='reflective')
# cos_theta_jn = Uniform(name='cos_theta_jn', minimum=-1, maximum=1, boundary=None)
psi = Uniform(name='psi', minimum=0, maximum=np.pi, boundary='periodic')
phase = Uniform(name='phase', minimum=0, maximum=2 * np.pi, boundary='periodic')
lambda_1 = Uniform(name='lambda_1', minimum=0, maximum=3000, boundary=None)
lambda_2 = Uniform(name='lambda_2', minimum=0, maximum=3000, boundary=None)
# lambda_tilde = Uniform(name='lambda_tilde', minimum=0, maximum=5000, boundary=None)
# delta_lambda = Uniform(name='delta_lambda', minimum=-5000, maximum=5000, boundary=None)
......@@ -72,8 +72,9 @@ def lal_binary_black_hole(
def lal_binary_neutron_star(
frequency_array, mass_1, mass_2, luminosity_distance, chi_1, chi_2,
theta_jn, phase, lambda_1, lambda_2, **kwargs):
frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, lambda_1, lambda_2,
""" A Binary Neutron Star waveform model using lalsimulation
......@@ -86,27 +87,27 @@ def lal_binary_neutron_star(
The mass of the lighter object in solar masses
luminosity_distance: float
The luminosity distance in megaparsec
chi_1: float
Dimensionless aligned spin
chi_2: float
Dimensionless aligned spin
a_1: float
Dimensionless primary spin magnitude
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
Azimuthal angle between the total binary angular momentum and the
orbital angular momentum
theta_jn: 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
......@@ -115,20 +116,16 @@ def lal_binary_neutron_star(
dict: A dictionary with the plus and cross polarisation strain modes
waveform_kwargs = dict(
waveform_approximant='TaylorF2', reference_frequency=50.0,
waveform_approximant='IMRPhenomPv2_NRTidal', reference_frequency=50.0,
minimum_frequency=20.0, maximum_frequency=frequency_array[-1],
pn_spin_order=-1, pn_tidal_order=-1, pn_phase_order=-1, pn_amplitude_order=0)
a_1 = abs(chi_1)
a_2 = abs(chi_2)
tilt_1 = np.arccos(np.sign(chi_1))
tilt_2 = np.arccos(np.sign(chi_2))
pn_spin_order=-1, pn_tidal_order=-1, pn_phase_order=-1,
return _base_lal_cbc_fd_waveform(
frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
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)
a_1=a_1, a_2=a_2, tilt_1=tilt_1, tilt_2=tilt_2, phi_12=phi_12,
phi_jl=phi_jl, lambda_1=lambda_1, lambda_2=lambda_2, **waveform_kwargs)
def lal_eccentric_binary_black_hole_no_spins(
......@@ -39,7 +39,7 @@ 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',
waveform_arguments = dict(waveform_approximant='IMRPhenomPv2_NRTidal',
reference_frequency=50., minimum_frequency=40.0)
# Create the waveform_generator using a LAL Binary Neutron Star source function
......@@ -64,6 +64,8 @@ interferometers.inject_signal(parameters=injection_parameters,
# Load the default prior for binary neutron stars.
# We're going to sample in chirp_mass, symmetric_mass_ratio, lambda_tilde, and
# delta_lambda rather than mass_1, mass_2, lambda_1, and lambda_2.
# BNS have aligned spins by default, if you want to allow precessing spins
# pass aligned_spin=False to the BNSPriorDict
priors =
for key in ['psi', 'geocent_time', 'ra', 'dec', 'chi_1', 'chi_2',
'theta_jn', 'luminosity_distance', 'phase']:
......@@ -81,7 +83,7 @@ priors['delta_lambda'] = bilby.core.prior.Uniform(
-5000, 5000, name='delta_lambda')
# Initialise the likelihood by passing in the interferometer data (IFOs)
# and the waveoform generator
# and the waveform generator
likelihood =
interferometers=interferometers, waveform_generator=waveform_generator,
time_marginalization=False, phase_marginalization=False,
......@@ -250,13 +250,6 @@ class TestConvertToLALParams(unittest.TestCase):
def test_lambda_1(self):
def test_bns_spherical_spin_to_aligned(self):
self.parameters['a_1'] = -1
self.parameters['tilt_1'] = np.arccos(0.5)
chi_1 = self.parameters['a_1'] * np.cos(self.parameters['tilt_1'])
self.assertEqual(self.parameters['chi_1'], chi_1)
class TestDistanceTransformations(unittest.TestCase):
......@@ -47,11 +47,12 @@ class TestLalBNS(unittest.TestCase):
def setUp(self):
self.parameters = dict(
mass_1=1.4, mass_2=1.4, luminosity_distance=400.0, chi_1=0.4,
chi_2=0.3, theta_jn=1.7, phase=0.0, lambda_1=100.0, lambda_2=100.0)
mass_1=1.4, mass_2=1.4, luminosity_distance=400.0, a_1=0.4,
a_2=0.3, tilt_1=0.2, tilt_2=1.7, phi_jl=0.2, phi_12=0.9,
theta_jn=1.7, phase=0.0, lambda_1=100.0, lambda_2=100.0)
self.waveform_kwargs = dict(
waveform_approximant='TaylorF2', reference_frequency=50.0,
reference_frequency=50.0, minimum_frequency=20.0)
self.frequency_array = bilby.core.utils.create_frequency_series(2048, 4)
def tearDown(self):
......@@ -59,7 +60,7 @@ class TestLalBNS(unittest.TestCase):
del self.waveform_kwargs
del self.frequency_array
def test_lal_bns_works_runs_valid_parameters(self):
def test_lal_bns_runs_with_valid_parameters(self):
......@@ -78,14 +79,6 @@ class TestLalBNS(unittest.TestCase):
self.frequency_array, **self.parameters)
def test_fails_without_aligned_spins(self):
with self.assertRaises(TypeError):
self.frequency_array, **self.parameters)
class TestEccentricLalBBH(unittest.TestCase):
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment