SEOBNRv5EHM: aligned spin eccentric model example
This notebook gives examples on how to run the SEOBNRv5EHM waveform generator.
[1]:
import matplotlib.pyplot as plt
import numpy as np
Using the pyseobnr interface
The following demonstrates how to use the internal SEOBNRv5EHM waveform generator. It takes in parameters in geometric units and returns results in geometric units as well. It also gives access to things beyond the waveform modes, such as the dynamics.
[2]:
import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio") # silence LAL warnings
from pyseobnr.generate_waveform import GenerateWaveform, generate_modes_opt
[3]:
# Input parameters
q = 5.3
chi_1 = 0.9
chi_2 = 0.3
omega_start = 0.0137 # This is the orbital frequency in geometric units with M=1
eccentricity = 0.4
rel_anomaly = 2.3
Default mode
We get a time array and a dictionary of complex modes. Note that for aligned spins, only positive m modes are returned
[4]:
t, modes = generate_modes_opt(
q,
chi_1,
chi_2,
omega_start,
eccentricity=eccentricity,
rel_anomaly=rel_anomaly,
approximant="SEOBNRv5EHM",
)
[5]:
modes.keys()
[5]:
dict_keys(['2,2', '2,1', '3,3', '3,2', '4,4', '4,3'])
[6]:
plt.figure()
plt.plot(t, modes["2,2"].real)
plt.xlabel("Time (M)")
plt.ylabel(r"$\Re[h_{22}]$")
plt.grid(True)
[7]:
plt.figure()
plt.plot(t, np.abs(modes["2,2"]))
plt.xlabel("Time (M)")
plt.ylabel(r"$|h_{22}|$")
plt.axvline(x=0, ls="--", color="k")
plt.xlim(-1000, 100)
plt.grid(True)
By default the aligned spin model sets \(t=0\) at the last peak of the frame-invariant amplitude, given a certain threshold
[8]:
for mode in modes.keys():
plt.semilogy(t, np.abs(modes[mode]), label=mode)
plt.xlabel("Time (M)")
plt.ylabel(r"$\|h_{\ell m}|$")
plt.xlim(-500, 150)
plt.ylim(1e-5, 1)
plt.legend(loc=3)
plt.grid(True)
Customising the default mode
Mode array
Suppose you want only the (2,2) and (2,1) modes. This can be done by adding a special argument to the settings dictionary, as follows:
[9]:
settings = dict(
# "EccIC" determines the prescription for the starting frequency.
# EccIC=0 means that omega_start corresponds to the *instantaneous* angular frequency, and
# EccIC=1 means that omega_start corresponds to the *orbit-averaged* angular frequency [default option]
EccIC=1,
# Specify which modes are to be returned,
return_modes=[(2, 2), (2, 1)],
)
# See the pyseobnr/generate_waveform.py file for a description of more settings
t_new, modes_new = generate_modes_opt(
q,
chi_1,
chi_2,
omega_start,
eccentricity=eccentricity,
rel_anomaly=rel_anomaly,
approximant="SEOBNRv5EHM",
settings=settings,
)
[10]:
modes_new.keys()
[10]:
dict_keys(['2,2', '2,1'])
[11]:
plt.figure()
plt.plot(t_new, np.abs(modes_new["2,1"]))
plt.xlabel("Time (M)")
plt.ylabel(r"$|h_{21}|$")
plt.axvline(x=0, ls="--", color="k")
plt.xlim(-1000, 100)
plt.grid(True)
Backwards evolution
One can also evolve the system backwards in time, by specifying a certain amount of time (in geometric units). This will evolve backwards the full equations of motion of the system.
This option is particularly useful for post-processing tasks, e.g. with gw-eccentricity.
This applies to the modes, but also to the BH dynamics (see below Expert mode).
Note that the warning message can be deactivated to avoid filling log files with these messages.
[12]:
t, modes = generate_modes_opt(
q,
chi_1,
chi_2,
omega_start,
eccentricity=eccentricity,
rel_anomaly=rel_anomaly,
approximant="SEOBNRv5EHM",
)
settings = dict(
t_backwards=1000, # Geometric units
# warning_bwd_int=False, # Setting this to False will avoid the warning message
)
t_back, modes_back = generate_modes_opt(
q,
chi_1,
chi_2,
omega_start,
eccentricity=eccentricity,
rel_anomaly=rel_anomaly,
approximant="SEOBNRv5EHM",
settings=settings,
)
WARNING:pyseobnr.models.SEOBNRv5EHM:Integrating backwards in time the full EOB equations of motion by the specified amount of time (|t| = 1000 M) with respect to the reference values. For a total mass of M = 50 M_sun, this corresponds to |t| = 0.24627454738206334 seconds. The starting parameters of the system will change accordingly.
[13]:
plt.figure()
plt.plot(t, np.abs(modes["2,2"]), label="Original system")
plt.plot(t_back, np.abs(modes_back["2,2"]), "--", label="Backwards integration")
plt.xlabel("Time (M)")
plt.ylabel(r"$|h_{2}|$")
plt.axvline(x=0, ls="--", color="k")
plt.legend()
plt.grid(True)
We note that the model has implemented by default a secular backwards integration which is activated whenever the starting separation of the binary is less than \(10M\). This avoids running into issues with the prescription for initial conditions in configurations where the BHs are very close to each other.
In the following, we see an example where this secular backwards integration is activated.
Note that the warning message can be deactivated to avoid filling log files with these messages.
[14]:
settings = dict(
# Setting this to False will avoid the warning message
# warning_secular_bwd_int=False,
)
t, modes = generate_modes_opt(
q,
chi_1,
chi_2,
omega_start,
eccentricity=0.5,
rel_anomaly=0.0,
approximant="SEOBNRv5EHM",
settings=settings,
)
WARNING:pyseobnr.models.SEOBNRv5EHM:The predicted starting separation of the system (r_start = 8.963155505950123) is below the minimum starting separation allowed by the model (r_start_min = 10.0). Waveform parameters are q = 5.3, M = 50, chi_1 = 0.9, chi_2 = 0.3, omega_avg = 0.0137, omega_inst = 0.04177693494914546, eccentricity = 0.5, rel_anomaly = 0.0. Integrating backwards in time a set of secular evolution equations to obtain a prediction for the starting separation larger than this minimum. This will change the starting input values of the system.
WARNING:pyseobnr.models.SEOBNRv5EHM:The new starting input values of the system are: eccentricity = 0.5006713470744493, rel_anomaly = 5.554768213259775, omega_avg = 0.013664469066738846.
[15]:
plt.figure()
plt.plot(t, np.abs(modes["2,2"]))
plt.xlabel("Time (M)")
plt.ylabel(r"$|h_{2}|$")
plt.axvline(x=0, ls="--", color="k")
plt.grid(True)
We can deactivate this secular backwards integration with the following flag: (now, the code will raise an error if the starting separation is below \(10M\))
[16]:
try:
settings = dict(
secular_bwd_int=False,
)
t, modes = generate_modes_opt(
q,
chi_1,
chi_2,
omega_start,
eccentricity=0.5,
rel_anomaly=0.0,
approximant="SEOBNRv5EHM",
settings=settings,
)
except Exception as e:
print(e)
ERROR:pyseobnr.eob.dynamics.initial_conditions_aligned_ecc_opt:Internal function call failed: Input domain error. The predicted post-Newtonian conservative initial separation (r0 = 8.963155505950123 M) is smaller than the minimum separation allowed by the model (r_min = 10.0 M). Aborting the waveform generation since the given input values go outside the validity regime of the model. Please, review the physical sense of the input parameters.
ERROR:pyseobnr.models.SEOBNRv5EHM:Waveform generation failed for q = 5.3, chi_1 = 0.9, chi_2 = 0.3, omega_avg = 0.0137, omega_inst = 0.04177693494914546, eccentricity = 0.5, rel_anomaly = 0.0.
Internal function call failed: Input domain error. The predicted post-Newtonian conservative initial separation (r0 = 8.963155505950123 M) is smaller than the minimum separation allowed by the model (r_min = 10.0 M). Aborting the waveform generation since the given input values go outside the validity regime of the model. Please, review the physical sense of the input parameters.
It is also possible to modify the minimum starting separation. In the following example, we set the minimum starting separation to \(7M\), and hence the secular backward integration is not activated
[17]:
settings = dict(
secular_bwd_int=True,
r_start_min=7.0,
)
t, modes = generate_modes_opt(
q,
chi_1,
chi_2,
omega_start,
eccentricity=0.5,
rel_anomaly=0.0,
approximant="SEOBNRv5EHM",
settings=settings,
)
One can explore combinations of the values of ‘secular_bwd_int’ and ‘r_start_min’. However, one must be careful while doing this, specially for high eccentricities or high frequencies.
However, there is a lower limit of minimum starting separation of \(6M\). Trying to initialize a binary with a separation below \(6M\) will very likely introduce some kind of error. Hence, it is no possible to go below \(6M\) with SEOBNRv5EHM.
Expert mode
In this mode we also get an SEOBNRv5EHM object back which contains more information than just the waveform modes.
[18]:
t, modes, model = generate_modes_opt(
q,
chi_1,
chi_2,
omega_start,
eccentricity=eccentricity,
rel_anomaly=rel_anomaly,
approximant="SEOBNRv5EHM",
debug=True,
)
[19]:
model
[19]:
<pyseobnr.models.SEOBNRv5EHM.SEOBNRv5EHM_opt at 0x7f199e32e1c0>
The model object does of course contain the waveform mode information:
[20]:
plt.figure()
plt.plot(model.t, model.waveform_modes["2,2"].real)
plt.xlabel("Time (M)")
plt.ylabel(r"$\Re[h_{22}]$")
plt.grid(True)
One also has access to the dynamics, stored as \((t, r, \phi, p_r, p_{\phi}, e, \zeta, x, H,\Omega)\), where
\(t\): time
\(r\): relative separation
\(\phi\): azimuthal orbital angle
\(p_r\): tortoise radial momentum
\(p_\phi\): angular momentum
\(e\): eccentricity
\(\zeta\): relativistic anomaly
\(x\): orbit-averaged frequency to the 2/3 power, i.e. \(x=< \omega >^{2/3}\) This variable presents oscillations since we compute it with a PN formula, the true orbit-averaged frequency does not have any oscillation
\(H\): EOB Hamiltonian
\(\Omega\): instantaneous angular frequency, i.e. \(\Omega=\dot\phi\)
[21]:
t, r, phi, pr, pphi, e, z, x, H, Omega = model.dynamics.T
[22]:
plt.figure()
plt.plot(t, r)
plt.xlabel("Time (M)")
plt.ylabel(r"$r$ (M)")
plt.grid(True)
For debugging purposes, other information is also provided, for example, one can examine directly the NQC coefficients:
[23]:
model.nqc_coeffs
[23]:
{(2, 2): {'a1': 73.84005709572764,
'a2': -172.29920201920385,
'a3': -7.882275216720489,
'b1': 47.01109822250898,
'b2': -1690.1118771835474,
'a3S': 0,
'a4': 0,
'a5': 0,
'b3': 0,
'b4': 0},
(2, 1): {'a1': -373.07566322128,
'a2': 3224.773538137499,
'a3': -3651.5331196281195,
'b1': 29.43676013297123,
'b2': -989.7728589584345,
'a3S': 0,
'a4': 0,
'a5': 0,
'b3': 0,
'b4': 0},
(3, 3): {'a1': -1329.7158219985527,
'a2': 13663.695560171982,
'a3': -16341.654099831576,
'b1': 70.97102069452713,
'b2': -3128.784443281338,
'a3S': 0,
'a4': 0,
'a5': 0,
'b3': 0,
'b4': 0},
(3, 2): {'a1': 350.6470921514586,
'a2': -2283.5418735461553,
'a3': 2213.4150852500247,
'b1': 40.069698077754076,
'b2': -1294.4181829522734,
'a3S': 0,
'a4': 0,
'a5': 0,
'b3': 0,
'b4': 0},
(4, 4): {'a1': -190.47818250804391,
'a2': 2932.0191720781963,
'a3': -3835.0912600612173,
'b1': 99.53680829001965,
'b2': -3477.164440899558,
'a3S': 0,
'a4': 0,
'a5': 0,
'b3': 0,
'b4': 0},
(4, 3): {'a1': -729.0930732957175,
'a2': 5640.563241918348,
'a3': -6107.389448194917,
'b1': 69.17144346608978,
'b2': -2388.878987949375,
'a3S': 0,
'a4': 0,
'a5': 0,
'b3': 0,
'b4': 0}}
A note on conventions
The internal SEOBNRv5EHM generator uses the same conventions as in the SEOBNRv5HM model.
Generate modes and polarizations in physical units with LAL conventions
The GenerateWaveform() class accepts a dictionary of parameters (example below) and from it, one can recover the gravitational modes dictionary with the right convention and physical scaling, the time-domain polarizations and the Fourier-domain polarizations
[24]:
# start with the usual parameter definitions
m1 = 50.0
m2 = 30.0
s1x = 0.0
s1y = 0.0
s1z = 0.0
s2x = 0.0
s2y = 0.0
s2z = 0.0
eccentricity = 0.3
rel_anomaly = np.pi
deltaT = 1.0 / 1024.0
f_min = 20.0
f_max = 512.0
distance = 1000.0
inclination = np.pi / 3.0
approximant = "SEOBNRv5EHM"
params_dict = {
"mass1": m1,
"mass2": m2,
"spin1x": s1x,
"spin1y": s1y,
"spin1z": s1z,
"spin2x": s2x,
"spin2y": s2y,
"spin2z": s2z,
"deltaT": deltaT,
"f22_start": f_min,
"distance": distance,
"inclination": inclination,
"f_max": f_max,
"return_modes": [
(2, 2),
(2, 1),
(3, 2),
(3, 3),
(4, 3),
(4, 4),
], # Specify which modes are to be returned
"approximant": approximant,
"eccentricity": eccentricity,
"rel_anomaly": rel_anomaly,
"EccIC": 1, # EccIC = 0 for instantaneous initial orbital frequency, and EccIC = 1 for orbit-averaged initial orbital frequency
}
[25]:
wfm_gen = GenerateWaveform(params_dict) # We call the generator with the parameters
[26]:
wfm_gen.parameters # Access the parameters
[26]:
{'EccIC': 1,
't_backwards': 0.0,
'secular_bwd_int': True,
'warning_bwd_int': True,
'warning_secular_bwd_int': True,
'spin1x': 0.0,
'spin1y': 0.0,
'spin1z': 0.0,
'spin2x': 0.0,
'spin2y': 0.0,
'spin2z': 0.0,
'eccentricity': 0.3,
'rel_anomaly': 3.141592653589793,
'distance': 1000.0,
'inclination': 1.0471975511965976,
'phi_ref': 0.0,
'f22_start': 20.0,
'deltaT': 0.0009765625,
'deltaF': 0.0,
'ModeArray': None,
'mode_array': None,
'approximant': 'SEOBNRv5EHM',
'conditioning': 1,
'polarizations_from_coprec': False,
'initial_conditions': 'adiabatic',
'initial_conditions_postadiabatic_type': 'analytic',
'postadiabatic': False,
'postadiabatic_type': 'analytic',
'r_size_input': 12,
'dA_dict': {'2,2': 0.0,
'2,1': 0.0,
'3,3': 0.0,
'3,2': 0.0,
'4,4': 0.0,
'4,3': 0.0,
'5,5': 0.0},
'dw_dict': {'2,2': 0.0,
'2,1': 0.0,
'3,3': 0.0,
'3,2': 0.0,
'4,4': 0.0,
'4,3': 0.0,
'5,5': 0.0},
'dTpeak': 0.0,
'da6': 0.0,
'ddSO': 0.0,
'domega_dict': {'2,2': 0.0,
'2,1': 0.0,
'3,3': 0.0,
'3,2': 0.0,
'4,4': 0.0,
'4,3': 0.0,
'5,5': 0.0},
'dtau_dict': {'2,2': 0.0,
'2,1': 0.0,
'3,3': 0.0,
'3,2': 0.0,
'4,4': 0.0,
'4,3': 0.0,
'5,5': 0.0},
'tol_PA': 1e-11,
'rtol_ode': 1e-11,
'atol_ode': 1e-12,
'deltaT_sampling': False,
'omega_prec_deviation': True,
'enable_antisymmetric_modes': False,
'antisymmetric_modes': [(2, 2)],
'antisymmetric_modes_hm': False,
'ivs_mrd': None,
'antisymmetric_fits_version': 240417,
'gwsignal_environment': False,
'mass1': 50.0,
'mass2': 30.0,
'f_max': 512.0,
'return_modes': [(2, 2), (2, 1), (3, 2), (3, 3), (4, 3), (4, 4)],
'f_ref': 20.0}
[27]:
# Generate mode dictionary
times, hlm = wfm_gen.generate_td_modes()
plt.figure()
plt.plot(times, hlm[(2, 2)].real)
plt.xlabel("Time (seconds)")
plt.ylabel(r"$\Re[h_{22}]$")
plt.grid(True)
plt.show()
plt.figure()
plt.plot(times, hlm[(3, 3)].imag)
plt.xlabel("Time (seconds)")
plt.ylabel(r"$\Im[h_{33}]$")
plt.grid(True)
plt.show()
[28]:
# Generate time-domain polarizations - As LAL REAL8TimeSeries
hp, hc = wfm_gen.generate_td_polarizations()
times = hp.deltaT * np.arange(hp.data.length) + hp.epoch
plt.figure()
plt.plot(times, hp.data.data, label=r"$h_+$")
plt.plot(times, hc.data.data, label=r"$h_{\times}$")
plt.xlabel("Time (seconds)")
plt.legend()
plt.grid(True)
plt.show()
[29]:
# Generate Fourier-domain polarizations - As LAL COMPLEX16FrequencySeries
hpf, hcf = wfm_gen.generate_fd_polarizations()
freqs = hpf.deltaF * np.arange(hpf.data.length)
plt.figure()
plt.plot(freqs, np.abs(hpf.data.data), label=r"$\tilde{h}_+$")
plt.plot(freqs, np.abs(hcf.data.data), label=r"$\tilde{h}_{\times}$")
plt.xlabel("f (Hz)")
plt.legend()
plt.xscale("log")
plt.yscale("log")
plt.grid(True)
plt.show()
Backwards evolution from GenerateWaveform
We also show here the backwards evolution with the GenerateWaveform class
[30]:
m1 = 50.0
m2 = 30.0
eccentricity = 0.3
rel_anomaly = np.pi
approximant = "SEOBNRv5EHM"
params_dict = {
"mass1": m1,
"mass2": m2,
"approximant": approximant,
"eccentricity": eccentricity,
"rel_anomaly": rel_anomaly,
}
params_dict_back = {
"mass1": m1,
"mass2": m2,
"approximant": approximant,
"eccentricity": eccentricity,
"rel_anomaly": rel_anomaly,
"t_backwards": 1000,
# 'warning_bwd_int' : False, # Setting this to False will avoid the warning message
}
wfm_gen = GenerateWaveform(params_dict) # We call the generator with the parameters
wfm_gen_back = GenerateWaveform(
params_dict_back
) # We call the generator with the parameters
[31]:
# Generate mode dictionary
times, hlm = wfm_gen.generate_td_modes()
times_back, hlm_back = wfm_gen_back.generate_td_modes()
WARNING:pyseobnr.models.SEOBNRv5EHM:Integrating backwards in time the full EOB equations of motion by the specified amount of time (|t| = 1000 M) with respect to the reference values. For a total mass of M = 80.0 M_sun, this corresponds to |t| = 0.3940392758113013 seconds. The starting parameters of the system will change accordingly.
[32]:
plt.figure()
plt.plot(times, np.abs(hlm[(2, 2)]), label="Original system")
plt.plot(times_back, np.abs(hlm_back[(2, 2)]), "--", label="Backwards integration")
plt.xlabel("Time (seconds)")
plt.ylabel(r"$\Re[h_{22}]$")
plt.grid(True)
plt.legend()
plt.show()
Note that we can also control the secular backwards integration behavior within the GenerateWaveform class.
[33]:
m1 = 50.0
m2 = 30.0
eccentricity = 0.4
rel_anomaly = 0
approximant = "SEOBNRv5EHM"
params_dict = {
"mass1": m1,
"mass2": m2,
"approximant": approximant,
"eccentricity": eccentricity,
"rel_anomaly": rel_anomaly,
# 'secular_bwd_int' : False, # Setting this to False will deactivate the secular bwd integration
# 'warning_secular_bwd_int' : False, # Setting this to False will avoid the warning message
}
wfm_gen = GenerateWaveform(params_dict)
[34]:
# Generate mode dictionary
times, hlm = wfm_gen.generate_td_modes()
WARNING:pyseobnr.models.SEOBNRv5EHM:The predicted starting separation of the system (r_start = 7.3481088186202745) is below the minimum starting separation allowed by the model (r_start_min = 10.0). Waveform parameters are q = 1.6666666666666667, M = 80.0, chi_1 = 0.0, chi_2 = 0.0, omega_avg = 0.024758217882292533, omega_inst = 0.0539898451109042, eccentricity = 0.4, rel_anomaly = 0. Integrating backwards in time a set of secular evolution equations to obtain a prediction for the starting separation larger than this minimum. This will change the starting input values of the system.
WARNING:pyseobnr.models.SEOBNRv5EHM:The new starting input values of the system are: eccentricity = 0.40773258208912455, rel_anomaly = 4.943155348511655, omega_avg = 0.024067213761958334.
[35]:
plt.figure()
plt.plot(times, np.abs(hlm[(2, 2)]))
plt.xlabel("Time (seconds)")
plt.ylabel(r"$\Re[h_{22}]$")
plt.grid(True)
plt.show()
Using the gwsignal interface (new waveform interface)
Unlike the internal generator, the interface can accept a much wider variety of inputs both in SI and so-called ‘cosmo’ units (where say masses are in solar masses). This interface also returns the modes and polarizations in SI units and follows LAL conventions.
[36]:
import astropy.units as u
from lalsimulation.gwsignal import (
GenerateFDWaveform,
GenerateTDModes,
GenerateTDWaveform,
)
from lalsimulation.gwsignal.models import gwsignal_get_waveform_generator
[37]:
# start with the usual parameter definitions
# Masses
m1 = 30.0
m2 = 30.0
# spings
s1x = 0.0
s1y = 0.0
s1z = 0.5
s2x = 0.0
s2y = 0.0
s2z = -0.2
# Extrinsic parameters
deltaT = 1.0 / 4096.0 / 2.0
f_min = 10.0
f_ref = 10.0
luminosity_distance = 1000.0
iota = np.pi / 3.0
phase = 0.0
# Eccentric parameters
eccentricity = 0.2
longitude_ascending_nodes = 0.0
meanPerAno = 0.5
# Conditioning
condition = 1
deltaF = 1.0 / 8.0 # /20.
# Create dict for gwsignal generator
gwsignal_dict = {
"mass1": m1 * u.solMass,
"mass2": m2 * u.solMass,
"spin1x": s1x * u.dimensionless_unscaled,
"spin1y": s1y * u.dimensionless_unscaled,
"spin1z": s1z * u.dimensionless_unscaled,
"spin2x": s2x * u.dimensionless_unscaled,
"spin2y": s2y * u.dimensionless_unscaled,
"spin2z": s2z * u.dimensionless_unscaled,
#'deltaF' : delta_frequency * u.Hz,
"deltaT": deltaT * u.s,
#'deltaF': deltaF*u.Hz,
"f22_start": f_min * u.Hz,
"f_max": f_max * u.Hz,
"f22_ref": f_ref * u.Hz,
"phi_ref": phase * u.rad,
"distance": luminosity_distance * u.Mpc,
"inclination": iota * u.rad,
"eccentricity": eccentricity * u.dimensionless_unscaled,
"longAscNodes": longitude_ascending_nodes * u.rad,
"meanPerAno": meanPerAno * u.rad,
# 'ModeArray': mode_array,
"condition": 0,
"lmax": 4,
"lmax_nyquist": 4,
# Some additional specific SEOBNRv5EHM parameters
"secular_bwd_int": True,
"warning_secular_bwd_int": False,
"t_backwards": 0,
"warning_bwd_int": False,
}
waveform_approximant = "SEOBNRv5EHM"
try:
wf_gen = gwsignal_get_waveform_generator(waveform_approximant)
except ValueError as e:
if str(e) != "Approximant not implemented in GWSignal!":
raise
wf_gen = None
print("SEOBNRv5EHM not supported by this version of lal")
[38]:
if wf_gen is not None:
hpc = GenerateTDWaveform(gwsignal_dict, wf_gen)
time = np.arange(len(hpc.hp)) * gwsignal_dict["deltaT"]
plt.plot(time, hpc.hp, label=r"$h_p$")
plt.plot(time, hpc.hc, label=r"$h_x$")
plt.xlabel("Time (seconds)")
plt.ylabel(r"$h_{x,p}$")
plt.show()
else:
print("Plot skipped")
[39]:
if wf_gen is not None:
# Generate TD modes
hlm = GenerateTDModes(gwsignal_dict, wf_gen)
l, m = 2, -1
plt.plot(time, np.real(hlm[(l, m)]), label=f"({l,m})")
plt.plot(time, np.real(hlm[(l, -m)]), ls="--", label=f"({l,-m})")
plt.xlabel("Time (seconds)")
plt.ylabel(r"$\Re[h_{22}]$")
plt.legend()
plt.show()
else:
print("Plot skipped")
[40]:
deltaF = 1.0 / 9
# Create dict for gwsignal generator
gwsignal_dict = {
"mass1": m1 * u.solMass,
"mass2": m2 * u.solMass,
"spin1x": s1x * u.dimensionless_unscaled,
"spin1y": s1y * u.dimensionless_unscaled,
"spin1z": s1z * u.dimensionless_unscaled,
"spin2x": s2x * u.dimensionless_unscaled,
"spin2y": s2y * u.dimensionless_unscaled,
"spin2z": s2z * u.dimensionless_unscaled,
#'deltaF' : delta_frequency * u.Hz,
#'deltaT': deltaT*u.s,
"deltaF": deltaF * u.Hz,
"f22_start": f_min * u.Hz,
"f_max": f_max * u.Hz,
"f22_ref": f_ref * u.Hz,
"phi_ref": phase * u.rad,
"distance": luminosity_distance * u.Mpc,
"inclination": iota * u.rad,
"eccentricity": eccentricity * u.dimensionless_unscaled,
"longAscNodes": longitude_ascending_nodes * u.rad,
"meanPerAno": meanPerAno * u.rad,
# 'ModeArray': mode_array,
"condition": 1,
"lmax": 4,
"lmax_nyquist": 4,
# Some additional specific SEOBNRv5EHM parameters
"secular_bwd_int": True,
"warning_bwd_int": False,
"warning_secular_bwd_int": False,
"t_backwards": 0,
}
waveform_approximant = "SEOBNRv5EHM"
try:
wf_gen = gwsignal_get_waveform_generator(waveform_approximant)
except ValueError as e:
if str(e) != "Approximant not implemented in GWSignal!":
raise
wf_gen = None
print("SEOBNRv5EHM not supported by this version of lal")
[41]:
if wf_gen is not None:
hpc = GenerateFDWaveform(gwsignal_dict, wf_gen)
freqs = np.arange(len(hpc.hp)) * gwsignal_dict["deltaF"]
plt.plot(freqs, np.abs(hpc.hp), label=r"$h_p$")
plt.plot(freqs, np.abs(hpc.hc), label=r"$h_x$")
plt.xlabel("Frequency (Hz)")
plt.ylabel(r"$h_{p,x}$")
plt.legend()
plt.xscale("log")
plt.yscale("log")
plt.grid(True)
plt.show()
else:
print("Plot skipped")
[ ]: