Commit d678a13e authored by Gregory Ashton's avatar Gregory Ashton

Merge branch 'bns-roq' into 'master'

Generalise ROQ

See merge request !490
parents a74961a8 c7c16d90
Pipeline #64113 passed with stages
in 5 minutes and 44 seconds
......@@ -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):
......
......@@ -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}
......@@ -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(
......
......@@ -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
......
......@@ -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')