Merge branch 'higher-order-modes' into 'master'

Higher order modes

See merge request !418

bilby/gw/source.py

@@ -62,7 +62,8 @@ def lal_binary_black_hole(
     """
     waveform_kwargs = dict(
         waveform_approximant='IMRPhenomPv2', reference_frequency=50.0,
-        minimum_frequency=20.0, maximum_frequency=frequency_array[-1])
+        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)
     waveform_kwargs.update(kwargs)
     return _base_lal_cbc_fd_waveform(
         frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
@@ -116,7 +117,9 @@ def lal_binary_neutron_star(
     """
     waveform_kwargs = dict(
         waveform_approximant='TaylorF2', reference_frequency=50.0,
-        minimum_frequency=20.0, maximum_frequency=frequency_array[-1])
+        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))
@@ -159,7 +162,8 @@ def lal_eccentric_binary_black_hole_no_spins(
     """
     waveform_kwargs = dict(
         waveform_approximant='EccentricFD', reference_frequency=10.0,
-        minimum_frequency=10.0, maximum_frequency=frequency_array[-1])
+        minimum_frequency=10.0, maximum_frequency=frequency_array[-1],
+        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,
@@ -216,6 +220,20 @@ def _base_lal_cbc_fd_waveform(
     reference_frequency = waveform_kwargs['reference_frequency']
     minimum_frequency = waveform_kwargs['minimum_frequency']
     maximum_frequency = waveform_kwargs['maximum_frequency']
+    pn_spin_order = waveform_kwargs['pn_spin_order']
+    pn_tidal_order = waveform_kwargs['pn_tidal_order']
+    pn_phase_order = waveform_kwargs['pn_phase_order']
+    pn_amplitude_order = waveform_kwargs['pn_amplitude_order']
+
+    approximant = lalsim_GetApproximantFromString(waveform_approximant)
+
+    if (pn_amplitude_order != 0):
+        start_frequency = lalsim.SimInspiralfLow2fStart(minimum_frequency,
+                                                        int(pn_amplitude_order),
+                                                        approximant)
+    else:
+        start_frequency = minimum_frequency
+
     delta_frequency = frequency_array[1] - frequency_array[0]
 
     frequency_bounds = ((frequency_array >= minimum_frequency) *
@@ -243,11 +261,13 @@ def _base_lal_cbc_fd_waveform(
     mean_per_ano = 0.0
 
     waveform_dictionary = lal.CreateDict()
+    lalsim.SimInspiralWaveformParamsInsertPNSpinOrder(waveform_dictionary, int(pn_spin_order))
+    lalsim.SimInspiralWaveformParamsInsertPNTidalOrder(waveform_dictionary, int(pn_tidal_order))
+    lalsim.SimInspiralWaveformParamsInsertPNPhaseOrder(waveform_dictionary, int(pn_phase_order))
+    lalsim.SimInspiralWaveformParamsInsertPNAmplitudeOrder(waveform_dictionary, int(pn_amplitude_order))
     lalsim_SimInspiralWaveformParamsInsertTidalLambda1(waveform_dictionary, lambda_1)
     lalsim_SimInspiralWaveformParamsInsertTidalLambda2(waveform_dictionary, lambda_2)
 
-    approximant = lalsim_GetApproximantFromString(waveform_approximant)
-
     if lalsim.SimInspiralImplementedFDApproximants(approximant):
         wf_func = lalsim_SimInspiralChooseFDWaveform
     else:
@@ -256,7 +276,7 @@ def _base_lal_cbc_fd_waveform(
         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,
+        start_frequency, maximum_frequency, reference_frequency,
         waveform_dictionary, approximant)
 
     h_plus = np.zeros_like(frequency_array, dtype=np.complex)