allow precessing BNS

bilby/gw/conversion.py

@@ -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 =\
         convert_to_lal_binary_black_hole_parameters(converted_parameters)
 
-    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] *
-                    np.cos(converted_parameters[tilt]))
-            else:
-                converted_parameters['chi_{}'.format(idx)] = (
-                    converted_parameters[mag])
-
-    # 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

bilby/gw/prior.py

@@ -336,10 +336,12 @@ class BNSPriorDict(PriorDict):
             By default this generates many additional parameters, see
             BNSPriorDict.default_conversion_function
         """
-        if not aligned_spin:
-            logger.warning('Non-aligned spins not yet supported for BNS.')
+        if aligned_spin:
+            default_file = 'binary_neutron_stars.prior'
+        else:
+            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)
             logger.info('No prior given, using default BNS priors in {}.'.format(filename))
         elif filename is not None:
             if not os.path.isfile(filename):

bilby/gw/prior_files/precessing_binary_neutron_stars.prior

+# 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}$', 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 = bilby.gw.prior.UniformSourceFrame(name='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)

bilby/gw/source.py

@@ -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,
+        **kwargs):
     """ A Binary Neutron Star waveform model using lalsimulation
 
     Parameters
@@ -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,
+        pn_amplitude_order=0)
     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, 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(

examples/injection_examples/binary_neutron_star_example.py

@@ -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 = bilby.gw.prior.BNSPriorDict()
 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 = bilby.gw.GravitationalWaveTransient(
     interferometers=interferometers, waveform_generator=waveform_generator,
     time_marginalization=False, phase_marginalization=False,

test/conversion_test.py

@@ -250,13 +250,6 @@ class TestConvertToLALParams(unittest.TestCase):
     def test_lambda_1(self):
         self._conversion_to_component_tidal(['lambda_1'])
 
-    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.bns_convert()
-        self.assertEqual(self.parameters['chi_1'], chi_1)
-
 
 class TestDistanceTransformations(unittest.TestCase):
 

test/gw_source_test.py

@@ -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,
-            minimum_frequency=20.0)
+            waveform_approximant='IMRPhenomPv2_NRTidal',
+            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):
         self.parameters.update(self.waveform_kwargs)
         self.assertIsInstance(
             bilby.gw.source.lal_binary_neutron_star(
@@ -78,14 +79,6 @@ class TestLalBNS(unittest.TestCase):
             bilby.gw.source.lal_binary_neutron_star(
                 self.frequency_array, **self.parameters)
 
-    def test_fails_without_aligned_spins(self):
-        self.parameters.pop('chi_1')
-        self.parameters.pop('chi_2')
-        self.parameters.update(self.waveform_kwargs)
-        with self.assertRaises(TypeError):
-            bilby.gw.source.lal_binary_neutron_star(
-                self.frequency_array, **self.parameters)
-
 
 class TestEccentricLalBBH(unittest.TestCase):