In [1]:
import NRSur7dq2

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

Load the NRSur7dq2 surrogate

In [2]:
sur = NRSur7dq2.NRSurrogate7dq2('NRSur7dq2.h5')

Read the surrogate documentation

In [3]:
sur?

Generate a waveform

In [4]:
q = 1.7 # Mass ratio. Should lie in the interval [1, 2].
chiA0 = np.array([0.6, 0.3, 0.1]) # Dimensionless spin on the larger black hole
chiB0 = np.array([0.5, -0.2, -0.2]) # Dimensionless spin on the smaller black hole

modes = sur(q, chiA0, chiB0)

# Plot the (2, 2) mode
ell, m = 2, 2
plt.plot(sur.t_coorb, modes[ell, m].real, label='real part of (%s, %s) mode'%(ell, m))
plt.plot(sur.t_coorb, modes[ell, m].imag, label='imaginary part of (%s, %s) mode'%(ell, m))
plt.legend(frameon=False, loc='upper left')
plt.ylabel("$rh/M$")
plt.xlabel("$t/M$")
plt.show()

# Plot the (2, 1) mode
ell, m = 2, 1
plt.plot(sur.t_coorb, modes[ell, m].real, label='real part of (%s, %s) mode'%(ell, m))
plt.plot(sur.t_coorb, modes[ell, m].imag, label='imaginary part of (%s, %s) mode'%(ell, m))
plt.legend(frameon=False, loc='upper left')
plt.ylabel("$rh/M$")
plt.xlabel("$t/M$")
plt.show()

print "Available modes:"
print sorted(modes.keys())
Available modes:
[(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)]

Evaluate an edge-on waveform at 500 megaparsecs for a system with a total mass of 70 solar masses in the detector frame:

In [5]:
# Evaluate the waveform on a uniformly-spaced array of times
dt = 1. / 4096. # in seconds
sample_times = np.arange(-0.3, 0.03, dt) # The peak amplitude occurs roughly at t=0

h = sur(q, chiA0, chiB0, theta=np.pi/2, phi=0., MTot=70.0, distance=500, t=sample_times)

h_plus = np.real(h)
h_cross = -1*np.imag(h)

plt.plot(sample_times, h_plus, label='$h_+$')
plt.plot(sample_times, h_cross, label='$h_x$')
plt.xlabel("time [seconds]")
plt.ylabel("strain")
plt.legend(frameon=False, loc='upper left')
plt.show()

The spin directions can be specified at a reference time or frequency, and time-dependent spins can be obtained:

In [6]:
modes, chiA, chiB = sur(q, chiA0, chiB0, t_ref=-1000, return_spins=True)

# The amplitude of the (2, 1) mode is small at t_ref, where the inertial frame is aligned
plt.plot(sur.t_coorb, abs(modes[2, 2]), label='Amplitude of (2, 2) mode')
plt.plot(sur.t_coorb, abs(modes[2, 1]), label='Amplitude of (2, 1) mode')
plt.legend(frameon=False, loc='upper left')
plt.ylabel("$rh/M$")
plt.xlabel("$t/M$")
plt.show()

# The spins agree with the provided spins at t_ref:
plt.plot(sur.t_coorb, chiA)
plt.plot([-1000], np.array([chiA0]), 'ko')
plt.ylabel("$\\vec{\chi}_A$")
plt.xlabel("$t/M$")
plt.show()

# For f_ref to be given in Hz, the total mass must be specified:
modes, chiA, chiB = sur(q, chiA0, chiB0, f_ref=20, MTot=70, return_spins=True)

# The amplitude of the (2, 1) mode is small at f_ref, where the inertial frame is aligned
plt.plot(sur.t_coorb, abs(modes[2, 2]), label='Amplitude of (2, 2) mode')
plt.plot(sur.t_coorb, abs(modes[2, 1]), label='Amplitude of (2, 1) mode')
plt.legend(frameon=False, loc='upper left')
plt.ylabel("$rh/M$")
plt.xlabel("$t/M$")
plt.show()
In [ ]: