Example usage of SEOBNRv4PHM surrogate model¶

Author: Bhooshan Gadre

In [1]:

import numpy as np
import matplotlib.pyplot as P
%matplotlib inline

import gwsurrogate

import numpy as np import matplotlib.pyplot as P %matplotlib inline import gwsurrogate

lal.MSUN_SI != Msun
__name__ = gwsurrogate.new.spline_evaluation
__package__= gwsurrogate.new

/home/balzani57/pool/anaconda3/envs/py311-gw-5-2024-dev/lib/python3.11/site-packages/gwtools/const.py:52: 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

Download surrogate data, this only needs to be done once¶

In [2]:

# This can take a few minutes
# gwsurrogate.catalog.pull('SEOBNRv4PHMSur')

# This can take a few minutes # gwsurrogate.catalog.pull('SEOBNRv4PHMSur')

Out[2]:

'/Users/vijay/src/gwsurrogate/gwsurrogate/surrogate_downloadsNRSur7dq4.h5'

Load the surrogate, this only needs to be done once at the start of a script¶

In [2]:

sur = gwsurrogate.LoadSurrogate('SEOBNRv4PHMSur')

sur = gwsurrogate.LoadSurrogate('SEOBNRv4PHMSur')

Loaded SEOBNRv4PHMSur model

In [3]:

sur.__dict__['_sur_dimless'].mode_list

sur.__dict__['_sur_dimless'].mode_list

Out[3]:

[(2, 1),
 (2, -1),
 (2, 2),
 (2, -2),
 (3, 3),
 (3, -3),
 (4, 4),
 (4, -4),
 (5, 5),
 (5, -5)]

Read the documentation¶

In [3]:

 help(sur)

help(sur)

Help on SEOBNRv4PHMSur in module gwsurrogate.surrogate object:

class SEOBNRv4PHMSur(SurrogateEvaluator)
 |  SEOBNRv4PHMSur(h5filename)
 |  
 |  A class for the SEOBNRv4PHM surrogate models.
 |  
 |  Evaluates gravitational waveforms generated by precessing binary black hole
 |  systems with generic mass ratios and spins.
 |  
 |  This model includes the following spin-weighted spherical harmonic modes:
 |  2<=ell<=5, -ell<=m<=ell in inertial frame.
 |  
 |  The parameter space of validity is:
 |  q \in [1, 20], and |chi1|,|chi2| \in [-0.99, 0.99], with generic directions.
 |  where q is the mass ratio and chi1/chi2 are the spin vectors of the
 |  heavier/lighter BH, respectively.
 |  
 |  The surrogate has been trained in the range
 |  q \in [1, 20] and |chi1|/|chi2| \in [-0.8, 0.8], but produces reasonable
 |  waveforms in the above range and has been tested against existing
 |  NR waveforms in that range.
 |  
 |  See the __call__ method on how to evaluate waveforms.
 |  In the __call__ method, x must have format x = [q, chi1, chi2].
 |  
 |  
 |  IMPORTANT NOTES:
 |  ===============
 |  
 |  The original SEOBNRv4PHM model (arXiv:2004.09442) parameterizes the 
 |  direction of the BH spins relative to the Newtonian orbital angular momentum
 |  vector of the binary in the co-precessing frame. The SEOBNRv4PHMSur model, 
 |  however, parameterizes the spins relative to the direction where the radiation
 |  is always strongest along the z-axis, and the (ell,m)= (2, ±2) modes are dominant.
 |  This is the same convention used in the NRSur7dq4 surrogate model.
 |  
 |  These two frames are different for non-precessing systems. Hence, specifying
 |  the spins in the SEOBNRv4PHM-frame (e.g. when called through LALSimulation)
 |  will not give a physically equivalent waveform when these same spin directions
 |  are used in the SEOBNRv4PHMSur surrogate model. Ref. arXiv:2203.00381 provides
 |  further discussion on this point.
 |  
 |  Method resolution order:
 |      SEOBNRv4PHMSur
 |      SurrogateEvaluator
 |      builtins.object
 |  
 |  Methods defined here:
 |  
 |  __init__(self, h5filename)
 |      name:           Name of the surrogate
 |      domain_type:    'Time' or 'Frequency'
 |      keywords:       keywords for this model. For allowed keys see
 |                      self._check_keywords_and_set_defaults.default_keywords.
 |                      If keywords['Precessing'] = False, will automatically
 |                      determine the m<0 modes from the m>0 modes.
 |      soft_param_lims: Parameter bounds beyond which a warning is raised.
 |      hard_param_lims: Parameter bounds beyond which an error is raised.
 |                       Should be in format [qMax, chimax]
 |                       Setting soft_param_lims/hard_param_lims to None will
 |                       skip that particular check.
 |  
 |  ----------------------------------------------------------------------
 |  Methods inherited from SurrogateEvaluator:
 |  
 |  __call__(self, q, chiA0, chiB0, M=None, dist_mpc=None, f_low=None, f_ref=None, dt=None, df=None, times=None, freqs=None, mode_list=None, ellMax=None, inclination=None, phi_ref=0, precessing_opts=None, tidal_opts=None, par_dict=None, units='dimensionless', skip_param_checks=False, taper_end_duration=None)
 |      INPUT
 |      =====
 |      q :         Mass ratio, mA/mB >= 1.
 |      chiA0:      Dimensionless spin vector of the heavier black hole at
 |                  reference epoch.
 |      chiB0:      Dimensionless spin vector of the lighter black hole at
 |                  reference epoch.
 |      
 |                  This follows the same convention as LAL, where the spin
 |                  components are defined as:
 |                  \chi_z = \chi \cdot \hat{L}, where L is the orbital angular
 |                      momentum vector at the epoch.
 |                  \chi_x = \chi \cdot \hat{n}, where n = body2 -> body1 is the
 |                      separation vector at the epoch. body1 is the heavier body.
 |                  \chi_y = \chi \cdot \hat{L \cross n}.
 |                  These spin components are frame-independent as they are
 |                  defined using vector inner products. This is equivalent to
 |                  specifying the spins in the coorbital frame used in the
 |                  surrogate papers.
 |      
 |      M, dist_mpc: Either specify both M and dist_mpc or neither.
 |          M        :  Total mass (solar masses). Default: None.
 |          dist_mpc :  Distance to binary system (MegaParsecs). Default: None.
 |      
 |      f_low :     Instantaneous initial frequency of the (2, 2) mode. In
 |                  practice, this is estimated to be twice the initial orbital
 |                  frequency in the coprecessing frame. Note: the coprecessing
 |                  frame is the minimal rotation frame of arXiv:1110.2965.
 |      
 |                  f_low should be in cycles/M if units = 'dimensionless',
 |                  should be in Hertz if units = 'mks'.
 |                  If 0, the entire waveform is returned.
 |                  Default: None, must be specified by user.
 |      
 |                  NOTE: For some models like NRSur7dq4, f_low=0 is recommended.
 |                  The role of f_low is only to truncate the lower frequencies
 |                  before returning the waveform. Since this model is already
 |                  very short, this truncation is not required. On the other hand,
 |                  f_ref is used to set the reference epoch, and can be freely
 |                  specified.
 |      
 |                  WARNING: Using f_low=0 with a small dt (like 0.1M) can lead to
 |                  very expensive evaluation for hybridized surrogates like
 |                  NRHybSur3dq8.
 |      
 |      f_ref:      Frequency used to set the reference epoch at which the
 |                  reference frame is defined and the spins are specified.
 |                  See below for definition of the reference frame.
 |                  Should be in cycles/M if units = 'dimensionless', should be
 |                  in Hertz if units = 'mks'.
 |                  Default: If f_ref is not given, we set f_ref = f_low. If
 |                  f_low is 0, this corresponds to the initial index.
 |      
 |                  For time domain models, f_ref is used to determine a t_ref,
 |                  such that the orbital frequency in the coprecessing frame
 |                  equals f_ref/2 at t=t_ref.
 |      
 |      dt, df :    Time/Frequency step size, specify at most one of dt/df,
 |                  depending on whether the surrogate is a time/frequency domain
 |                  surrogate.
 |                  Default: None. If None, the internal domain of the surrogate is
 |                  used, which can be nonuniformly sampled.
 |                  dt (df) Should be in M (cycles/M) if units = 'dimensionless',
 |                  should be in seconds (Hertz) if units = 'mks'. Do not specify
 |                  times/freqs if using dt/df.
 |      
 |      
 |      times, freqs:
 |                  Array of time/frequency samples at which to evaluate the
 |                  waveform, depending on whether the surrogate is a
 |                  time/frequency domain surrogate. time (freqs) should be in
 |                  M (cycles/M) if units = 'dimensionless', should be in
 |                  seconds (Hertz) if units = 'mks'. Do not specify dt/df if
 |                  using times/freqs. Default None.
 |      
 |      ellMax:     Maximum ell index for modes to include. All available m
 |                  indicies for each ell will be included automatically. The 
 |                  m<0 modes will automatically be included for nonprecessing
 |                  models. 
 |                  Default: None, in which case all available ells will be
 |                  included. 
 |      
 |      mode_list : A list of (ell, m) modes tuples to be included. Valid only
 |                  for nonprecessing models.
 |      
 |                  Example: mode_list = [(2,2),(2,1)].
 |                  Default: None, in which case all available modes are included.
 |      
 |                  At most one of ellMax and mode_list can be specified.
 |      
 |                  Note: mode_list is allowed only for nonprecessing models; for
 |                  precessing models use ellMax. For precessing systems, all m
 |                  indices of a given ell index mix with each other, so there is
 |                  no clear hierarchy. To get the individual modes just don't
 |                  specify inclination and a dictionary of modes will be returned.
 |      
 |                  Note: When the inclination is set, the m<0 modes are 
 |                  automatically included. For example, passing mode_list = [(2,2)] 
 |                  will include the (2,2) and (2,-2) modes in the computation of
 |                  the strain.
 |      
 |                  Note: When the inclination is None, the m<0 modes are
 |                  automatically generated.
 |      
 |      inclination : Inclination angle between the orbital angular momentum
 |                  direction at the reference epoch and the line-of-sight to the
 |                  observer. If inclination is None, the mode data is returned
 |                  as a dictionary.
 |                  Default: None.
 |      
 |      phi_ref :   The azimuthal angle on the sky of the source frame following
 |                  the LAL convention.
 |                  Default: 0.
 |      
 |                  If inclination/phi_ref are specified, the complex strain (h =
 |                  hplus -i hcross) evaluated at (inclination, pi/2 - phi_ref) on
 |                  the sky of the reference frame is returned. This follows the
 |                  same convention as LAL. See below for definition of the
 |                  reference frame.
 |      
 |      precessing_opts:
 |                  A dictionary containing optional parameters for a precessing
 |                  surrogate model. Default: None.
 |                  Allowed keys are:
 |                  init_orbphase: The orbital phase in the coprecessing frame
 |                      at the reference epoch.
 |                      Default: 0, in which case the coorbital frame and
 |                      coprecessing frame are the same.
 |                  init_quat: The unit quaternion (length 4 vector) giving the
 |                      rotation from the coprecessing frame to the inertial frame
 |                      at the reference epoch.
 |                      Default: None, in which case the coprecessing frame is the
 |                      same as the inertial frame.
 |                  return_dynamics:
 |                      Return the frame dynamics and spin evolution along with
 |                      the waveform. Default: False.
 |                  Example: precessing_opts = {
 |                                      'init_orbphase': 0,
 |                                      'init_quat': [1,0,0,0],
 |                                      'return_dynamics': True
 |                                      }
 |      
 |      tidal_opts:
 |                  A dictionary containing optional parameters for a tidal
 |                  surrogate model. Default: None.
 |                  Allowed keys are:
 |                  Lambda1: The tidal deformability parameter for the heavier
 |                      object.
 |                  Lambda2: The tidal deformability parameter for the lighter
 |                      object.
 |                  Example: tidal_opts = {'Lambda1': 200, 'Lambda2': 300}
 |      
 |      
 |      par_dict:   A dictionary containing any additional parameters needed for a
 |                  particular surrogate model. Default: None.
 |      
 |      units:      'dimensionless' or 'mks'. Default: 'dimensionless'.
 |                  If 'dimensionless': Any of f_low, f_ref, dt, df, times and
 |                      freqs, if specified, must be in dimensionless units. That
 |                      is, dt/times should be in units of M, while f_ref, f_low
 |                      and df/freqs should be in units of cycles/M.
 |                      M and dist_mpc must be None. The waveform and domain are
 |                      returned as dimensionless quantities as well.
 |                  If 'mks': Any of f_low, f_ref, dt, df, times and freqs, if
 |                      specified, must be in MKS units. That is, dt/times should
 |                      be in seconds, while f_ref, f_low and df/freqs should be
 |                      in Hz. M and dist_mpc must be specified. The waveform and
 |                      domain are returned in MKS units as well.
 |      
 |      
 |      skip_param_checks :
 |                  Skip sanity checks for inputs. Use this if you want to
 |                  extrapolate outside allowed range. Default: False.
 |      
 |      taper_end_duration:
 |                  Taper the last TAPER_END_DURATION (M) of a time-domain waveform
 |                  in units of M. For exmple, passing 40 will taper the last 40M.
 |                  When set to None, no taper is applied
 |                  Default: None.
 |      
 |      RETURNS
 |      =====
 |      
 |      domain, h, dynamics
 |      
 |      
 |      domain :    Array of time/frequency samples corresponding to h and
 |                  dynamics, depending on whether the surrogate is a
 |                  time/frequency domain model. This is the same as times/freqs
 |                  if times/freqs are given as an inputs.
 |                  For time domain models the time is set to 0 at the peak of
 |                  the waveform. The time (frequency) values are in M (cycles/M)
 |                  if units = 'dimensionless', they are in seconds (Hertz) if
 |                  units = 'mks'
 |      
 |      h :         The waveform.
 |                      If inclination is specified, the complex strain (h = hplus
 |                      -i hcross) evaluated at (inclination, pi/2 - phi_ref) on
 |                      the sky of the reference frame is returned. This follows
 |                      the LAL convention, see below for details.  This includes
 |                      all modes given in the ellMax/mode_list argument. For
 |                      nonprecessing systems the m<0 modes are automatically
 |                      deduced from the m>0 modes. To see if a model is precessing
 |                      check self.keywords.
 |      
 |                      Else, h is a dictionary of available modes with (l, m)
 |                      tuples as keys. For example, h22 = h[(2,2)].
 |      
 |                      If M and dist_mpc are given, the physical waveform
 |                      at that distance is returned. Else, it is returned in
 |                      code units: r*h/M extrapolated to future null-infinity.
 |      
 |      dynamics:   A dict containing the frame dynamics and spin evolution. This
 |                  is None for nonprecessing models. This is also None if
 |                  return_dynamics in precessing_opts is False (Default).
 |      
 |                  The dynamics include (L=len(domain)):
 |      
 |                  q_copr = dynamics['q_copr']
 |                      The quaternion representing the coprecessing frame with
 |                      shape (4, L)
 |                  orbphase = dynamics['orbphase']
 |                      The orbital phase in the coprecessing frame with length L.
 |                  chiA = dynamics['chiA']
 |                      The inertial frame chiA with shape (L, 3)
 |                  chiB = dynamics['chiB']
 |                      The inertial frame chiB with shape (L, 3)
 |      
 |      
 |      IMPORTANT NOTES:
 |      ===============
 |      
 |      The reference frame (or inertial frame) is defined as follows:
 |          The +ve z-axis is along the orbital angular momentum at the reference
 |          epoch. The separation vector from the lighter BH to the heavier BH at
 |          the reference epoch is along the +ve x-axis. The y-axis completes the
 |          right-handed triad. The reference epoch is set using f_ref.
 |      
 |          Now, if inclination/phi_ref are given, the waveform is evaluated at
 |          (inclination, pi/2 - phi_ref) in the reference frame. This agrees with
 |          the LAL convention. See LIGO DCC document T1800226 for the LAL frame
 |          diagram.
 |  
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from SurrogateEvaluator:
 |  
 |  __dict__
 |      dictionary for instance variables
 |  
 |  __weakref__
 |      list of weak references to the object

Evaluate the waveform¶

Evaluate waveform modes in dimensionless units (default)¶

In [5]:

q = 4
chiA = [-0.2, 0.4, 0.1]
chiB = [-0.5, 0.2, -0.4]
dt = 0.1        # timestep size, Units of M
f_low = 0       # initial frequency, f_low=0 returns the full surrogate
t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low)   # dyn stands for dynamics, do dyn.keys() to see contents

q = 4 chiA = [-0.2, 0.4, 0.1] chiB = [-0.5, 0.2, -0.4] dt = 0.1 # timestep size, Units of M f_low = 0 # initial frequency, f_low=0 returns the full surrogate t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low) # dyn stands for dynamics, do dyn.keys() to see contents

In [6]:

# Let's see all available modes
print( sorted(h.keys()) )

# Let's see all available modes print( sorted(h.keys()) )

[(2, -2), (2, -1), (2, 0), (2, 1), (2, 2), (3, -3), (3, -2), (3, -1), (3, 0), (3, 1), (3, 2), (3, 3), (4, -4), (4, -3), (4, -2), (4, -1), (4, 0), (4, 1), (4, 2), (4, 3), (4, 4), (5, -5), (5, -4), (5, -3), (5, -2), (5, -1), (5, 0), (5, 1), (5, 2), (5, 3), (5, 4), (5, 5)]

In [7]:

P.plot(t, h[(2,2)].real, label='l2m2 real')
P.plot(t, h[(2,1)].real, label='l2m1 real')
P.plot(t, h[(3,3)].real, label='l3m3 real')
P.ylabel('Re[$h_{lm}$]', fontsize=18)
P.xlabel('t [M]', fontsize=18)
P.legend()

P.plot(t, h[(2,2)].real, label='l2m2 real') P.plot(t, h[(2,1)].real, label='l2m1 real') P.plot(t, h[(3,3)].real, label='l3m3 real') P.ylabel('Re[$h_{lm}$]', fontsize=18) P.xlabel('t [M]', fontsize=18) P.legend()

Out[7]:

<matplotlib.legend.Legend at 0x7a165d7a4290>

Evaluate waveform modes in physical units¶

In [8]:

q = 4
chiA = [-0.2, 0.4, 0.1]
chiB = [-0.5, 0.2, -0.4]
f_ref = 20         # Reference frequecny in Hz. The spins are assumed to specified at this frequency
f_low = 0          # initial frequency, f_low=0 returns the full surrogate
M = 70             # Total masss in solar masses
dist_mpc = 100     # distance in megaparsecs
dt = 1./4096       # step size in seconds
ellMax = 4         # Highest ell index for modes to use

# dyn stands for dynamics, do dyn.keys() to see contents
t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc, units='mks')

P.plot(t, h[(2,2)].real, label='l2m2 real')
P.plot(t, h[(2,1)].real, label='l2m1 real')
P.ylabel('Re[$h_{lm}$]', fontsize=18)
P.xlabel('t [s]', fontsize=18)
P.legend()

q = 4 chiA = [-0.2, 0.4, 0.1] chiB = [-0.5, 0.2, -0.4] f_ref = 20 # Reference frequecny in Hz. The spins are assumed to specified at this frequency f_low = 0 # initial frequency, f_low=0 returns the full surrogate M = 70 # Total masss in solar masses dist_mpc = 100 # distance in megaparsecs dt = 1./4096 # step size in seconds ellMax = 4 # Highest ell index for modes to use # dyn stands for dynamics, do dyn.keys() to see contents t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc, units='mks') P.plot(t, h[(2,2)].real, label='l2m2 real') P.plot(t, h[(2,1)].real, label='l2m1 real') P.ylabel('Re[$h_{lm}$]', fontsize=18) P.xlabel('t [s]', fontsize=18) P.legend()

Out[8]:

<matplotlib.legend.Legend at 0x7a165d6a3150>

Evaluate waveform at a point on the sky¶

In [9]:

q = 4
chiA = [-0.2, 0.4, 0.1]
chiB = [-0.5, 0.2, -0.4]
f_ref = 20         # Reference frequecny in Hz. The spins are assumed to specified at this frequency
f_low = 0          # initial frequency, f_low=0 returns the full surrogate
M = 70             # Total masss in solar masses
dist_mpc = 100     # distance in megaparsecs
dt = 1./4096       # step size in seconds
ellMax = 4         # Highest ell index for modes to use
inclination = np.pi/4
phi_ref = np.pi/5

# Will only include modes ell<=ellMax
# Returns h = h_+ -i h_x at (inclination, phi_ref) in the sky of the source frame
# dyn stands for dynamics, do dyn.keys() to see contents
t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc, 
           inclination=inclination, phi_ref=phi_ref, units='mks')

P.plot(t, h.real)
P.ylabel('$h_{+}$ $(\iota, \phi_{ref})$', fontsize=18)
P.xlabel('t [s]', fontsize=18)

q = 4 chiA = [-0.2, 0.4, 0.1] chiB = [-0.5, 0.2, -0.4] f_ref = 20 # Reference frequecny in Hz. The spins are assumed to specified at this frequency f_low = 0 # initial frequency, f_low=0 returns the full surrogate M = 70 # Total masss in solar masses dist_mpc = 100 # distance in megaparsecs dt = 1./4096 # step size in seconds ellMax = 4 # Highest ell index for modes to use inclination = np.pi/4 phi_ref = np.pi/5 # Will only include modes ell<=ellMax # Returns h = h_+ -i h_x at (inclination, phi_ref) in the sky of the source frame # dyn stands for dynamics, do dyn.keys() to see contents t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc, inclination=inclination, phi_ref=phi_ref, units='mks') P.plot(t, h.real) P.ylabel('$h_{+}$ $(\iota, \phi_{ref})$', fontsize=18) P.xlabel('t [s]', fontsize=18)

Out[9]:

Text(0.5, 0, 't [s]')

Dynamics Surrogate¶

In [10]:

q = 4
chiA = [-0.2, 0.4, 0.1]  # unless f_ref is given the spins are assumed to be given at the start of the waveform
chiB = [-0.5, 0.2, -0.4]
dt = 0.1        # step size, Units of M
f_low = 0       # initial frequency, f_low=0 returns the full surrogate
t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low, precessing_opts={'return_dynamics': True})        # dyn stands for dynamics, do dyn.keys() to see contents

q = 4 chiA = [-0.2, 0.4, 0.1] # unless f_ref is given the spins are assumed to be given at the start of the waveform chiB = [-0.5, 0.2, -0.4] dt = 0.1 # step size, Units of M f_low = 0 # initial frequency, f_low=0 returns the full surrogate t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low, precessing_opts={'return_dynamics': True}) # dyn stands for dynamics, do dyn.keys() to see contents

In [11]:

# See all available dynamics data
dyn.keys()

# See all available dynamics data dyn.keys()

Out[11]:

dict_keys(['chiA', 'chiB', 'chiA_copr', 'chiB_copr', 'q_copr', 'orbphase'])

In [12]:

P.figure(1)
P.plot(t, dyn['chiA'][:,0], label='$\chi_{Ax}$')
P.plot(t, dyn['chiA'][:,1], label='$\chi_{Ay}$')
P.plot(t, dyn['chiA'][:,2], label='$\chi_{Az}$')
P.ylabel('$\chi_A$', fontsize=18)
P.xlabel('t', fontsize=18)
P.title('Spin of heavier BH')
P.legend(fontsize=14)

P.figure(2)
P.plot(t, dyn['chiB'][:,0], label='$\chi_{Bx}$')
P.plot(t, dyn['chiB'][:,1], label='$\chi_{By}$')
P.plot(t, dyn['chiB'][:,2], label='$\chi_{Bz}$')
P.ylabel('$\chi_B$', fontsize=18)
P.xlabel('t', fontsize=18)
P.title('Spin of lighter BH')
P.legend(fontsize=14)

P.figure(3)
P.plot(t, dyn['orbphase'][:,])
P.ylabel('$\phi_{\mathrm{orb}}$', fontsize=18)
P.xlabel('t', fontsize=18)
P.title('Orbital phase')

P.figure(4)
P.plot(t, dyn['q_copr'][0,:], label='$\hat{Q}_0$')
P.plot(t, dyn['q_copr'][1,:], label='$\hat{Q}_1$')
P.plot(t, dyn['q_copr'][2,:], label='$\hat{Q}_2$')
P.plot(t, dyn['q_copr'][3,:], label='$\hat{Q}_3$')
P.ylabel('$\hat{Q}$', fontsize=18)
P.xlabel('t', fontsize=18)
P.title('Coprecessing frame quaternions')
P.legend(fontsize=14)

P.figure(1) P.plot(t, dyn['chiA'][:,0], label='$\chi_{Ax}$') P.plot(t, dyn['chiA'][:,1], label='$\chi_{Ay}$') P.plot(t, dyn['chiA'][:,2], label='$\chi_{Az}$') P.ylabel('$\chi_A$', fontsize=18) P.xlabel('t', fontsize=18) P.title('Spin of heavier BH') P.legend(fontsize=14) P.figure(2) P.plot(t, dyn['chiB'][:,0], label='$\chi_{Bx}$') P.plot(t, dyn['chiB'][:,1], label='$\chi_{By}$') P.plot(t, dyn['chiB'][:,2], label='$\chi_{Bz}$') P.ylabel('$\chi_B$', fontsize=18) P.xlabel('t', fontsize=18) P.title('Spin of lighter BH') P.legend(fontsize=14) P.figure(3) P.plot(t, dyn['orbphase'][:,]) P.ylabel('$\phi_{\mathrm{orb}}$', fontsize=18) P.xlabel('t', fontsize=18) P.title('Orbital phase') P.figure(4) P.plot(t, dyn['q_copr'][0,:], label='$\hat{Q}_0$') P.plot(t, dyn['q_copr'][1,:], label='$\hat{Q}_1$') P.plot(t, dyn['q_copr'][2,:], label='$\hat{Q}_2$') P.plot(t, dyn['q_copr'][3,:], label='$\hat{Q}_3$') P.ylabel('$\hat{Q}$', fontsize=18) P.xlabel('t', fontsize=18) P.title('Coprecessing frame quaternions') P.legend(fontsize=14)

Out[12]:

<matplotlib.legend.Legend at 0x7a165c58dcd0>

Going Further: Comparison between NRSur7dq4 and SEOBNRv4PHMSur¶

NOTE: to compute matches, we need PyCBC. It can installed via pip using

>>> pip install pycbc

In [13]:

import seaborn as sns
sns.set_theme(style="whitegrid")
sns.color_palette("colorblind")
from seaborn.palettes import SEABORN_PALETTES
colors = SEABORN_PALETTES["colorblind"]

from gwtools.gwtools import minimize_norm_error_multi_mode

import matplotlib.pyplot as P
import matplotlib.pylab as plt
%matplotlib inline

import seaborn as sns sns.set_theme(style="whitegrid") sns.color_palette("colorblind") from seaborn.palettes import SEABORN_PALETTES colors = SEABORN_PALETTES["colorblind"] from gwtools.gwtools import minimize_norm_error_multi_mode import matplotlib.pyplot as P import matplotlib.pylab as plt %matplotlib inline

In [14]:

eob_sur = gwsurrogate.LoadSurrogate('SEOBNRv4PHMSur')
nr_sur = gwsurrogate.LoadSurrogate('NRSur7dq4')

eob_sur = gwsurrogate.LoadSurrogate('SEOBNRv4PHMSur') nr_sur = gwsurrogate.LoadSurrogate('NRSur7dq4')

Loaded SEOBNRv4PHMSur model
Loaded NRSur7dq4 model

Aligned-spin case¶

Here we compare the EOB and NR surrogate models for an aligned-spin system. We will show a time shift and z-axis roation are needed to match the two waveforms

In [15]:

q = 1.5
chiA = [0.0, 0.0, 0.0]
chiB = [0.0, 0.0, 0.0]

dt = 0.1        # timestep size, Units of M
f_low = 0       # initial frequency, f_low=0 returns the full surrogate
f_ref = 0.008

t_nr, h_nr, dyn_nr = nr_sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, precessing_opts = {'return_dynamics': True} )
t_eob, h_eob, dyn_eob = eob_sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, precessing_opts = {'return_dynamics': True}) 

# Requested list of modes to match over
mode_list = [(2, 2), (2, 1), (3, 3), (4, 4)]

hnr_subset   = {mode: h_nr[mode] for mode in mode_list}
heob_subset  = {mode: h_eob[mode] for mode in mode_list}

[errors_before_min, errors_after_min], [tc, phic], [common_times,h_eob_aligned,h_nr_aligned] = \
    minimize_norm_error_multi_mode(t_eob, heob_subset, t_nr, hnr_subset, mode_list=None, ell_m_match=(2,2),\
    t_low_adj=0,t_up_adj=0,method='nelder-mead')
print("time-domain L2 difference = %e"%errors_after_min)

q = 1.5 chiA = [0.0, 0.0, 0.0] chiB = [0.0, 0.0, 0.0] dt = 0.1 # timestep size, Units of M f_low = 0 # initial frequency, f_low=0 returns the full surrogate f_ref = 0.008 t_nr, h_nr, dyn_nr = nr_sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, precessing_opts = {'return_dynamics': True} ) t_eob, h_eob, dyn_eob = eob_sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, precessing_opts = {'return_dynamics': True}) # Requested list of modes to match over mode_list = [(2, 2), (2, 1), (3, 3), (4, 4)] hnr_subset = {mode: h_nr[mode] for mode in mode_list} heob_subset = {mode: h_eob[mode] for mode in mode_list} [errors_before_min, errors_after_min], [tc, phic], [common_times,h_eob_aligned,h_nr_aligned] = \ minimize_norm_error_multi_mode(t_eob, heob_subset, t_nr, hnr_subset, mode_list=None, ell_m_match=(2,2),\ t_low_adj=0,t_up_adj=0,method='nelder-mead') print("time-domain L2 difference = %e"%errors_after_min)

time-domain L2 difference = 9.033986e-03

In [16]:

P.figure(1)
P.title("Before time/phase alginment")
P.plot(t_eob,heob_subset[(2,2)],'b')
P.plot(t_nr,hnr_subset[(2,2)],'r--')
P.xlim([-800, 100])

P.figure(2)
P.title("After time/phase alginment")
P.plot(common_times,h_eob_aligned[(2,2)],'b')
P.plot(common_times,h_nr_aligned[(2,2)],'r--')
P.xlim([-800, 100])

P.figure(1) P.title("Before time/phase alginment") P.plot(t_eob,heob_subset[(2,2)],'b') P.plot(t_nr,hnr_subset[(2,2)],'r--') P.xlim([-800, 100]) P.figure(2) P.title("After time/phase alginment") P.plot(common_times,h_eob_aligned[(2,2)],'b') P.plot(common_times,h_nr_aligned[(2,2)],'r--') P.xlim([-800, 100])

/home/balzani57/pool/anaconda3/envs/py311-gw-5-2024-dev/lib/python3.11/site-packages/matplotlib/cbook.py:1699: ComplexWarning: Casting complex values to real discards the imaginary part
  return math.isfinite(val)
/home/balzani57/pool/anaconda3/envs/py311-gw-5-2024-dev/lib/python3.11/site-packages/matplotlib/cbook.py:1345: ComplexWarning: Casting complex values to real discards the imaginary part
  return np.asarray(x, float)

Out[16]:

(-800.0, 100.0)

Evaluate waveform modes in dimensionless units (default)¶

Now we consider precessing systems

In [30]:

# Some parameters to match waveforms for

# Remember:
## NRsur extrapolated for 4 < q < 6
## Both models extrapolate for |chi_1,2| > 0.8
## Uncomment a single case to be used

# Within the domain of both models
q = 3.9
chiA = [-0.2, 0.4, 0.3]
chiB = [-0.1, 0.5, -0.1]

## Extrapolate q
# q = 5.9
# chiA = [-0.2, 0.4, 0.1]
# chiB = [-0.5, 0.2, -0.2]

## Extrapolate spins
# q = 3.9
# chiA = [-0.2, 0.6, 0.7]
# chiB = [-0.5, 0.2, -0.7]

##extrapolate spins and q
# q = 5.9
# chiA = [-0.2, 0.6, 0.7]
# chiB = [-0.5, 0.2, -0.7]

# Some parameters to match waveforms for # Remember: ## NRsur extrapolated for 4 < q < 6 ## Both models extrapolate for |chi_1,2| > 0.8 ## Uncomment a single case to be used # Within the domain of both models q = 3.9 chiA = [-0.2, 0.4, 0.3] chiB = [-0.1, 0.5, -0.1] ## Extrapolate q # q = 5.9 # chiA = [-0.2, 0.4, 0.1] # chiB = [-0.5, 0.2, -0.2] ## Extrapolate spins # q = 3.9 # chiA = [-0.2, 0.6, 0.7] # chiB = [-0.5, 0.2, -0.7] ##extrapolate spins and q # q = 5.9 # chiA = [-0.2, 0.6, 0.7] # chiB = [-0.5, 0.2, -0.7]

In [31]:

dt = 0.1        # timestep size, Units of M
f_low = 0       # initial frequency, f_low=0 returns the full surrogate
f_ref = 0.008

t_nr, h_nr, dyn_nr = nr_sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, precessing_opts = {'return_dynamics': True} )
t_eob, h_eob, dyn_eob = eob_sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, precessing_opts = {'return_dynamics': True}) 

dt = 0.1 # timestep size, Units of M f_low = 0 # initial frequency, f_low=0 returns the full surrogate f_ref = 0.008 t_nr, h_nr, dyn_nr = nr_sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, precessing_opts = {'return_dynamics': True} ) t_eob, h_eob, dyn_eob = eob_sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, precessing_opts = {'return_dynamics': True})

In [32]:

plt.figure(figsize=(20, 6))

plt.subplot(1,2,1)
P.plot(t_nr, h_nr[(2,2)].real, label='l2m2 real: NR',c=colors[0])
P.plot(t_nr, h_nr[(2,1)].real, label='l2m1 real: NR',c=colors[1])
# P.plot(t_nr, h_nr[(3,3)].real, label='l3m3 real: NR',c=colors[2])

P.plot(t_eob, h_eob[(2,2)].real, label='l2m2 real: EOB', c=colors[0], ls='dashed')
P.plot(t_eob, h_eob[(2,1)].real, label='l2m1 real: EOB',c=colors[1], ls='dashed')
# P.plot(t_eob, h_eob[(3,3)].real, label='l3m3 real: EOB',c=colors[2], ls='dashed')
P.ylabel('Re[$h_{lm}$]', fontsize=18)
P.xlabel('t [M]', fontsize=18)


plt.subplot(1,2,2)
P.plot(t_nr, h_nr[(2,2)].real, label='l2m2 real: NR',c=colors[0])
P.plot(t_nr, h_nr[(2,1)].real, label='l2m1 real: NR',c=colors[1])
# P.plot(t_nr, h_nr[(3,3)].real, label='l3m3 real: NR',c=colors[2])

P.plot(t_eob, h_eob[(2,2)].real, label='l2m2 real: EOB', c=colors[0], ls='dashed')
P.plot(t_eob, h_eob[(2,1)].real, label='l2m1 real: EOB',c=colors[1], ls='dashed')
# P.plot(t_eob, h_eob[(3,3)].real, label='l3m3 real: EOB',c=colors[2], ls='dashed')
P.xlabel('t [M]', fontsize=18)
P.xlim(-100, 50)

P.legend()
plt.show()

plt.figure(figsize=(20, 6)) plt.subplot(1,2,1) P.plot(t_nr, h_nr[(2,2)].real, label='l2m2 real: NR',c=colors[0]) P.plot(t_nr, h_nr[(2,1)].real, label='l2m1 real: NR',c=colors[1]) # P.plot(t_nr, h_nr[(3,3)].real, label='l3m3 real: NR',c=colors[2]) P.plot(t_eob, h_eob[(2,2)].real, label='l2m2 real: EOB', c=colors[0], ls='dashed') P.plot(t_eob, h_eob[(2,1)].real, label='l2m1 real: EOB',c=colors[1], ls='dashed') # P.plot(t_eob, h_eob[(3,3)].real, label='l3m3 real: EOB',c=colors[2], ls='dashed') P.ylabel('Re[$h_{lm}$]', fontsize=18) P.xlabel('t [M]', fontsize=18) plt.subplot(1,2,2) P.plot(t_nr, h_nr[(2,2)].real, label='l2m2 real: NR',c=colors[0]) P.plot(t_nr, h_nr[(2,1)].real, label='l2m1 real: NR',c=colors[1]) # P.plot(t_nr, h_nr[(3,3)].real, label='l3m3 real: NR',c=colors[2]) P.plot(t_eob, h_eob[(2,2)].real, label='l2m2 real: EOB', c=colors[0], ls='dashed') P.plot(t_eob, h_eob[(2,1)].real, label='l2m1 real: EOB',c=colors[1], ls='dashed') # P.plot(t_eob, h_eob[(3,3)].real, label='l3m3 real: EOB',c=colors[2], ls='dashed') P.xlabel('t [M]', fontsize=18) P.xlim(-100, 50) P.legend() plt.show()

In [33]:

for k in h_nr:
    plt.figure(figsize=(20, 6))
    plt.subplot(1,2,1)
    P.plot(t_nr, h_nr[k].real, label='l{}m{} real: NR'.format(k[0], k[1]), c=colors[0])
    P.plot(t_nr, h_nr[k].imag, label='l{}m{} imag: NR'.format(k[0], k[1]), c=colors[0], ls='dashed')

    P.plot(t_eob, h_eob[k].real, label='l{}m{} real: EOB'.format(k[0], k[1]), c=colors[1])
    P.plot(t_eob, h_eob[k].imag, label='l{}m{} imag: EOB'.format(k[0], k[1]), c=colors[1], ls='dashed')
   
    P.ylabel('Re[$h_{lm}$]', fontsize=18)
    P.xlabel('t [M]', fontsize=18)


    plt.subplot(1,2,2)
    P.plot(t_nr, h_nr[k].real, label='l{}m{} real: NR'.format(k[0], k[1]), c=colors[0])
    P.plot(t_nr, h_nr[k].imag, label='l{}m{} imag: NR'.format(k[0], k[1]), c=colors[0], ls='dashed')

    P.plot(t_eob, h_eob[k].real, label='l{}m{} real: EOB'.format(k[0], k[1]), c=colors[1])
    P.plot(t_eob, h_eob[k].imag, label='l{}m{} imag: EOB'.format(k[0], k[1]), c=colors[1], ls='dashed')
   
    P.ylabel('Re[$h_{lm}$]', fontsize=18)
    P.xlabel('t [M]', fontsize=18)
    P.xlim(-100, 50)

    P.legend()
    plt.show()

for k in h_nr: plt.figure(figsize=(20, 6)) plt.subplot(1,2,1) P.plot(t_nr, h_nr[k].real, label='l{}m{} real: NR'.format(k[0], k[1]), c=colors[0]) P.plot(t_nr, h_nr[k].imag, label='l{}m{} imag: NR'.format(k[0], k[1]), c=colors[0], ls='dashed') P.plot(t_eob, h_eob[k].real, label='l{}m{} real: EOB'.format(k[0], k[1]), c=colors[1]) P.plot(t_eob, h_eob[k].imag, label='l{}m{} imag: EOB'.format(k[0], k[1]), c=colors[1], ls='dashed') P.ylabel('Re[$h_{lm}$]', fontsize=18) P.xlabel('t [M]', fontsize=18) plt.subplot(1,2,2) P.plot(t_nr, h_nr[k].real, label='l{}m{} real: NR'.format(k[0], k[1]), c=colors[0]) P.plot(t_nr, h_nr[k].imag, label='l{}m{} imag: NR'.format(k[0], k[1]), c=colors[0], ls='dashed') P.plot(t_eob, h_eob[k].real, label='l{}m{} real: EOB'.format(k[0], k[1]), c=colors[1]) P.plot(t_eob, h_eob[k].imag, label='l{}m{} imag: EOB'.format(k[0], k[1]), c=colors[1], ls='dashed') P.ylabel('Re[$h_{lm}$]', fontsize=18) P.xlabel('t [M]', fontsize=18) P.xlim(-100, 50) P.legend() plt.show()

In [34]:

P.figure(1,  figsize=(10, 6))
P.plot(t_nr, dyn_nr['chiA'][:,0], label='$\chi_{Ax}$: NR', c=colors[0])
P.plot(t_nr, dyn_nr['chiA'][:,1], label='$\chi_{Ay}$: NR', c=colors[1])
P.plot(t_nr, dyn_nr['chiA'][:,2], label='$\chi_{Az}$: NR', c=colors[2])

P.plot(t_eob, dyn_eob['chiA'][:,0], label='$\chi_{Ax}$: EOB', c=colors[0], ls='dashed')
P.plot(t_eob, dyn_eob['chiA'][:,1], label='$\chi_{Ay}$: EOB', c=colors[1], ls='dashed')
P.plot(t_eob, dyn_eob['chiA'][:,2], label='$\chi_{Az}$: EOB', c=colors[2], ls='dashed')

P.ylabel('$\chi_A$', fontsize=18)
P.xlabel('t', fontsize=18)
P.title('Spin of heavier BH')
P.legend(fontsize=14)

P.figure(2,figsize=(10, 6))
P.plot(t_nr, dyn_nr['chiB'][:,0], label='$\chi_{Bx}$: NR',c=colors[3])
P.plot(t_nr, dyn_nr['chiB'][:,1], label='$\chi_{By}$: NR',c=colors[4])
P.plot(t_nr, dyn_nr['chiB'][:,2], label='$\chi_{Bz}$: NR',c=colors[5])

P.plot(t_eob, dyn_eob['chiB'][:,0], label='$\chi_{Bx}$: EOB',c=colors[3], ls='dashed')
P.plot(t_eob, dyn_eob['chiB'][:,1], label='$\chi_{By}$: EOB',c=colors[4], ls='dashed')
P.plot(t_eob, dyn_eob['chiB'][:,2], label='$\chi_{Bz}$: EOB',c=colors[5], ls='dashed')

P.ylabel('$\chi_B$', fontsize=18)
P.xlabel('t', fontsize=18)
P.title('Spin of lighter BH')
P.legend(fontsize=14)

P.figure(3,figsize=(10, 6))
P.plot(t_nr, dyn_nr['orbphase'][:,],c=colors[0],label='NR')
P.plot(t_eob, dyn_eob['orbphase'][:,],c=colors[0],ls='dashed', label='EOB')
P.ylabel('$\phi_{\mathrm{orb}}$', fontsize=18)
P.xlabel('t', fontsize=18)
P.title('Orbital phase')

P.figure(4,figsize=(10, 6))
P.plot(t_nr, dyn_nr['q_copr'][0,:], label='$\hat{Q}_0$: NR',c=colors[6])
P.plot(t_nr, dyn_nr['q_copr'][1,:], label='$\hat{Q}_1$: NR',c=colors[7])
P.plot(t_nr, dyn_nr['q_copr'][2,:], label='$\hat{Q}_2$: NR',c=colors[8])
P.plot(t_nr, dyn_nr['q_copr'][3,:], label='$\hat{Q}_3$: NR',c=colors[9])

P.plot(t_eob, dyn_eob['q_copr'][0,:], label='$\hat{Q}_0$: EOB',c=colors[6], ls='dashed')
P.plot(t_eob, dyn_eob['q_copr'][1,:], label='$\hat{Q}_1$: EOB',c=colors[7], ls='dashed')
P.plot(t_eob, dyn_eob['q_copr'][2,:], label='$\hat{Q}_2$: EOB',c=colors[8], ls='dashed')
P.plot(t_eob, dyn_eob['q_copr'][3,:], label='$\hat{Q}_3$: EOB',c=colors[9], ls='dashed')
P.ylabel('$\hat{Q}$', fontsize=18)
P.xlabel('t', fontsize=18)
P.title('Coprecessing frame quaternions')
P.legend(fontsize=14)

P.figure(1, figsize=(10, 6)) P.plot(t_nr, dyn_nr['chiA'][:,0], label='$\chi_{Ax}$: NR', c=colors[0]) P.plot(t_nr, dyn_nr['chiA'][:,1], label='$\chi_{Ay}$: NR', c=colors[1]) P.plot(t_nr, dyn_nr['chiA'][:,2], label='$\chi_{Az}$: NR', c=colors[2]) P.plot(t_eob, dyn_eob['chiA'][:,0], label='$\chi_{Ax}$: EOB', c=colors[0], ls='dashed') P.plot(t_eob, dyn_eob['chiA'][:,1], label='$\chi_{Ay}$: EOB', c=colors[1], ls='dashed') P.plot(t_eob, dyn_eob['chiA'][:,2], label='$\chi_{Az}$: EOB', c=colors[2], ls='dashed') P.ylabel('$\chi_A$', fontsize=18) P.xlabel('t', fontsize=18) P.title('Spin of heavier BH') P.legend(fontsize=14) P.figure(2,figsize=(10, 6)) P.plot(t_nr, dyn_nr['chiB'][:,0], label='$\chi_{Bx}$: NR',c=colors[3]) P.plot(t_nr, dyn_nr['chiB'][:,1], label='$\chi_{By}$: NR',c=colors[4]) P.plot(t_nr, dyn_nr['chiB'][:,2], label='$\chi_{Bz}$: NR',c=colors[5]) P.plot(t_eob, dyn_eob['chiB'][:,0], label='$\chi_{Bx}$: EOB',c=colors[3], ls='dashed') P.plot(t_eob, dyn_eob['chiB'][:,1], label='$\chi_{By}$: EOB',c=colors[4], ls='dashed') P.plot(t_eob, dyn_eob['chiB'][:,2], label='$\chi_{Bz}$: EOB',c=colors[5], ls='dashed') P.ylabel('$\chi_B$', fontsize=18) P.xlabel('t', fontsize=18) P.title('Spin of lighter BH') P.legend(fontsize=14) P.figure(3,figsize=(10, 6)) P.plot(t_nr, dyn_nr['orbphase'][:,],c=colors[0],label='NR') P.plot(t_eob, dyn_eob['orbphase'][:,],c=colors[0],ls='dashed', label='EOB') P.ylabel('$\phi_{\mathrm{orb}}$', fontsize=18) P.xlabel('t', fontsize=18) P.title('Orbital phase') P.figure(4,figsize=(10, 6)) P.plot(t_nr, dyn_nr['q_copr'][0,:], label='$\hat{Q}_0$: NR',c=colors[6]) P.plot(t_nr, dyn_nr['q_copr'][1,:], label='$\hat{Q}_1$: NR',c=colors[7]) P.plot(t_nr, dyn_nr['q_copr'][2,:], label='$\hat{Q}_2$: NR',c=colors[8]) P.plot(t_nr, dyn_nr['q_copr'][3,:], label='$\hat{Q}_3$: NR',c=colors[9]) P.plot(t_eob, dyn_eob['q_copr'][0,:], label='$\hat{Q}_0$: EOB',c=colors[6], ls='dashed') P.plot(t_eob, dyn_eob['q_copr'][1,:], label='$\hat{Q}_1$: EOB',c=colors[7], ls='dashed') P.plot(t_eob, dyn_eob['q_copr'][2,:], label='$\hat{Q}_2$: EOB',c=colors[8], ls='dashed') P.plot(t_eob, dyn_eob['q_copr'][3,:], label='$\hat{Q}_3$: EOB',c=colors[9], ls='dashed') P.ylabel('$\hat{Q}$', fontsize=18) P.xlabel('t', fontsize=18) P.title('Coprecessing frame quaternions') P.legend(fontsize=14)

Out[34]:

<matplotlib.legend.Legend at 0x7a1639dd8a10>

Evaluate waveform modes in physical units¶

In [35]:

f_ref = 25        # Reference frequecny in Hz. The spins are assumed to specified at this frequency
f_low = 25          # initial frequency, f_low=0 returns the full surrogate
M = 70             # Total masss in solar masses
dist_mpc = 100     # distance in megaparsecs
dt = 1./4096       # step size in seconds
ellMax = 4         # Highest ell index for modes to use

# dyn stands for dynamics, do dyn.keys() to see contents
# t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc, units='mks')

t_nr, h_nr, dyn_nr = nr_sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc, units='mks')
t_eob, h_eob, dyn_eob = eob_sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc, units='mks') 

P.figure(figsize=(20, 6))

plt.subplot(1,2,1)
P.plot(t_nr, h_nr[(2,2)].real, label='l2m2 real', c=colors[0])
P.plot(t_nr, h_nr[(2,1)].real, label='l2m1 real', c=colors[1])

P.plot(t_eob, h_eob[(2,2)].real, label='l2m2 real', ls='dashed', c=colors[0])
P.plot(t_eob, h_eob[(2,1)].real, label='l2m1 real', ls='dashed', c=colors[1])

P.ylabel('Re[$h_{lm}$]', fontsize=18)
P.xlabel('t [s]', fontsize=18)
P.legend()

plt.subplot(1,2,2)
P.plot(t_nr, h_nr[(2,2)].real, label='l2m2 real', c=colors[0])
P.plot(t_nr, h_nr[(2,1)].real, label='l2m1 real', c=colors[1])

P.plot(t_eob, h_eob[(2,2)].real, label='l2m2 real', ls='dashed', c=colors[0])
P.plot(t_eob, h_eob[(2,1)].real, label='l2m1 real', ls='dashed', c=colors[1])

P.ylabel('Re[$h_{lm}$]', fontsize=18)
P.xlabel('t [s]', fontsize=18)
P.legend()
P.xlim(-0.1,0.05)

f_ref = 25 # Reference frequecny in Hz. The spins are assumed to specified at this frequency f_low = 25 # initial frequency, f_low=0 returns the full surrogate M = 70 # Total masss in solar masses dist_mpc = 100 # distance in megaparsecs dt = 1./4096 # step size in seconds ellMax = 4 # Highest ell index for modes to use # dyn stands for dynamics, do dyn.keys() to see contents # t, h, dyn = sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc, units='mks') t_nr, h_nr, dyn_nr = nr_sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc, units='mks') t_eob, h_eob, dyn_eob = eob_sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc, units='mks') P.figure(figsize=(20, 6)) plt.subplot(1,2,1) P.plot(t_nr, h_nr[(2,2)].real, label='l2m2 real', c=colors[0]) P.plot(t_nr, h_nr[(2,1)].real, label='l2m1 real', c=colors[1]) P.plot(t_eob, h_eob[(2,2)].real, label='l2m2 real', ls='dashed', c=colors[0]) P.plot(t_eob, h_eob[(2,1)].real, label='l2m1 real', ls='dashed', c=colors[1]) P.ylabel('Re[$h_{lm}$]', fontsize=18) P.xlabel('t [s]', fontsize=18) P.legend() plt.subplot(1,2,2) P.plot(t_nr, h_nr[(2,2)].real, label='l2m2 real', c=colors[0]) P.plot(t_nr, h_nr[(2,1)].real, label='l2m1 real', c=colors[1]) P.plot(t_eob, h_eob[(2,2)].real, label='l2m2 real', ls='dashed', c=colors[0]) P.plot(t_eob, h_eob[(2,1)].real, label='l2m1 real', ls='dashed', c=colors[1]) P.ylabel('Re[$h_{lm}$]', fontsize=18) P.xlabel('t [s]', fontsize=18) P.legend() P.xlim(-0.1,0.05)

Out[35]:

(-0.1, 0.05)

Evaluate waveform at a point on the sky¶

In [36]:

M = 70             # Total masss in solar masses
dist_mpc = 100     # distance in megaparsecs
dt = 1./4096       # step size in seconds
ellMax = 4         # Highest ell index for modes to use
inclination = 2.*np.pi/3
phi_ref = np.pi/5

# Will only include modes ell<=ellMax
# Returns h = h_+ -i h_x at (inclination, phi_ref) in the sky of the source frame
# dyn stands for dynamics, do dyn.keys() to see contents
t_nr, h_nr, dyn = nr_sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc, 
           inclination=inclination, phi_ref=phi_ref, units='mks')

t_eob, h_eob, dyn = eob_sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc, 
           inclination=inclination, phi_ref=phi_ref, units='mks')

P.figure(figsize=(20, 6))

plt.subplot(1,2,1)
P.plot(t_nr, h_nr.real, c=colors[0], label='NR')
P.plot(t_eob, h_eob.real, c=colors[1], ls='dashed', label='EOB')
P.ylabel('$h_{+}$ $(\iota, \phi_{ref})$', fontsize=18)
P.xlabel('t [s]', fontsize=18)
P.legend()

plt.subplot(1,2,2)
P.plot(t_nr, h_nr.real, c=colors[0], label='NR')
P.plot(t_eob, h_eob.real, c=colors[1], ls='dashed', label='EOB')
P.ylabel('$h_{+}$ $(\iota, \phi_{ref})$', fontsize=18)
P.xlabel('t [s]', fontsize=18)
P.legend()
P.xlim(-0.1,0.05)

M = 70 # Total masss in solar masses dist_mpc = 100 # distance in megaparsecs dt = 1./4096 # step size in seconds ellMax = 4 # Highest ell index for modes to use inclination = 2.*np.pi/3 phi_ref = np.pi/5 # Will only include modes ell<=ellMax # Returns h = h_+ -i h_x at (inclination, phi_ref) in the sky of the source frame # dyn stands for dynamics, do dyn.keys() to see contents t_nr, h_nr, dyn = nr_sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc, inclination=inclination, phi_ref=phi_ref, units='mks') t_eob, h_eob, dyn = eob_sur(q, chiA, chiB, dt=dt, f_low=f_low, f_ref=f_ref, ellMax=ellMax, M=M, dist_mpc=dist_mpc, inclination=inclination, phi_ref=phi_ref, units='mks') P.figure(figsize=(20, 6)) plt.subplot(1,2,1) P.plot(t_nr, h_nr.real, c=colors[0], label='NR') P.plot(t_eob, h_eob.real, c=colors[1], ls='dashed', label='EOB') P.ylabel('$h_{+}$ $(\iota, \phi_{ref})$', fontsize=18) P.xlabel('t [s]', fontsize=18) P.legend() plt.subplot(1,2,2) P.plot(t_nr, h_nr.real, c=colors[0], label='NR') P.plot(t_eob, h_eob.real, c=colors[1], ls='dashed', label='EOB') P.ylabel('$h_{+}$ $(\iota, \phi_{ref})$', fontsize=18) P.xlabel('t [s]', fontsize=18) P.legend() P.xlim(-0.1,0.05)

Out[36]:

(-0.1, 0.05)

Check matches¶

In [37]:

from pycbc.types import TimeSeries, FrequencySeries
from pycbc import filter as _filter
import pycbc.waveform.utils as wutils
import pycbc.types as pt
import pycbc.waveform as pw
import pycbc.psd as pp

from pycbc.types import TimeSeries, FrequencySeries from pycbc import filter as _filter import pycbc.waveform.utils as wutils import pycbc.types as pt import pycbc.waveform as pw import pycbc.psd as pp

In [38]:

def waveform_basis(hp, hc, psd=None, flow=None, ffinal=None):
    '''Given (hp, hc), the function gives orthonormal basis vectors.
    From appendix B (eqns B3, B4) of
    https://journals.aps.org/prd/pdf/10.1103/PhysRevD.95.024010
    '''
    hptilde = _filter.make_frequency_series(hp)
    hctilde = _filter.make_frequency_series(hc)
    sig1 = _filter.sigma(hp, psd, flow, ffinal)
    sig2 = _filter.sigma(hc, psd, flow, ffinal)

    proj = simple_inner(hptilde, hctilde, psd, flow, ffinal, sig1, sig2).real

    if isinstance(hp, pt.TimeSeries):
        hpp = pt.TimeSeries(hp.data/sig1, delta_t=hc.delta_t, epoch=hc._epoch)
        hper = pt.TimeSeries(hc.data/sig2, delta_t=hc.delta_t, epoch=hc._epoch)
    elif isinstance(hp, pt.FrequencySeries):
        hpp = pt.FrequencySeries(hp.data/sig1, delta_f=hc.delta_f, epoch=hc._epoch)
        hper = pt.FrequencySeries(hc.data/sig2, delta_f=hc.delta_f, epoch=hc._epoch)

    hper.data = (hper.data - proj*hpp.data)/np.sqrt(1-proj*proj)
    hper.data /= _filter.sigma(hper, psd, flow, ffinal)

    return hpp, hper

def simple_inner(htilde, stilde, psd=None, flow=None, fhigh=None, norm1=None, norm2=None):
    kmin, kmax = _filter.get_cutoff_indices(flow, fhigh, htilde.delta_f, (len(htilde)-1) * 2)
    indices = slice(kmin, kmax)
    if norm1 and norm2:
        norm = norm1*norm2
    else:
        norm1 = _filter.sigma(htilde, psd, flow, fhigh)
        norm2 = _filter.sigma(stilde, psd, flow, fhigh)
        norm = norm1*norm2

    if psd:
        return (np.conjugate(htilde.data[indices])*stilde.data[indices]/psd.data[indices]).sum()*4.0*psd.delta_f / norm
    else:
        return (np.conjugate(htilde.data[indices])*stilde.data[indices]).sum()*4.0*htilde.delta_f / norm

def minmax_match_with_basis(basis1, basis2, psd=None, flow=None, fhigh=None, norm1=None, norm2=None):
    '''Return phases min-max matches from orthonormal bases.
    From appendix B (eqns B10 - B14) of
    https://journals.aps.org/prd/pdf/10.1103/PhysRevD.95.024010
    '''
    assert len(basis1[0]) == len(basis2[0]), "Length of both basis do not match"
    match11 = _filter.matched_filter(basis1[0], basis2[0], psd, flow, fhigh, norm1)
    match12 = _filter.matched_filter(basis1[0], basis2[1], psd, flow, fhigh, norm2)
    match21 = _filter.matched_filter(basis1[1], basis2[0], psd, flow, fhigh, norm1)
    match22 = _filter.matched_filter(basis1[1], basis2[1], psd, flow, fhigh, norm2)

    a = match11.real().data*match11.real().data + match21.real().data*match21.real().data

    b = match12.real().data*match12.real().data + match22.real().data*match22.real().data

    c = match11.real().data*match12.real().data + match21.real().data*match22.real().data

    delta = np.sqrt((a - b)*(a - b) + 4*c*c)
    min_max = np.sqrt((a+b-delta)/2.0)
    max_max = np.sqrt((a+b+delta)/2.0)
    return min_max, max_max


def get_pycbc_timeseries(data, dt, duration, taper='start'):
    ts = wutils.taper_timeseries(TimeSeries(data, delta_t=dt), tapermethod=taper)
    ts.resize(int(duration/dt))
    return ts

def waveform_basis(hp, hc, psd=None, flow=None, ffinal=None): '''Given (hp, hc), the function gives orthonormal basis vectors. From appendix B (eqns B3, B4) of https://journals.aps.org/prd/pdf/10.1103/PhysRevD.95.024010 ''' hptilde = _filter.make_frequency_series(hp) hctilde = _filter.make_frequency_series(hc) sig1 = _filter.sigma(hp, psd, flow, ffinal) sig2 = _filter.sigma(hc, psd, flow, ffinal) proj = simple_inner(hptilde, hctilde, psd, flow, ffinal, sig1, sig2).real if isinstance(hp, pt.TimeSeries): hpp = pt.TimeSeries(hp.data/sig1, delta_t=hc.delta_t, epoch=hc._epoch) hper = pt.TimeSeries(hc.data/sig2, delta_t=hc.delta_t, epoch=hc._epoch) elif isinstance(hp, pt.FrequencySeries): hpp = pt.FrequencySeries(hp.data/sig1, delta_f=hc.delta_f, epoch=hc._epoch) hper = pt.FrequencySeries(hc.data/sig2, delta_f=hc.delta_f, epoch=hc._epoch) hper.data = (hper.data - proj*hpp.data)/np.sqrt(1-proj*proj) hper.data /= _filter.sigma(hper, psd, flow, ffinal) return hpp, hper def simple_inner(htilde, stilde, psd=None, flow=None, fhigh=None, norm1=None, norm2=None): kmin, kmax = _filter.get_cutoff_indices(flow, fhigh, htilde.delta_f, (len(htilde)-1) * 2) indices = slice(kmin, kmax) if norm1 and norm2: norm = norm1*norm2 else: norm1 = _filter.sigma(htilde, psd, flow, fhigh) norm2 = _filter.sigma(stilde, psd, flow, fhigh) norm = norm1*norm2 if psd: return (np.conjugate(htilde.data[indices])*stilde.data[indices]/psd.data[indices]).sum()*4.0*psd.delta_f / norm else: return (np.conjugate(htilde.data[indices])*stilde.data[indices]).sum()*4.0*htilde.delta_f / norm def minmax_match_with_basis(basis1, basis2, psd=None, flow=None, fhigh=None, norm1=None, norm2=None): '''Return phases min-max matches from orthonormal bases. From appendix B (eqns B10 - B14) of https://journals.aps.org/prd/pdf/10.1103/PhysRevD.95.024010 ''' assert len(basis1[0]) == len(basis2[0]), "Length of both basis do not match" match11 = _filter.matched_filter(basis1[0], basis2[0], psd, flow, fhigh, norm1) match12 = _filter.matched_filter(basis1[0], basis2[1], psd, flow, fhigh, norm2) match21 = _filter.matched_filter(basis1[1], basis2[0], psd, flow, fhigh, norm1) match22 = _filter.matched_filter(basis1[1], basis2[1], psd, flow, fhigh, norm2) a = match11.real().data*match11.real().data + match21.real().data*match21.real().data b = match12.real().data*match12.real().data + match22.real().data*match22.real().data c = match11.real().data*match12.real().data + match21.real().data*match22.real().data delta = np.sqrt((a - b)*(a - b) + 4*c*c) min_max = np.sqrt((a+b-delta)/2.0) max_max = np.sqrt((a+b+delta)/2.0) return min_max, max_max def get_pycbc_timeseries(data, dt, duration, taper='start'): ts = wutils.taper_timeseries(TimeSeries(data, delta_t=dt), tapermethod=taper) ts.resize(int(duration/dt)) return ts

In [39]:

## Fixing duration so as to get correct delta_f and PSD

duration = 16.0
df = 1./duration
flen = int(duration/dt) // 2 + 1

psd = pp.aLIGOZeroDetHighPowerGWINC(flen, df,  10)
psd.data[-1] = psd.data[-2]

## Fixing duration so as to get correct delta_f and PSD duration = 16.0 df = 1./duration flen = int(duration/dt) // 2 + 1 psd = pp.aLIGOZeroDetHighPowerGWINC(flen, df, 10) psd.data[-1] = psd.data[-2]

In [40]:

hp_eob = get_pycbc_timeseries(h_eob.real, dt, duration)
hc_eob = get_pycbc_timeseries(h_eob.imag, dt, duration)

hp_nr = get_pycbc_timeseries(h_nr.real, dt, duration)
hc_nr = get_pycbc_timeseries(h_nr.imag, dt, duration)

plt.figure(figsize=(10, 6))
plt.plot(hp_eob.sample_times, hp_eob, ls='dashed', label='EOB')
plt.plot(hp_nr.sample_times, hp_nr, label='NR')
plt.xlim(0, 2)

hp_eob = get_pycbc_timeseries(h_eob.real, dt, duration) hc_eob = get_pycbc_timeseries(h_eob.imag, dt, duration) hp_nr = get_pycbc_timeseries(h_nr.real, dt, duration) hc_nr = get_pycbc_timeseries(h_nr.imag, dt, duration) plt.figure(figsize=(10, 6)) plt.plot(hp_eob.sample_times, hp_eob, ls='dashed', label='EOB') plt.plot(hp_nr.sample_times, hp_nr, label='NR') plt.xlim(0, 2)

Out[40]:

(0.0, 2.0)

In [41]:

mat = _filter.match(hp_eob, hp_nr, psd, f_low)

basis = waveform_basis(hp_nr, hc_nr, psd, f_low)
basis1 = waveform_basis(hp_eob, hc_eob, psd, f_low)
min_max, max_max = minmax_match_with_basis(basis, basis1, psd, f_low)

print("matches with Adv LIGO PSD")
print(mat, min_max.max(), max_max.max())

mat = _filter.match(hp_eob, hp_nr, psd, f_low) basis = waveform_basis(hp_nr, hc_nr, psd, f_low) basis1 = waveform_basis(hp_eob, hc_eob, psd, f_low) min_max, max_max = minmax_match_with_basis(basis, basis1, psd, f_low) print("matches with Adv LIGO PSD") print(mat, min_max.max(), max_max.max())

matches with Adv LIGO PSD
(0.9862393750457854, 65522) 0.9245868078815853 0.9524529103830763

In [42]:

mat = _filter.match(hp_eob, hp_nr, None, f_low)

basis = waveform_basis(hp_nr, hc_nr, None, f_low)
basis1 = waveform_basis(hp_eob, hc_eob, None, f_low)
min_max, max_max = minmax_match_with_basis(basis, basis1, None, f_low)

print("matches with flat PSD or without PSD")
print(mat, min_max.max(), max_max.max())

mat = _filter.match(hp_eob, hp_nr, None, f_low) basis = waveform_basis(hp_nr, hc_nr, None, f_low) basis1 = waveform_basis(hp_eob, hc_eob, None, f_low) min_max, max_max = minmax_match_with_basis(basis, basis1, None, f_low) print("matches with flat PSD or without PSD") print(mat, min_max.max(), max_max.max())

matches with flat PSD or without PSD
(0.991265397844026, 65522) 0.9195209992794233 0.9371959363728289

In [ ]: