Using gwsignal
import lalsimulation.gwsignal as gws
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import copy
from phenomxpy.gwsignal_wrapper import *
/home/cecilio.garcia-quiros/.conda/envs/phenomxpy/lib/python3.10/site-packages/lalsimulation/lalsimulation.py:8: UserWarning: Wswiglal-redir-stdio:
SWIGLAL standard output/error redirection is enabled in IPython.
This may lead to performance penalties. To disable locally, use:
with lal.no_swig_redirect_standard_output_error():
...
To disable globally, use:
lal.swig_redirect_standard_output_error(False)
Note however that this will likely lead to error messages from
LAL functions being either misdirected or lost when called from
Jupyter notebooks.
To suppress this warning, use:
import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")
import lal
import lal
/home/cecilio.garcia-quiros/.conda/envs/phenomxpy/lib/python3.10/site-packages/cupyx/jit/_interface.py:173: FutureWarning: cupyx.jit.rawkernel is experimental. The interface can change in the future.
cupy._util.experimental('cupyx.jit.rawkernel')
# start with the usual parameter definitions
m1 = 30.0 * u.solMass
m2 = 10.0 * u.solMass
s1x = 0.5 * u.dimensionless_unscaled
s1y = 0.2 * u.dimensionless_unscaled
s1z = 0.0 * u.dimensionless_unscaled
s2x = 0.5 * u.dimensionless_unscaled
s2y = -0.2 * u.dimensionless_unscaled
s2z = 0.0 * u.dimensionless_unscaled
deltaT = 1.0 / 8048.0 * u.s
f_min = 20.0 * u.Hz
f_ref = 20.0 * u.Hz
distance = 1000.0 * u.Mpc
inclination = 1.5 * u.rad
phiRef = 0.0 * u.rad
eccentricity = 0.0 * u.dimensionless_unscaled
longAscNodes = 0.0 * u.rad
meanPerAno = 0.0 * u.rad
params = {
"mass1": m1,
"mass2": m2,
"spin1x": s1x,
"spin1y": s1y,
"spin1z": s1z,
"spin2x": s2x,
"spin2y": s2y,
"spin2z": s2z,
"deltaT": deltaT,
"f22_start": f_min,
"f22_ref": f_ref,
"phi_ref": phiRef,
"distance": distance,
"inclination": inclination,
"eccentricity": eccentricity,
"longAscNodes": longAscNodes,
"meanPerAno": meanPerAno,
#'approximant': 'IMRPhenomTPHM',
"condition": 1,
}
Standard method, one generator per approximant
gen = PyIMRPhenomTPHM() # PyIMRPhenomTHM() PyIMRPhenomTP(), PyIMRPhenomTPHM()
gen.metadata
{'type': 'precessing',
'f_ref_spin': True,
'modes': True,
'polarizations': True,
'implemented_domain': 'time',
'approximant': 'IMRPhenomTPHM',
'implementation': '',
'conditioning_routines': ''}
hp, hc = gws.core.waveform.GenerateTDWaveform(params, gen)
plt.plot(hp, label=gen.metadata["approximant"])
plt.xlabel(hp.times.unit)
plt.ylabel(hp.unit)
plt.legend()
plt.show()
hp, hc = gws.core.waveform.GenerateFDWaveform(params, gen)
plt.loglog(np.abs(hp), label=gen.metadata["approximant"])
plt.xlabel(hp.frequencies.unit)
plt.ylabel(hp.unit)
plt.legend()
plt.show()
Non-default model
params_nondefault = copy.copy(params)
hp, hc = gws.core.waveform.GenerateFDWaveform(params_nondefault, gen)
plt.loglog(np.abs(hp), label="Default (Numerical)")
params_nondefault["prec_version"] = "NNLO"
hp, hc = gws.core.waveform.GenerateFDWaveform(params_nondefault, gen)
plt.loglog(np.abs(hp), label="NNLO")
params_nondefault["final_spin_version"] = "1"
hp, hc = gws.core.waveform.GenerateFDWaveform(params_nondefault, gen)
plt.loglog(np.abs(hp), label="NNLO FS1")
plt.xlabel(hp.frequencies.unit)
plt.ylabel(hp.unit)
plt.legend()
plt.show()
Using a generator for the whole family
gen = PyIMRPhenomTFamily()
gen.metadata
{'type': 'aligned_spin, precessing',
'f_ref_spin': True,
'modes': True,
'polarizations': True,
'implemented_domain': 'time',
'approximant': 'IMRPhenomTFamily',
'implementation': '',
'conditioning_routines': ''}
for approximant in ["IMRPhenomT", "IMRPhenomTHM", "IMRPhenomTP", "IMRPhenomTPHM"]:
params["approximant"] = approximant
hp, hc = gws.core.waveform.GenerateFDWaveform(params, gen)
plt.loglog(np.abs(hp), label=approximant)
plt.legend()
plt.xlabel(hp.frequencies.unit)
plt.ylabel(hp.unit)
plt.show()
Using non-default model
Extra options can be passed to the params dictionary as before, but this generator also supports passing options in the approximant name
for approximant in ["IMRPhenomTPHM_nnlo", "IMRPhenomTPHM_MSA", "IMRPhenomTPHM_SpinTaylor"]:
params["approximant"] = approximant
hp, hc = gws.core.waveform.GenerateTDWaveform(params, gen)
plt.plot(hp, label=approximant)
del params["approximant"]
plt.legend()
plt.xlabel(hp.times.unit)
plt.ylabel(hp.unit)
plt.xlim(-0.1, 0.025)
plt.show()
Comparison with LAL
from lalsimulation.gwsignal import gwsignal_get_waveform_generator
fig, axs = plt.subplots(2, 2, figsize=(10, 6))
params_as = copy.copy(params)
params_as["spin1x"], params_as["spin1y"], params_as["spin2x"], params_as["spin2y"] = (
np.array([0, 0, 0, 0]) * u.dimensionless_unscaled
)
gen_lal = gwsignal_get_waveform_generator("IMRPhenomT")
gen_py = PyIMRPhenomT()
hp, hc = gws.core.waveform.GenerateTDWaveform(params_as, gen_lal)
axs[0, 0].plot(np.real(hp), label="IMRPhenomT LAL")
hp, hc = gws.core.waveform.GenerateTDWaveform(params_as, gen_py)
axs[0, 0].plot(np.real(hp), "--", label="IMtRPhenomT python")
axs[0, 0].legend()
gen_lal = gwsignal_get_waveform_generator("IMRPhenomTHM")
gen_py = PyIMRPhenomTHM()
hp, hc = gws.core.waveform.GenerateTDWaveform(params_as, gen_lal)
axs[0, 1].plot(np.real(hp), label="IMRPhenomTHM LAL")
hp, hc = gws.core.waveform.GenerateTDWaveform(params_as, gen_py)
axs[0, 1].plot(np.real(hp), "--", label="IMRPhenomTHM python")
axs[0, 1].legend()
gen_lal = gwsignal_get_waveform_generator("IMRPhenomTP")
gen_py = PyIMRPhenomTP()
hp, hc = gws.core.waveform.GenerateTDWaveform(params, gen_lal)
axs[1, 0].plot(np.real(hp), label="IMRPhenomTP LAL")
hp, hc = gws.core.waveform.GenerateTDWaveform(params, gen_py)
axs[1, 0].plot(np.real(hp), "--", label="IMRPhenomTP python")
axs[1, 0].legend()
gen_lal = gwsignal_get_waveform_generator("IMRPhenomTPHM")
gen_py = PyIMRPhenomTPHM()
hp, hc = gws.core.waveform.GenerateTDWaveform(params, gen_lal)
axs[1, 1].plot(np.real(hp), label="IMRPhenomTPHM LAL")
hp, hc = gws.core.waveform.GenerateTDWaveform(params, gen_py)
axs[1, 1].plot(np.real(hp), "--", label="IMRPhenomTPHM python")
axs[1, 1].legend()
for ax in axs.flat:
ax.set_xlim(-0.1, 0.05)
plt.show()
/home/cecilio.garcia-quiros/.conda/envs/phenomxpy/lib/python3.10/site-packages/lalsimulation/gwsignal/core/waveform.py:226: UserWarning: This code is currently UNREVIEWED, use with caution!
warnings.warn("This code is currently UNREVIEWED, use with caution!")
fig, axs = plt.subplots(2, 2, figsize=(10, 6))
params["deltaF"] = 0.0625 * u.Hz # NOTE gwsignal uses this value if deltaF not provided
# params["deltaF"] = 0 * u.Hz
params_as = copy.copy(params)
params_as["spin1x"], params_as["spin1y"], params_as["spin2x"], params_as["spin2y"] = (
np.array([0, 0, 0, 0]) * u.dimensionless_unscaled
)
gen_lal = gwsignal_get_waveform_generator("IMRPhenomT")
gen_py = PyIMRPhenomT()
hp, hc = gws.core.waveform.GenerateFDWaveform(params_as, gen_lal)
axs[0, 0].loglog(np.abs(hp), label="IMRPhenomT LAL")
hp, hc = gws.core.waveform.GenerateFDWaveform(params_as, gen_py)
axs[0, 0].loglog(np.abs(hp), label="IMRPhenomT python")
axs[0, 0].legend()
gen_lal = gwsignal_get_waveform_generator("IMRPhenomTHM")
gen_py = PyIMRPhenomTHM()
hp, hc = gws.core.waveform.GenerateFDWaveform(params_as, gen_lal)
axs[0, 1].loglog(np.abs(hp), label="IMRPhenomTHM LAL")
hp, hc = gws.core.waveform.GenerateFDWaveform(params_as, gen_py)
axs[0, 1].loglog(np.abs(hp), label="IMRPhenomTHM python")
axs[0, 1].legend()
gen_lal = gwsignal_get_waveform_generator("IMRPhenomTP")
gen_py = PyIMRPhenomTP()
hp, hc = gws.core.waveform.GenerateFDWaveform(params, gen_lal)
axs[1, 0].loglog(np.abs(hp), label="IMRPhenomTP LAL")
hp, hc = gws.core.waveform.GenerateFDWaveform(params, gen_py)
axs[1, 0].loglog(np.abs(hp), label="IMRPhenomTP python")
axs[1, 0].legend()
gen_lal = gwsignal_get_waveform_generator("IMRPhenomTPHM")
gen_py = PyIMRPhenomTPHM()
hp, hc = gws.core.waveform.GenerateFDWaveform(params, gen_lal)
axs[1, 1].loglog(np.abs(hp), label="IMRPhenomTPHM LAL")
hp, hc = gws.core.waveform.GenerateFDWaveform(params, gen_py)
axs[1, 1].loglog(np.abs(hp), label="IMRPhenomTPHM python")
axs[1, 1].legend()
plt.show()
gen_lal = gwsignal_get_waveform_generator("IMRPhenomTPHM")
gen_py = PyIMRPhenomTPHM()
hp, hc = gws.core.waveform.GenerateFDWaveform(params, gen_lal)
plt.plot(np.real(hp), label="IMRPhenomTPHM LAL")
hp2, hc2 = gws.core.waveform.GenerateFDWaveform(params, gen_py)
plt.plot(np.real(hp2), "--", label="IMRPhenomTPHM python")
plt.legend()
plt.xlim(20, 30)
(20.0, 30.0)
