lalsuite's ESIGMA branch should be installed (and sourced!) $\rightarrow$ required for the inspiral piece.NRSur7dq4's data file should be downloaded, with the path to the directory containing it been set to the bash environment variable LAL_DATA_PATH $\rightarrow$ required for the merger-ringdown piece.import numpy as np
from gwnr.waveform import esigma_utils as eu
import matplotlib.pyplot as plt
The polarizations $h_+$ and $h_\times$ can be generated via the get_imr_esigma_waveform function. They are returned as PyCBC TimeSeries objects.
from gwnr.waveform import esigma_utils as eu
m1 = 20. # masses (in solar masses)
m2 = 30.
spin1z = 0.5 # dimensionless spins
spin2z = -0.3
eccentricity = 0.15 # starting eccentricity
mean_anomaly = 60 * np.pi/180. # starting mean anomaly
modes_to_use = [(2,2), (2,-2)] # Only using the (2,2) mode
distance = 400. # source luminosity distance (in Mpc)
inclination = 30 * np.pi/180. # orbital inclination with line-of-sight
f_low = 20. # starting frequency (in Hz)
delta_t = 1/2**12 # time grid-spacing (in s)
hp, hc = eu.get_imr_esigma_waveform(mass1=m1,
mass2=m2,
spin1z=spin1z,
spin2z=spin2z,
eccentricity=eccentricity,
mean_anomaly=mean_anomaly,
distance=distance,
inclination=inclination,
f_lower=f_low,
delta_t=delta_t,
modes_to_use=modes_to_use)
plt.figure(figsize=(10,4))
plt.title(fr"""$m_1={m1} M_\odot, m_2={m2} M_\odot, \chi_{{1z}}={spin1z}, \chi_{{2z}}={spin2z}$
$e_0={eccentricity}, l_0={mean_anomaly:.2f}, f_{{\rm{{low}}}}={f_low}\rm{{Hz}}$
$D_L={distance}\rm{{Mpc}}, \iota={inclination:.2f}$""")
hp.plot(label=r"$h_+$")
hc.plot(label=r"$h_\times$")
plt.xlabel(r"$t (s)$")
plt.ylabel(r"$h$")
plt.legend()
<matplotlib.legend.Legend at 0x7fd576ecca00>
One can also specify the eccentricity and mean anomaly at some reference frequency $f_\rm{ref}$ different from the starting frequency $f_{\rm{low}}$ of the waveform. However, currently we only support $f_{\rm{ref}} \leq f_{\rm{low}}$.
f_ref = 12.
f_low = 20.
# This time, the following eccentricity and mean anomaly values are
# defined at f_ref, and not at f_low
eccentricity = 0.15
mean_anomaly = 30. * np.pi/180.
hp, hc = eu.get_imr_esigma_waveform(mass1=m1,
mass2=m2,
spin1z=spin1z,
spin2z=spin2z,
eccentricity=eccentricity,
mean_anomaly=mean_anomaly,
distance=distance,
inclination=inclination,
f_ref=f_ref,
f_lower=f_low,
delta_t=delta_t,
modes_to_use=modes_to_use)
plt.figure(figsize=(10,4))
plt.title(fr"""$m_1={m1} M_\odot, m_2={m2} M_\odot, \chi_{{1z}}={spin1z}, \chi_{{2z}}={spin2z}$
$e_{{\rm{{ref}}}}={eccentricity}, l_{{\rm{{ref}}}}={mean_anomaly:.2f}, f_{{\rm{{ref}}}}={f_ref}\rm{{Hz}}, f_{{\rm{{low}}}}={f_low}\rm{{Hz}}$
$D_L={distance}\rm{{Mpc}}, \iota={inclination:.2f}$""")
hp.plot(label=r"$h_+$")
hc.plot(label=r"$h_\times$")
plt.xlabel(r"$t (s)$")
plt.ylabel(r"$h$")
plt.legend()
<matplotlib.legend.Legend at 0x7fd575cfb580>
One can also generate the spin-weighted spherical harmonic modes via the get_imr_esigma_modes function.
m1 = 10. # masses (in solar masses)
m2 = 40.
spin1z = 0.7 # dimensionless spins
spin2z = 0.6
eccentricity = 0.15 # starting eccentricity
mean_anomaly = 60 * np.pi/180. # starting mean anomaly
distance = 400. # source luminosity distance (in Mpc)
inclination = 30 * np.pi/180. # orbital inclination with line-of-sight
f_low = 20. # starting frequency (in Hz)
delta_t = 1/2**12 # time grid-spacing (in s)
modes_to_use=[(2,2), (2,1), (3,3), (3,2), (4,4), (4,3)]
modes = eu.get_imr_esigma_modes(mass1=m1,
mass2=m2,
spin1z=spin1z,
spin2z=spin2z,
eccentricity=eccentricity,
mean_anomaly=mean_anomaly,
distance=distance,
f_lower=f_low,
delta_t=delta_t,
modes_to_use=modes_to_use,
include_conjugate_modes=False)
fig, axs = plt.subplots(len(modes_to_use), sharex=True, figsize=(10, 10))
axs[0].set_title(fr"""$m_1={m1} M_\odot, m_2={m2} M_\odot, \chi_{{1z}}={spin1z}, \chi_{{2z}}={spin2z}$
$e_0={eccentricity}, l_0={mean_anomaly:.2f}, f_{{\rm{{low}}}}={f_low}\rm{{Hz}}, D_L={distance}\rm{{Mpc}}$""")
for i, mode_name in enumerate(modes_to_use):
ell, m = mode_name
axs[i].plot(modes[mode_name].sample_times.data, modes[mode_name].real().data, label=rf"$\Re(h_{{{ell} {m}}})$")
axs[i].plot(modes[mode_name].sample_times.data, modes[mode_name].imag().data, label=rf"$\Im(h_{{{ell} {m}}})$")
axs[i].legend(loc=2)
plt.xlabel(r"$t (s)$")
Text(0.5, 0, '$t (s)$')
Other options related to hybridization settings, merger-ringdown model choice, etc. are available as well. One can check them in the respective docstrings of the functions like so:
help(eu.get_imr_esigma_waveform)
Help on function get_imr_esigma_waveform in module gwnr.waveform.esigma_utils:
get_imr_esigma_waveform(mass1, mass2, f_lower, delta_t, f_ref=None, spin1z=0.0, spin2z=0.0, eccentricity=0.0, mean_anomaly=0.0, coa_phase=0.0, inclination=0.0, distance=1.0, modes_to_use=[(2, 2), (3, 3), (4, 4)], mode_to_align_by=(2, 2), f_mr_transition=None, f_window_mr_transition=None, num_hyb_orbits=0.25, hybridize_using_avg_orbital_frequency=True, hybridize_aligning_merger_to_inspiral=True, keep_f_mr_transition_at_center=False, merger_ringdown_approximant='NRSur7dq4', return_hybridization_info=False, return_orbital_params=False, failsafe=True, verbose=False, **kwargs)
Returns IMR GW polarizations constructed using IMR ESIGMA modes
Parameters:
-----------
mass1, mass2 -- Binary's component masses (in solar masses)
f_lower -- Starting frequency of the waveform (in Hz)
f_ref -- Reference frequency at which to define the
waveform parameters. We require that
`f_ref <= f_lower`.
`f_ref = f_lower` by default.
delta_t -- Waveform's time grid-spacing (in s)
spin1z, spin2z -- z-components of component dimensionless
spins (lies in [0,1))
eccentricity -- Initial eccentricity
mean_anomaly -- Mean anomaly of the periastron (in rad)
inclination -- Inclination (in rad), defined as the angle
between the orbital angular momentum L and
the line-of-sight
coa_phase -- Coalesence phase of the binary (in rad)
distance -- Luminosity distance to the binary (in Mpc)
modes_to_use -- GW modes to use. List of tuples (l, |m|)
mode_to_align_by -- GW mode to use to align inspiral and merger
in phase and time
f_mr_transition -- Inspiral to merger transition GW frequency
(Hz).
Defaults to the minimum of the Kerr and
Schwarzschild ISCO frequency
f_window_mr_transition -- Hybridization frequency window (in Hz).
Disabled by the default value (None). In
such a case, the hybridization proceeds
over a window of `num_hyb_orbits` orbital
cycles (1 orbital cycle ~ 2 GW cycles)
that ends at the frequency value given by
`f_mr_transition`.
Also see `keep_f_mr_transition_at_center`
to choose the position of `f_mr_transition`
within this window.
num_hyb_orbits -- number of orbital cycles to hybridize over.
Only used if f_window_mr_transition is not
specified.
hybridize_using_avg_orbital_frequency -- If True, the orbit averaged
frequency during the inspiral is used to
hybridize modes, instead of the modes'
frequency.
keep_f_mr_transition_at_center -- If True, `f_mr_transition` is kept at
the center of the hybridization window.
Otherwise, it's kept at the end of the
window (default).
merger_ringdown_approximant -- Choose merger-ringdown model. Tested
choices: [NRSur7dq4, SEOBNRv4PHM]
return_hybridization_info -- If True, returns hybridization related data
return_orbital_params -- If True, returns the orbital evolution of
all the orbital elements (in
geometrized units). Can also be a list of
orbital variable names to return
only those specific variables. Available
orbital variables names are:
['x', 'e', 'l', 'phi', 'phidot', 'r', 'rdot'].
Note that these are available only for the
inspiral portion of the waveform!
failsafe -- If True, we make reasonable choices for the
user, if the inputs to this method lead
into exceptions.
verbose -- Verbosity level. Available values are: 0, 1, 2
Returns:
--------
hp, hc -- Plus and cross IMR GW polarizations PyCBC TimeSeries
orbital_vars_dict -- Dictionary of evolution of orbital elements.
Returned only if return_orbital_params is specified
retval -- Hybridization related data.
Returned only if return_hybridization_info is True
One can also generate the inspiral-only wavefroms and modes via the get_inspiral_esigma_waveform and get_inspiral_esigma_modes functions respectively. Unlike the full IMR waveform which uses a quasi-circular merger-ringdown model and thus forbids using large starting eccentricities, the inspiral waveforms can be generated for arbitrarily eccentric (bounded) orbits (i.e. with $0 \leq e < 1$).
m1 = 10. # masses (in solar masses)
m2 = 50.
spin1z = 0.5 # dimensionless spins
spin2z = 0.6
eccentricity = 0.55 # starting eccentricity
mean_anomaly = 60 * np.pi/180. # starting mean anomaly
modes_to_use = [(2,2), (2,-2)] # Only using the (2,2) mode
distance = 400. # source luminosity distance (in Mpc)
inclination = 30 * np.pi/180. # orbital inclination with line-of-sight
f_low = 8 # starting frequency (in Hz)
delta_t = 1/2**12 # time grid-spacing (in s)
hp, hc = eu.get_inspiral_esigma_waveform(mass1=m1,
mass2=m2,
spin1z=spin1z,
spin2z=spin2z,
eccentricity=eccentricity,
mean_anomaly=mean_anomaly,
distance=distance,
inclination=inclination,
f_lower=f_low,
delta_t=delta_t,
modes_to_use=modes_to_use
)
plt.figure(figsize=(10,4))
plt.title(fr"""$m_1={m1} M_\odot, m_2={m2} M_\odot, \chi_{{1z}}={spin1z}, \chi_{{2z}}={spin2z}$
$e_0={eccentricity}, l_0={mean_anomaly:.2f}, f_{{\rm{{low}}}}={f_low}\rm{{Hz}}$
$D_L={distance}\rm{{Mpc}}, \iota={inclination:.2f}$""")
hp.plot(label=r"$h_+$")
hc.plot(label=r"$h_\times$")
plt.xlabel(r"$t (s)$")
plt.ylabel(r"$h$")
plt.legend()
<matplotlib.legend.Legend at 0x7fd5759eb190>
modes_to_use=[(2,2), (2,1), (3,3), (3,2), (4,4), (4,3)]
modes = eu.get_inspiral_esigma_modes(mass1=m1,
mass2=m2,
spin1z=spin1z,
spin2z=spin2z,
eccentricity=eccentricity,
mean_anomaly=mean_anomaly,
distance=distance,
f_lower=f_low,
delta_t=delta_t,
modes_to_use=modes_to_use,
include_conjugate_modes=False)
fig, axs = plt.subplots(len(modes_to_use), sharex=True, figsize=(10, 10))
axs[0].set_title(fr"""$m_1={m1} M_\odot, m_2={m2} M_\odot, \chi_{{1z}}={spin1z}, \chi_{{2z}}={spin2z}$
$e_0={eccentricity}, l_0={mean_anomaly:.2f}, f_{{\rm{{low}}}}={f_low}\rm{{Hz}}, D_L={distance}\rm{{Mpc}}$""")
for i, mode_name in enumerate(modes_to_use):
ell, m = mode_name
axs[i].plot(modes[mode_name].sample_times.data, modes[mode_name].real().data, label=rf"$\Re(h_{{{ell} {m}}})$")
axs[i].plot(modes[mode_name].sample_times.data, modes[mode_name].imag().data, label=rf"$\Im(h_{{{ell} {m}}})$")
axs[i].legend(loc=2)
plt.xlabel(r"$t (s)$")
Text(0.5, 0, '$t (s)$')
The evolution of binary's orbital elements can also be accessed, but only for the inspiral part of the dynamics. These can be accessed via the argument return_orbital_params in all of the above discussed waveform/mode functions. The available orbital elements are
All of these are returned in geometric units ($G=c=1$).
m1 = 10. # masses (in solar masses)
m2 = 50.
spin1z = 0.5 # dimensionless spins
spin2z = 0.6
eccentricity = 0.15 # starting eccentricity
mean_anomaly = 60 * np.pi/180. # starting mean anomaly
modes_to_use = [(2,2), (2,-2)] # Only using the (2,2) mode
orb_params_list = ['x', 'e', 'l', 'r', 'rdot', 'phi', 'phidot' ]
distance = 400. # source luminosity distance (in Mpc)
inclination = 30 * np.pi/180. # orbital inclination with line-of-sight
f_low = 20 # starting frequency (in Hz)
delta_t = 1/2**12 # time grid-spacing (in s)
orb_vars, hp, hc = eu.get_inspiral_esigma_waveform(mass1=m1,
mass2=m2,
spin1z=spin1z,
spin2z=spin2z,
eccentricity=eccentricity,
mean_anomaly=mean_anomaly,
distance=distance,
inclination=inclination,
f_lower=f_low,
delta_t=delta_t,
modes_to_use=modes_to_use,
return_orbital_params=orb_params_list)
fig, axs = plt.subplots(len(orb_params_list), sharex=True, figsize=(10, 10))
axs[0].set_title(fr"$m_1={m1} M_\odot, m_2={m2} M_\odot, \chi_{{1z}}={spin1z}, \chi_{{2z}}={spin2z}, e_0={eccentricity}, l_0={mean_anomaly:.2f}, f_{{\rm{{low}}}}={f_low}\rm{{Hz}}$")
for i, orb_params_name in enumerate(orb_params_list):
axs[i].plot(orb_vars[orb_params_name].sample_times.data, orb_vars[orb_params_name].data, label=rf"{orb_params_name}")
axs[i].legend(loc=2)
plt.xlabel(r"$t (s)$")
Text(0.5, 0, '$t (s)$')