From c7c16d90269aa388352b8fcbd8a88b53d4f3172d Mon Sep 17 00:00:00 2001
From: Colm Talbot <colm.talbot@ligo.org>
Date: Tue, 28 May 2019 00:36:05 -0500
Subject: [PATCH] Generalise ROQ

---
 bilby/gw/conversion.py                     |  19 ++
 bilby/gw/source.py                         | 270 +++++++++++----------
 bilby/gw/utils.py                          | 141 ++++++++---
 examples/injection_examples/roq_example.py |  11 +-
 test/gw_utils_test.py                      |  23 +-
 5 files changed, 292 insertions(+), 172 deletions(-)

diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py
index 0dea42aa3..c2b842e83 100644
--- a/bilby/gw/conversion.py
+++ b/bilby/gw/conversion.py
@@ -50,6 +50,25 @@ def luminosity_distance_to_comoving_distance(distance, cosmology=None):
     return redshift_to_comoving_distance(redshift, cosmology)
 
 
+def bilby_to_lalsimulation_spins(
+        theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1, mass_2,
+        reference_frequency, phase):
+    if tilt_1 in [0, np.pi] and tilt_2 in [0, np.pi]:
+        spin_1x = 0
+        spin_1y = 0
+        spin_1z = a_1 * np.cos(tilt_1)
+        spin_2x = 0
+        spin_2y = 0
+        spin_2z = a_2 * np.cos(tilt_2)
+        iota = theta_jn
+    else:
+        iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = \
+            transform_precessing_spins(
+                theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1,
+                mass_2, reference_frequency, phase)
+    return iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z
+
+
 @np.vectorize
 def transform_precessing_spins(theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1,
                                a_2, mass_1, mass_2, reference_frequency, phase):
diff --git a/bilby/gw/source.py b/bilby/gw/source.py
index bdb5fcf56..5a90b5be2 100644
--- a/bilby/gw/source.py
+++ b/bilby/gw/source.py
@@ -4,14 +4,13 @@ import numpy as np
 
 from ..core import utils
 from ..core.utils import logger
-from .conversion import transform_precessing_spins
+from .conversion import bilby_to_lalsimulation_spins
 from .utils import (lalsim_GetApproximantFromString,
                     lalsim_SimInspiralFD,
                     lalsim_SimInspiralChooseFDWaveform,
                     lalsim_SimInspiralWaveformParamsInsertTidalLambda1,
                     lalsim_SimInspiralWaveformParamsInsertTidalLambda2,
-                    lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame,
-                    lalsim_SimIMRPhenomPFrequencySequence)
+                    lalsim_SimInspiralChooseFDWaveformSequence)
 
 try:
     import lal
@@ -227,10 +226,9 @@ def _base_lal_cbc_fd_waveform(
 
     approximant = lalsim_GetApproximantFromString(waveform_approximant)
 
-    if (pn_amplitude_order != 0):
-        start_frequency = lalsim.SimInspiralfLow2fStart(minimum_frequency,
-                                                        int(pn_amplitude_order),
-                                                        approximant)
+    if pn_amplitude_order != 0:
+        start_frequency = lalsim.SimInspiralfLow2fStart(
+            minimum_frequency, int(pn_amplitude_order), approximant)
     else:
         start_frequency = minimum_frequency
 
@@ -243,30 +241,27 @@ def _base_lal_cbc_fd_waveform(
     mass_1 = mass_1 * utils.solar_mass
     mass_2 = mass_2 * utils.solar_mass
 
-    if tilt_1 in [0.0, np.pi] and tilt_2 in [0, np.pi]:
-        spin_1x = 0.0
-        spin_1y = 0.0
-        spin_1z = a_1 * np.cos(tilt_1)
-        spin_2x = 0.0
-        spin_2y = 0.0
-        spin_2z = a_2 * np.cos(tilt_2)
-        iota = theta_jn
-    else:
-        iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = (
-            transform_precessing_spins(
-                theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1,
-                mass_2, reference_frequency, phase))
+    iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = bilby_to_lalsimulation_spins(
+        theta_jn=theta_jn, phi_jl=phi_jl, tilt_1=tilt_1, tilt_2=tilt_2,
+        phi_12=phi_12, a_1=a_1, a_2=a_2, mass_1=mass_1, mass_2=mass_2,
+        reference_frequency=reference_frequency, phase=phase)
 
     longitude_ascending_nodes = 0.0
     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)
+    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)
 
     if lalsim.SimInspiralImplementedFDApproximants(approximant):
         wf_func = lalsim_SimInspiralChooseFDWaveform
@@ -291,70 +286,53 @@ def _base_lal_cbc_fd_waveform(
     h_plus *= frequency_bounds
     h_cross *= frequency_bounds
 
-    return {'plus': h_plus, 'cross': h_cross}
-
-
-def sinegaussian(frequency_array, hrss, Q, frequency, **kwargs):
-    tau = Q / (np.sqrt(2.0) * np.pi * frequency)
-    temp = Q / (4.0 * np.sqrt(np.pi) * frequency)
-    fm = frequency_array - frequency
-    fp = frequency_array + frequency
-
-    h_plus = ((hrss / np.sqrt(temp * (1 + np.exp(-Q**2)))) *
-              ((np.sqrt(np.pi) * tau) / 2.0) *
-              (np.exp(-fm**2 * np.pi**2 * tau**2) +
-              np.exp(-fp**2 * np.pi**2 * tau**2)))
-
-    h_cross = (-1j * (hrss / np.sqrt(temp * (1 - np.exp(-Q**2)))) *
-               ((np.sqrt(np.pi) * tau) / 2.0) *
-               (np.exp(-fm**2 * np.pi**2 * tau**2) -
-               np.exp(-fp**2 * np.pi**2 * tau**2)))
-
-    return{'plus': h_plus, 'cross': h_cross}
-
-
-def supernova(
-        frequency_array, realPCs, imagPCs, file_path, luminosity_distance, **kwargs):
-    """ A supernova NR simulation for injections """
-
-    realhplus, imaghplus, realhcross, imaghcross = np.loadtxt(
-        file_path, usecols=(0, 1, 2, 3), unpack=True)
-
-    # waveform in file at 10kpc
-    scaling = 1e-3 * (10.0 / luminosity_distance)
-
-    h_plus = scaling * (realhplus + 1.0j * imaghplus)
-    h_cross = scaling * (realhcross + 1.0j * imaghcross)
-    return {'plus': h_plus, 'cross': h_cross}
-
+    return dict(plus=h_plus, cross=h_cross)
 
-def supernova_pca_model(
-        frequency_array, pc_coeff1, pc_coeff2, pc_coeff3, pc_coeff4, pc_coeff5,
-        luminosity_distance, **kwargs):
-    """ Supernova signal model """
 
-    realPCs = kwargs['realPCs']
-    imagPCs = kwargs['imagPCs']
+def roq(
+        frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
+        phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **waveform_arguments):
+    logger.warning(
+        'bilby.gw.source.roq has been deprecated and will be removed. Use '
+        'bilby.gw.source.binary_black_hole_roq instead.')
+    return binary_black_hole_roq(
+        frequency_array=frequency_array, mass_1=mass_1, mass_2=mass_2,
+        luminosity_distance=luminosity_distance, a_1=a_1, tilt_1=tilt_1,
+        phi_12=phi_12, a_2=a_2, tilt_2=tilt_2, phi_jl=phi_jl,
+        theta_jn=theta_jn, phase=phase, **waveform_arguments)
 
-    pc1 = realPCs[:, 0] + 1.0j * imagPCs[:, 0]
-    pc2 = realPCs[:, 1] + 1.0j * imagPCs[:, 1]
-    pc3 = realPCs[:, 2] + 1.0j * imagPCs[:, 2]
-    pc4 = realPCs[:, 3] + 1.0j * imagPCs[:, 3]
-    pc5 = realPCs[:, 4] + 1.0j * imagPCs[:, 5]
 
-    # file at 10kpc
-    scaling = 1e-23 * (10.0 / luminosity_distance)
+def binary_black_hole_roq(
+        frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
+        phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **waveform_arguments):
+    waveform_kwargs = dict(
+        waveform_approximant='IMRPhenomPv2', reference_frequency=20.0)
+    waveform_kwargs.update(waveform_arguments)
+    return _base_roq_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, phi_jl=phi_jl,
+        phi_12=phi_12, lambda_1=0.0, lambda_2=0.0, **waveform_kwargs)
 
-    h_plus = scaling * (pc_coeff1 * pc1 + pc_coeff2 * pc2 + pc_coeff3 * pc3 +
-                        pc_coeff4 * pc4 + pc_coeff5 * pc5)
-    h_cross = scaling * (pc_coeff1 * pc1 + pc_coeff2 * pc2 + pc_coeff3 * pc3 +
-                         pc_coeff4 * pc4 + pc_coeff5 * pc5)
 
-    return {'plus': h_plus, 'cross': h_cross}
+def binary_neutron_star_roq(
+        frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
+        phi_12, a_2, tilt_2, phi_jl, lambda_1, lambda_2, theta_jn, phase,
+        **waveform_arguments):
+    waveform_kwargs = dict(
+        waveform_approximant='IMRPhenomD_NRTidal', reference_frequency=20.0)
+    waveform_kwargs.update(waveform_arguments)
+    return _base_roq_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, phi_jl=phi_jl,
+        phi_12=phi_12, lambda_1=lambda_1, lambda_2=lambda_2, **waveform_kwargs)
 
 
-def roq(frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
-        phi_12, a_2, tilt_2, phi_jl, theta_jn, phase, **waveform_arguments):
+def _base_roq_waveform(
+        frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
+        phi_12, a_2, tilt_2, lambda_1, lambda_2, phi_jl, theta_jn, phase,
+        **waveform_arguments):
     """
     See https://git.ligo.org/lscsoft/lalsuite/blob/master/lalsimulation/src/LALSimInspiral.c#L1460
 
@@ -391,7 +369,7 @@ def roq(frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
     frequency_nodes_linear: np.array
     frequency_nodes_quadratic: np.array
     reference_frequency: float
-    version: str
+    approximant: str
 
     Note: for the frequency_nodes_linear and frequency_nodes_quadratic arguments,
     if using data from https://git.ligo.org/lscsoft/ROQ_data, this should be
@@ -402,59 +380,101 @@ def roq(frequency_array, mass_1, mass_2, luminosity_distance, a_1, tilt_1,
     waveform_polarizations: dict
         Dict containing plus and cross modes evaluated at the linear and
         quadratic frequency nodes.
-
     """
     frequency_nodes_linear = waveform_arguments['frequency_nodes_linear']
     frequency_nodes_quadratic = waveform_arguments['frequency_nodes_quadratic']
-    reference_frequency = getattr(waveform_arguments,
-                                  'reference_frequency', 20.0)
-    versions = dict(IMRPhenomPv2=lalsim.IMRPhenomPv2_V)
-    version = versions[getattr(waveform_arguments, 'version', 'IMRPhenomPv2')]
+    reference_frequency = waveform_arguments['reference_frequency']
+    approximant = lalsim_GetApproximantFromString(
+        waveform_arguments['approximant'])
 
     luminosity_distance = luminosity_distance * 1e6 * utils.parsec
     mass_1 = mass_1 * utils.solar_mass
     mass_2 = mass_2 * utils.solar_mass
 
-    if tilt_1 == 0 and tilt_2 == 0:
-        spin_1x = 0
-        spin_1y = 0
-        spin_1z = a_1
-        spin_2x = 0
-        spin_2y = 0
-        spin_2z = a_2
-        iota = theta_jn
-    else:
-        iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = \
-            transform_precessing_spins(
-                theta_jn, phi_jl, tilt_1, tilt_2, phi_12, a_1, a_2, mass_1,
-                mass_2, reference_frequency, phase)
-
-    chi_1_l, chi_2_l, chi_p, theta_jn, alpha, phase_aligned, zeta =\
-        lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame(
-            mass_1, mass_2, reference_frequency, phase, iota, spin_1x,
-            spin_1y, spin_1z, spin_2x, spin_2y, spin_2z, version)
+    waveform_dictionary = lal.CreateDict()
+    lalsim_SimInspiralWaveformParamsInsertTidalLambda1(
+        waveform_dictionary, lambda_1)
+    lalsim_SimInspiralWaveformParamsInsertTidalLambda2(
+        waveform_dictionary, lambda_2)
+
+    iota, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z = bilby_to_lalsimulation_spins(
+        theta_jn=theta_jn, phi_jl=phi_jl, tilt_1=tilt_1, tilt_2=tilt_2,
+        phi_12=phi_12, a_1=a_1, a_2=a_2, mass_1=mass_1, mass_2=mass_2,
+        reference_frequency=reference_frequency, phase=phase)
+
+    h_linear_plus, h_linear_cross = lalsim_SimInspiralChooseFDWaveformSequence(
+        phase, mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y,
+        spin_2z, reference_frequency, luminosity_distance, iota,
+        waveform_dictionary, approximant, frequency_nodes_linear)
+
+    h_quadratic_plus, h_quadratic_cross = lalsim_SimInspiralChooseFDWaveformSequence(
+        phase, mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y,
+        spin_2z, reference_frequency, luminosity_distance, iota,
+        waveform_dictionary, approximant, frequency_nodes_quadratic)
 
     waveform_polarizations = dict()
-
-    h_linear_plus, h_linear_cross = lalsim_SimIMRPhenomPFrequencySequence(
-        frequency_nodes_linear, chi_1_l, chi_2_l, chi_p, theta_jn,
-        mass_1, mass_2, luminosity_distance,
-        alpha, phase_aligned, reference_frequency, version)
-    h_quadratic_plus, h_quadratic_cross = lalsim_SimIMRPhenomPFrequencySequence(
-        frequency_nodes_quadratic, chi_1_l, chi_2_l, chi_p, theta_jn,
-        mass_1, mass_2, luminosity_distance,
-        alpha, phase_aligned, reference_frequency, version)
-
     waveform_polarizations['linear'] = dict(
-        plus=(np.cos(2 * zeta) * h_linear_plus.data.data +
-              np.sin(2 * zeta) * h_linear_cross.data.data),
-        cross=(np.cos(2 * zeta) * h_linear_cross.data.data -
-               np.sin(2 * zeta) * h_linear_plus.data.data))
-
+        plus=h_linear_plus.data.data, cross=h_linear_cross.data.data)
     waveform_polarizations['quadratic'] = dict(
-        plus=(np.cos(2 * zeta) * h_quadratic_plus.data.data +
-              np.sin(2 * zeta) * h_quadratic_cross.data.data),
-        cross=(np.cos(2 * zeta) * h_quadratic_cross.data.data -
-               np.sin(2 * zeta) * h_quadratic_plus.data.data))
+        plus=h_quadratic_plus.data.data, cross=h_quadratic_cross.data.data)
 
     return waveform_polarizations
+
+
+def sinegaussian(frequency_array, hrss, Q, frequency, **kwargs):
+    tau = Q / (np.sqrt(2.0) * np.pi * frequency)
+    temp = Q / (4.0 * np.sqrt(np.pi) * frequency)
+    fm = frequency_array - frequency
+    fp = frequency_array + frequency
+
+    h_plus = ((hrss / np.sqrt(temp * (1 + np.exp(-Q**2)))) *
+              ((np.sqrt(np.pi) * tau) / 2.0) *
+              (np.exp(-fm**2 * np.pi**2 * tau**2) +
+              np.exp(-fp**2 * np.pi**2 * tau**2)))
+
+    h_cross = (-1j * (hrss / np.sqrt(temp * (1 - np.exp(-Q**2)))) *
+               ((np.sqrt(np.pi) * tau) / 2.0) *
+               (np.exp(-fm**2 * np.pi**2 * tau**2) -
+               np.exp(-fp**2 * np.pi**2 * tau**2)))
+
+    return{'plus': h_plus, 'cross': h_cross}
+
+
+def supernova(
+        frequency_array, realPCs, imagPCs, file_path, luminosity_distance, **kwargs):
+    """ A supernova NR simulation for injections """
+
+    realhplus, imaghplus, realhcross, imaghcross = np.loadtxt(
+        file_path, usecols=(0, 1, 2, 3), unpack=True)
+
+    # waveform in file at 10kpc
+    scaling = 1e-3 * (10.0 / luminosity_distance)
+
+    h_plus = scaling * (realhplus + 1.0j * imaghplus)
+    h_cross = scaling * (realhcross + 1.0j * imaghcross)
+    return {'plus': h_plus, 'cross': h_cross}
+
+
+def supernova_pca_model(
+        frequency_array, pc_coeff1, pc_coeff2, pc_coeff3, pc_coeff4, pc_coeff5,
+        luminosity_distance, **kwargs):
+    """ Supernova signal model """
+
+    realPCs = kwargs['realPCs']
+    imagPCs = kwargs['imagPCs']
+
+    pc1 = realPCs[:, 0] + 1.0j * imagPCs[:, 0]
+    pc2 = realPCs[:, 1] + 1.0j * imagPCs[:, 1]
+    pc3 = realPCs[:, 2] + 1.0j * imagPCs[:, 2]
+    pc4 = realPCs[:, 3] + 1.0j * imagPCs[:, 3]
+    pc5 = realPCs[:, 4] + 1.0j * imagPCs[:, 5]
+
+    # file at 10kpc
+    scaling = 1e-23 * (10.0 / luminosity_distance)
+
+    h_plus = scaling * (pc_coeff1 * pc1 + pc_coeff2 * pc2 + pc_coeff3 * pc3 +
+                        pc_coeff4 * pc4 + pc_coeff5 * pc5)
+    h_cross = scaling * (pc_coeff1 * pc1 + pc_coeff2 * pc2 + pc_coeff3 * pc3 +
+                         pc_coeff4 * pc4 + pc_coeff5 * pc5)
+
+    return {'plus': h_plus, 'cross': h_cross}
diff --git a/bilby/gw/utils.py b/bilby/gw/utils.py
index d7b56d453..16937ebae 100644
--- a/bilby/gw/utils.py
+++ b/bilby/gw/utils.py
@@ -18,6 +18,7 @@ except ImportError:
                    " not be able to use some of the prebuilt functions.")
 
 try:
+    import lal
     import lalsimulation as lalsim
 except ImportError:
     logger.warning("You do not have lalsuite installed currently. You will"
@@ -688,8 +689,33 @@ def lalsim_SimInspiralFD(
         longitude_ascending_nodes, eccentricity, mean_per_ano, delta_frequency,
         minimum_frequency, maximum_frequency, reference_frequency,
         waveform_dictionary, approximant):
+    """
+    Safely call lalsimulation.SimInspiralFD
+
+    Parameters
+    ----------
+    phase: float, int
+    mass_1: float, int
+    mass_2: float, int
+    spin_1x: float, int
+    spin_1y: float, int
+    spin_1z: float, int
+    spin_2x: float, int
+    spin_2y: float, int
+    spin_2z: float, int
+    reference_frequency: float, int
+    luminosity_distance: float, int
+    iota: float, int
+    longitude_ascending_nodes: float, int
+    eccentricity: float, int
+    mean_per_ano: float, int
+    delta_frequency: float, int
+    minimum_frequency: float, int
+    maximum_frequency: float, int
+    waveform_dictionary: None, lal.Dict
+    approximant: int, str
+    """
 
-    # Convert values to floats
     [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,
@@ -699,8 +725,11 @@ def lalsim_SimInspiralFD(
         eccentricity, mean_per_ano, delta_frequency, minimum_frequency,
         maximum_frequency, reference_frequency)
 
-    # Note, this is the approximant number returns by GetApproximantFromString
-    if isinstance(approximant, int) is False:
+    if isinstance(approximant, int):
+        pass
+    elif isinstance(approximant, str):
+        approximant = lalsim_GetApproximantFromString(approximant)
+    else:
         raise ValueError("approximant not an int")
 
     return lalsim.SimInspiralFD(
@@ -717,8 +746,33 @@ def lalsim_SimInspiralChooseFDWaveform(
         longitude_ascending_nodes, eccentricity, mean_per_ano, delta_frequency,
         minimum_frequency, maximum_frequency, reference_frequency,
         waveform_dictionary, approximant):
+    """
+    Safely call lalsimulation.SimInspiralChooseFDWaveform
+
+    Parameters
+    ----------
+    phase: float, int
+    mass_1: float, int
+    mass_2: float, int
+    spin_1x: float, int
+    spin_1y: float, int
+    spin_1z: float, int
+    spin_2x: float, int
+    spin_2y: float, int
+    spin_2z: float, int
+    reference_frequency: float, int
+    luminosity_distance: float, int
+    iota: float, int
+    longitude_ascending_nodes: float, int
+    eccentricity: float, int
+    mean_per_ano: float, int
+    delta_frequency: float, int
+    minimum_frequency: float, int
+    maximum_frequency: float, int
+    waveform_dictionary: None, lal.Dict
+    approximant: int, str
+    """
 
-    # Convert values to floats
     [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,
@@ -728,8 +782,11 @@ def lalsim_SimInspiralChooseFDWaveform(
         eccentricity, mean_per_ano, delta_frequency, minimum_frequency,
         maximum_frequency, reference_frequency)
 
-    # Note, this is the approximant number returns by GetApproximantFromString
-    if isinstance(approximant, int) is False:
+    if isinstance(approximant, int):
+        pass
+    elif isinstance(approximant, str):
+        approximant = lalsim_GetApproximantFromString(approximant)
+    else:
         raise ValueError("approximant not an int")
 
     return lalsim.SimInspiralChooseFDWaveform(
@@ -740,31 +797,53 @@ def lalsim_SimInspiralChooseFDWaveform(
         waveform_dictionary, approximant)
 
 
-def lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame(
-        mass_1, mass_2, reference_frequency, phase, iota, spin_1x,
-        spin_1y, spin_1z, spin_2x, spin_2y, spin_2z, version):
-    [mass_1, mass_2, reference_frequency, phase, iota, spin_1x,
-     spin_1y, spin_1z, spin_2x, spin_2y, spin_2z] = convert_args_list_to_float(
-        mass_1, mass_2, reference_frequency, phase, iota, spin_1x,
-        spin_1y, spin_1z, spin_2x, spin_2y, spin_2z)
-    return lalsim.SimIMRPhenomPCalculateModelParametersFromSourceFrame(
-        mass_1, mass_2, reference_frequency, phase, iota, spin_1x,
-        spin_1y, spin_1z, spin_2x, spin_2y, spin_2z, version)
-
-
-def lalsim_SimIMRPhenomPFrequencySequence(
-        frequency_nodes, chi_1_l, chi_2_l, chi_p, theta_jn,
-        mass_1, mass_2, luminosity_distance,
-        alpha, phase_aligned, reference_frequency, version):
-    [chi_1_l, chi_2_l, chi_p, theta_jn, mass_1, mass_2, luminosity_distance,
-     alpha, phase_aligned, reference_frequency] = convert_args_list_to_float(
-        chi_1_l, chi_2_l, chi_p, theta_jn, mass_1, mass_2, luminosity_distance,
-        alpha, phase_aligned, reference_frequency)
-
-    return lalsim.SimIMRPhenomPFrequencySequence(
-        frequency_nodes, chi_1_l, chi_2_l, chi_p, theta_jn, mass_1, mass_2,
-        luminosity_distance, alpha, phase_aligned, reference_frequency, version,
-        None)
+def lalsim_SimInspiralChooseFDWaveformSequence(
+        phase, mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y,
+        spin_2z, reference_frequency, luminosity_distance, iota,
+        waveform_dictionary, approximant, frequency_array):
+    """
+    Safely call lalsimulation.SimInspiralChooseFDWaveformSequence
+
+    Parameters
+    ----------
+    phase: float, int
+    mass_1: float, int
+    mass_2: float, int
+    spin_1x: float, int
+    spin_1y: float, int
+    spin_1z: float, int
+    spin_2x: float, int
+    spin_2y: float, int
+    spin_2z: float, int
+    reference_frequency: float, int
+    luminosity_distance: float, int
+    iota: float, int
+    waveform_dictionary: None, lal.Dict
+    approximant: int, str
+    frequency_array: np.ndarray, lal.REAL8Vector
+    """
+
+    [mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z,
+     luminosity_distance, iota, phase, reference_frequency] = convert_args_list_to_float(
+        mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y, spin_2z,
+        luminosity_distance, iota, phase, reference_frequency)
+
+    if isinstance(approximant, int):
+        pass
+    elif isinstance(approximant, str):
+        approximant = lalsim_GetApproximantFromString(approximant)
+    else:
+        raise ValueError("approximant not an int")
+
+    if not isinstance(frequency_array, lal.REAL8Vector):
+        old_frequency_array = frequency_array.copy()
+        frequency_array = lal.CreateREAL8Vector(len(old_frequency_array))
+        frequency_array.data = old_frequency_array
+
+    return lalsim.SimInspiralChooseFDWaveformSequence(
+        phase, mass_1, mass_2, spin_1x, spin_1y, spin_1z, spin_2x, spin_2y,
+        spin_2z, reference_frequency, luminosity_distance, iota,
+        waveform_dictionary, approximant, frequency_array)
 
 
 def lalsim_SimInspiralWaveformParamsInsertTidalLambda1(
diff --git a/examples/injection_examples/roq_example.py b/examples/injection_examples/roq_example.py
index b98f4f2bf..cf6a28250 100644
--- a/examples/injection_examples/roq_example.py
+++ b/examples/injection_examples/roq_example.py
@@ -69,12 +69,11 @@ ifos.inject_signal(waveform_generator=waveform_generator,
 # make ROQ waveform generator
 search_waveform_generator = bilby.gw.waveform_generator.WaveformGenerator(
     duration=duration, sampling_frequency=sampling_frequency,
-    frequency_domain_source_model=bilby.gw.source.roq,
-    waveform_arguments=dict(frequency_nodes_linear=freq_nodes_linear,
-                            frequency_nodes_quadratic=freq_nodes_quadratic,
-                            reference_frequency=20. * scale_factor,
-                            minimum_frequency=20. * scale_factor,
-                            approximant='IMRPhenomPv2'),
+    frequency_domain_source_model=bilby.gw.source.binary_black_hole_roq,
+    waveform_arguments=dict(
+        frequency_nodes_linear=freq_nodes_linear,
+        frequency_nodes_quadratic=freq_nodes_quadratic,
+        reference_frequency=20. * scale_factor, approximant='IMRPhenomPv2'),
     parameter_conversion=bilby.gw.conversion.convert_to_lal_binary_black_hole_parameters)
 
 # Here we add constraints on chirp mass and mass ratio to the prior, these are
diff --git a/test/gw_utils_test.py b/test/gw_utils_test.py
index 90e67fd0d..b8d1838ca 100644
--- a/test/gw_utils_test.py
+++ b/test/gw_utils_test.py
@@ -38,6 +38,10 @@ class TestGWUtils(unittest.TestCase):
         self.assertTrue(np.all(psd == (freq_data * 2 * df ** 0.5)**2))
 
     def test_time_delay_from_geocenter(self):
+        """
+        The difference in the two detector case is due to rounding error.
+        Different hardware gives different numbers in the last decimal place.
+        """
         det1 = np.array([0.1, 0.2, 0.3])
         det2 = np.array([0.1, 0.2, 0.5])
         ra = 0.5
@@ -45,9 +49,9 @@ class TestGWUtils(unittest.TestCase):
         time = 10
         self.assertEqual(
             gwutils.time_delay_geocentric(det1, det1, ra, dec, time), 0)
-        self.assertEqual(
+        self.assertAlmostEqual(
             gwutils.time_delay_geocentric(det1, det2, ra, dec, time),
-            1.3253791114055397e-10)
+            1.3253791114055397e-10, 14)
 
     def test_get_polarization_tensor(self):
         ra = 1
@@ -171,16 +175,15 @@ class TestGWUtils(unittest.TestCase):
         self.assertEqual(type(a[0]), lal.COMPLEX16FrequencySeries)
         self.assertEqual(type(a[1]), lal.COMPLEX16FrequencySeries)
 
-        with self.assertRaises(ValueError):
-            a = gwutils.lalsim_SimInspiralChooseFDWaveform(
+        with self.assertRaises(RuntimeError):
+            _ = gwutils.lalsim_SimInspiralChooseFDWaveform(
                 35.2, 20.4, 0.1, 0.2, 0.2, 0.2, 0.2, 0.1, 1000, 2, 2.3,
-                45., 0.1, 10, 0.01, 10, 1000, 20, None, 'IMRPhenomPV2')
+                45., 0.1, 10, 0.01, 10, 1000, 20, None, 'Fail')
 
-    def test_lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame(self):
-        version = lalsim.IMRPhenomPv2_V
-        a = gwutils.lalsim_SimIMRPhenomPCalculateModelParametersFromSourceFrame(
-            25.6, 12.2, 50, 0.2, 0, 0, 0, 0, 0, 0, 0, version)
-        self.assertEqual(len(a), 7)
+        with self.assertRaises(ValueError):
+            _ = gwutils.lalsim_SimInspiralChooseFDWaveform(
+                35.2, 20.4, 0.1, 0.2, 0.2, 0.2, 0.2, 0.1, 1000, 2, 2.3,
+                45., 0.1, 10, 0.01, 10, 1000, 20, None, 1.5)
 
 
 if __name__ == '__main__':
-- 
GitLab