PySEOBNR introduction

The notebook shows how to get started with pyseobnr.

import warnings

# silence warnings coming from LAL
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")

import numpy as np
from matplotlib import pyplot as plt

# import the library
from pyseobnr.generate_waveform import GenerateWaveform

# Start with the usual parameter definitions
# Masses in solar masses
m1 = 50.0
m2 = 30.0
s1x, s1y, s1z = 0.0, 0.0, 0.5
s2x, s2y, s2z = 0.0, 0.0, 0.8

deltaT = 1.0 / 2048.0
f_min = 20.0
f_max = 1024.0

distance = 1000.0  # Mpc
inclination = np.pi / 3.0
phiRef = 0.0
approximant = "SEOBNRv5HM"

params_dict = {
    "mass1": m1,
    "mass2": m2,
    "spin1x": s1x,
    "spin1y": s1y,
    "spin1z": s1z,
    "spin2x": s2x,
    "spin2y": s2y,
    "spin2z": s2z,
    "deltaT": deltaT,
    "f22_start": f_min,
    "phi_ref": phiRef,
    "distance": distance,
    "inclination": inclination,
    "f_max": f_max,
    "approximant": approximant,
}

# We call the generator with the parameters
wfm_gen = GenerateWaveform(params_dict)

# Generate mode dictionary
times, hlm = wfm_gen.generate_td_modes()

Plot some 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()

Expert mode

from pyseobnr.generate_waveform import generate_modes_opt

q = 5.3
chi_1 = 0.9
chi_2 = 0.3
omega0 = 0.0137  # This is the orbital frequency in geometric units with M=1

_, _, model = generate_modes_opt(q, chi_1, chi_2, omega0, debug=True)
model

<pyseobnr.models.SEOBNRv5HM.SEOBNRv5HM_opt at 0x7fd58c57fbb0>

t, r, phi, pr, pphi, H, Omega, _ = model.dynamics.T

import pandas as pd

frame_dynamics = pd.DataFrame(
    data=model.dynamics,
    columns="t, r, phi, pr, pphi, H, Omega, Omega_circular".replace(" ", "").split(","),
)
frame_dynamics

	t	r	phi	pr	pphi	H	Omega	Omega_circular
0	0.000000	17.371216	0.000000	-0.000325	4.436354	7.460756	0.013700	0.013700
1	91.725333	17.343818	1.258109	-0.000325	4.433223	7.460713	0.013732	0.013732
2	183.450666	17.316357	2.519172	-0.000327	4.430076	7.460670	0.013764	0.013764
3	275.175999	17.288693	3.783213	-0.000330	4.426911	7.460626	0.013797	0.013797
4	366.901332	17.260762	5.050276	-0.000333	4.423729	7.460582	0.013830	0.013830
...	...	...	...	...	...	...	...	...
2805	14876.512636	2.088023	326.953577	-0.231948	2.171192	7.334665	0.207106	0.214091
2806	14876.612636	2.078014	326.974325	-0.235755	2.169841	7.334335	0.207820	0.215069
2807	14876.712636	2.067742	326.995145	-0.239688	2.168453	7.333993	0.208560	0.216090
2808	14876.812636	2.057193	327.016039	-0.243751	2.167025	7.333637	0.209326	0.217158
2809	14876.912636	2.046353	327.037009	-0.247950	2.165557	7.333268	0.210121	0.218275

2810 rows × 8 columns

Selecting different convention for reference phase and \(t=0\)

By default, pyseobnr will generate the waveforms with the following conventions for phase and time:

• at the reference frequency, the orbital phase will be set to 0,

• the origin of the time axis, \(t=0\), is set to coincide with the peak of the frame invariant amplitude \(A\):

  \[A(t) = \sqrt{\sum_{\ell m}|h_{\ell m}(t)|}\newline\]
  \[{\rm argmax}_t A = 0\]

Other models, in particular the IMRPhenomT models, follow a different convention for setting the reference phase and the time origin:

• at the reference frequency, the phase of the co-precessing dominant \((2,2)\) mode is set to 0, while the other modes are rotated consistenly to preserve the correct relative phase with the dominant mode,

• the time origin is set to coincide with the peak amplitude of the co-precessing \((2,2)\) mode.

One can set this second convention in pyseobnr via the following options:

• convention_coprocessing_phase22_set_to_0_at_reference_time to True

• convention_t0_set_to_0_at_coprocessing_amplitude22_peak to True

Note that when convention_coprocessing_phase22_set_to_0_at_reference_time is True as well as (for precessing systems) polarizations_from_coprec is True, the phase rotations will be applied directly to the spin weighted spherical harmonics.

# Start with the usual parameter definitions
# Masses in solar masses
m1 = 50.0
m2 = 10.0
s1x, s1y, s1z = 0.0, 0.0, 0.5
s2x, s2y, s2z = 0.0, 0.0, 0.8

deltaT = 1.0 / 2048.0
f_min = 20.0
f_max = 1024.0

distance = 1000.0  # Mpc
inclination = np.pi / 3.0
phiRef = 0.0
approximant = "SEOBNRv5HM"

params_dict = {
    "mass1": m1,
    "mass2": m2,
    "spin1x": s1x,
    "spin1y": s1y,
    "spin1z": s1z,
    "spin2x": s2x,
    "spin2y": s2y,
    "spin2z": s2z,
    "deltaT": deltaT,
    "f22_start": f_min,
    "phi_ref": phiRef,
    "distance": distance,
    "inclination": inclination,
    "f_max": f_max,
    "approximant": approximant,
}

# We call the generator with the parameters
wfm_gen = GenerateWaveform(params_dict)

# We update with the new convention
wfm_gen_2 = GenerateWaveform(
    params_dict
    | {
        "convention_coprocessing_phase22_set_to_0_at_reference_time": True,
        "convention_t0_set_to_0_at_coprocessing_amplitude22_peak": True,
    }
)

# Generate mode dictionary
times_0, hlms_0 = wfm_gen.generate_td_modes()
times_1, hlms_1 = wfm_gen_2.generate_td_modes()


plt.figure()
plt.plot(times_0, hlms_0[(2, 2)].real, label="default pyseobnr convention")
plt.plot(times_1, hlms_1[(2, 2)].real, label="IMRPhenomT convention", ls="dashed")
plt.xlabel("Time (seconds)")
plt.ylabel(r"$\Re[h_{22}]$")
plt.legend()
plt.xlim([-0.05, 0.01])
plt.grid(True)
plt.tight_layout()