Example notebook to use BHPTNRSur1dq1e4-gwNRHME model

BHPTNRSur1dq1e4-gwNRHME (multi-modal eccentric nonspinning model) = SEOBNRv5EHM (22 mode eccentric model) + BHPTNRSur1dq1e4 (multi-modal circular nonspinning model) + gwNRHME framework

[1]:

import sys
# install gwModels in editable mode from the parent directory
!{sys.executable} -m pip install -e ../ --no-deps --quiet

import matplotlib.pyplot as plt
import numpy as np
import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")

import gwModels
gwModels.utils.set_rcparams()

lal.MSUN_SI != Msun

[2]:

import gwsurrogate
import os, contextlib

# setup model for BHPTNRSur1dq1e4 — update path_to_surrogate for your environment
path_to_gws = gwsurrogate.__path__[0] + '/'
path_to_surrogate = path_to_gws + 'surrogate_downloads/BHPTNRSur1dq1e4.h5'

with contextlib.redirect_stdout(open(os.devnull, 'w')):
    model_obj = gwsurrogate.EvaluateSurrogate(path_to_surrogate,
                    ell_m=[(2,2),(2,1),(3,1),(3,2),(3,3),(4,2),(4,3),(4,4)])

1. Generate waveforms by combining BHPTNRSur1dq1e4 (multi-mode) and SEOBNRv5EHM (22 mode)

[3]:

# Set the binary parameters
params = {"q": 2.6,       # mass ratio
          "omega0": 0.02, # initial dimensionless orbital frequency
          "e0": 0.1,      # initial eccentricity
          "l0": 0}        # initial relativistic anomaly

[4]:

# generate eccentric (2,2) mode using SEOBNRv5EHM
seob = gwModels.waveforms.genSEOBNRv5EHM()
tecc, hecc_dict = seob.generate_SEOBNRv5EHM(params)
hecc = hecc_dict['h_l2m2']

# generate circular higher modes using BHPTNRSur1dq1e4
bhpt = gwModels.waveforms.genBHPTNRSur1dq1e4(model_obj)
tsur, hsur = bhpt.generate_BHPTNRSur1dq1e4(model_obj, params)

# combine using gwNRHME framework
gwnrhme = gwModels.frameworks.NRHME(t_ecc=tecc,
                                     h_ecc_dict={'h_l2m2': hecc},
                                     t_cir=tsur,
                                     h_cir_dict=hsur)

tNRE = gwnrhme.t_common
hNRE = gwnrhme.hNRE

[5]:

# print available modes
print(hNRE.keys())

dict_keys(['h_l2m2', 'h_l2m1', 'h_l3m3', 'h_l3m2', 'h_l4m4', 'h_l4m3'])

[6]:

# plot waveform
plt.figure(figsize=(10,9))
plt.subplot(411)
plt.plot(tNRE, hNRE['h_l2m2'].real, '-', lw=2, color='C0', label='(2,2)')
plt.xlabel('$t/M$')
plt.ylabel('$rh_{22}/M$')
plt.subplot(412)
plt.plot(tNRE, hNRE['h_l2m1'].real, '-', lw=2, color='C1', label='(2,1)')
plt.xlabel('$t/M$')
plt.ylabel('$rh_{21}/M$')
plt.subplot(413)
plt.plot(tNRE, hNRE['h_l3m3'].real, '-', lw=2, color='C2', label='(3,3)')
plt.xlabel('$t/M$')
plt.ylabel('$rh_{33}/M$')
plt.subplot(414)
plt.plot(tNRE, hNRE['h_l4m4'].real, '-', lw=2, color='C3', label='(4,4)')
plt.xlabel('$t/M$')
plt.ylabel('$rh_{44}/M$')
plt.tight_layout()
plt.show()

../_images/notebooks_2_2_BHPTNRSur1dq1e4-gwNRHME_example_7_0.png

[7]:

# plot amplitudes
plt.figure(figsize=(10,9))
plt.subplot(411)
plt.plot(tNRE, abs(hNRE['h_l2m2']), '-', lw=2, color='C0', label='(2,2)')
plt.xlabel('$t/M$')
plt.ylabel('$|rh_{22}/M|$')
plt.subplot(412)
plt.plot(tNRE, abs(hNRE['h_l2m1']), '-', lw=2, color='C1', label='(2,1)')
plt.xlabel('$t/M$')
plt.ylabel('$|rh_{21}/M|$')
plt.subplot(413)
plt.plot(tNRE, abs(hNRE['h_l3m3']), '-', lw=2, color='C2', label='(3,3)')
plt.xlabel('$t/M$')
plt.ylabel('$|rh_{33}/M|$')
plt.subplot(414)
plt.plot(tNRE, abs(hNRE['h_l4m4']), '-', lw=2, color='C3', label='(4,4)')
plt.xlabel('$t/M$')
plt.ylabel('$|rh_{44}/M|$')
plt.tight_layout()
plt.show()

../_images/notebooks_2_2_BHPTNRSur1dq1e4-gwNRHME_example_8_0.png

[8]:

# plot frequencies
from gwModels.utils.features import get_frequency

plt.figure(figsize=(10,9))
plt.subplot(411)
plt.plot(tNRE, get_frequency(tNRE, hNRE['h_l2m2']), '-', lw=2, color='C0', label='(2,2)')
plt.xlabel('$t/M$')
plt.ylabel('$M\\omega_{22}$')
plt.subplot(412)
plt.plot(tNRE, get_frequency(tNRE, hNRE['h_l2m1']), '-', lw=2, color='C1', label='(2,1)')
plt.xlabel('$t/M$')
plt.ylabel('$M\\omega_{21}$')
plt.subplot(413)
plt.plot(tNRE, get_frequency(tNRE, hNRE['h_l3m3']), '-', lw=2, color='C2', label='(3,3)')
plt.xlabel('$t/M$')
plt.ylabel('$M\\omega_{33}$')
plt.subplot(414)
plt.plot(tNRE, get_frequency(tNRE, hNRE['h_l4m4']), '-', lw=2, color='C3', label='(4,4)')
plt.xlabel('$t/M$')
plt.ylabel('$M\\omega_{44}$')
plt.tight_layout()
plt.show()

../_images/notebooks_2_2_BHPTNRSur1dq1e4-gwNRHME_example_9_0.png

2. Generate waveforms from combined model using BHPTNRSur1dq1e4_gwNRHME

This is a convenience wrapper that does all three steps (SEOBNRv5EHM + BHPTNRSur1dq1e4 + gwNRHME) in a single call.

[9]:

# instantiate and generate waveform in one go
wf = gwModels.waveforms.BHPTNRSur1dq1e4_gwNRHME(model_obj=model_obj)
tNRE2, hNRE2 = wf.generate_waveform(params)
print(hNRE2.keys())

dict_keys(['h_l2m2', 'h_l2m1', 'h_l3m3', 'h_l3m2', 'h_l4m4', 'h_l4m3'])

[10]:

# plot waveform
plt.figure(figsize=(10,9))
plt.subplot(411)
plt.plot(tNRE2, hNRE2['h_l2m2'].real, '-', lw=2, color='C0', label='(2,2)')
plt.xlabel('$t/M$')
plt.ylabel('$rh_{22}/M$')
plt.subplot(412)
plt.plot(tNRE2, hNRE2['h_l2m1'].real, '-', lw=2, color='C1', label='(2,1)')
plt.xlabel('$t/M$')
plt.ylabel('$rh_{21}/M$')
plt.subplot(413)
plt.plot(tNRE2, hNRE2['h_l3m3'].real, '-', lw=2, color='C2', label='(3,3)')
plt.xlabel('$t/M$')
plt.ylabel('$rh_{33}/M$')
plt.subplot(414)
plt.plot(tNRE2, hNRE2['h_l4m4'].real, '-', lw=2, color='C3', label='(4,4)')
plt.xlabel('$t/M$')
plt.ylabel('$rh_{44}/M$')
plt.tight_layout()
plt.show()

../_images/notebooks_2_2_BHPTNRSur1dq1e4-gwNRHME_example_12_0.png

[11]:

# plot amplitudes
plt.figure(figsize=(10,9))
plt.subplot(411)
plt.plot(tNRE2, abs(hNRE2['h_l2m2']), '-', lw=2, color='C0', label='(2,2)')
plt.xlabel('$t/M$')
plt.ylabel('$|rh_{22}/M|$')
plt.subplot(412)
plt.plot(tNRE2, abs(hNRE2['h_l2m1']), '-', lw=2, color='C1', label='(2,1)')
plt.xlabel('$t/M$')
plt.ylabel('$|rh_{21}/M|$')
plt.subplot(413)
plt.plot(tNRE2, abs(hNRE2['h_l3m3']), '-', lw=2, color='C2', label='(3,3)')
plt.xlabel('$t/M$')
plt.ylabel('$|rh_{33}/M|$')
plt.subplot(414)
plt.plot(tNRE2, abs(hNRE2['h_l4m4']), '-', lw=2, color='C3', label='(4,4)')
plt.xlabel('$t/M$')
plt.ylabel('$|rh_{44}/M|$')
plt.tight_layout()
plt.show()

../_images/notebooks_2_2_BHPTNRSur1dq1e4-gwNRHME_example_13_0.png

[12]:

# plot frequencies
plt.figure(figsize=(10,9))
plt.subplot(411)
plt.plot(tNRE2, get_frequency(tNRE2, hNRE2['h_l2m2']), '-', lw=2, color='C0', label='(2,2)')
plt.xlabel('$t/M$')
plt.ylabel('$M\\omega_{22}$')
plt.subplot(412)
plt.plot(tNRE2, get_frequency(tNRE2, hNRE2['h_l2m1']), '-', lw=2, color='C1', label='(2,1)')
plt.xlabel('$t/M$')
plt.ylabel('$M\\omega_{21}$')
plt.subplot(413)
plt.plot(tNRE2, get_frequency(tNRE2, hNRE2['h_l3m3']), '-', lw=2, color='C2', label='(3,3)')
plt.xlabel('$t/M$')
plt.ylabel('$M\\omega_{33}$')
plt.subplot(414)
plt.plot(tNRE2, get_frequency(tNRE2, hNRE2['h_l4m4']), '-', lw=2, color='C3', label='(4,4)')
plt.xlabel('$t/M$')
plt.ylabel('$M\\omega_{44}$')
plt.tight_layout()
plt.show()

../_images/notebooks_2_2_BHPTNRSur1dq1e4-gwNRHME_example_14_0.png

[ ]: