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 [ ]: