In [9]:
%matplotlib inline
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

fig_width_pt = 3*246.0  # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27               # Convert pt to inch
golden_mean = (np.sqrt(5)-1.0)/2.0         # Aesthetic ratio
fig_width = fig_width_pt*inches_per_pt  # width in inches
fig_height = fig_width*golden_mean      # height in inches
fig_size =  [fig_width,fig_height]

params = { 'axes.labelsize': 24,
          'font.family': 'serif',
          'font.serif': 'Computer Modern Raman',
          'font.size': 24,
          'legend.fontsize': 20,
          'xtick.labelsize': 24,
          'ytick.labelsize': 24,
          'axes.grid' : True,
          'text.usetex': True,
          'savefig.dpi' : 100,
          'lines.markersize' : 14,
          'figure.figsize': fig_size}

mpl.rcParams.update(params)
In [10]:
import lal
import lalsimulation as lalsim
import numpy as np
from pycbc.waveform import get_td_waveform

Calculate the phase offset and align the waveforms

In [11]:
from scipy import interpolate
def coa_phase_offset(anotherwf,**kwargs):
    tidx = []
    wfinterp = []
    for apx in [anotherwf,'SEOBNRE']:
        hp, hc = get_td_waveform(approximant=apx,
                                 **kwargs)
        wfinterp.append(interpolate.interp1d(hp.sample_times, np.unwrap(np.angle(np.vectorize(complex)(hp,-hc)))))
        tidx.append(hp.sample_times[0])
    offset = wfinterp[0](max(tidx)) - wfinterp[1](max(tidx))
    while offset < 0:
        offset = offset + 2*np.pi
    return offset/2
In [12]:
anotherwf = 'SEOBNRv1'

kwargs = {'mass1':10,
          'mass2':10,
          'spin1z':0.1,
          'spin2z':0.1,
          'distance':400,
          'delta_t':1./4096,
          'eccentricity':0,
          'coa_phase':100,
          'f_lower':20,
          'long_asc_nodes':700}

for apx in [anotherwf,'SEOBNRE']:
    if apx == 'SEOBNRE':
        ls = '-.'
    else:
        ls = '-'
    if apx == 'SEOBNRE':
        offset = coa_phase_offset(anotherwf,**kwargs)
        kwargs['coa_phase'] = kwargs['coa_phase'] - offset
    hp, hc = get_td_waveform(approximant=apx,
                                 **kwargs)
    plt.plot(hp.sample_times, hp, ls=ls,label='hplus:'+apx)
    #plt.plot(hc.sample_times, hc, ls=ls,label='hcross:'+apx)
plt.legend()
plt.xlim(-0.05,0.03)
plt.xlabel('Time [s]')
plt.ylabel('Strain')
Out[12]:
Text(0, 0.5, 'Strain')
In [13]:
anotherwf = 'SEOBNRv2'

for apx in [anotherwf,'SEOBNRE']:
    if apx == 'SEOBNRE':
        ls = '-.'
    else:
        ls = '-'
    if apx == 'SEOBNRE':
        offset = coa_phase_offset(anotherwf,**kwargs)
        kwargs['coa_phase'] = kwargs['coa_phase'] - offset
    hp, hc = get_td_waveform(approximant=apx,
                                 **kwargs)
    plt.plot(hp.sample_times, hp, ls=ls,label='hplus:'+apx)
    #plt.plot(hc.sample_times, hc, ls=ls,label='hcross:'+apx)
plt.legend()
plt.xlim(-0.05,0.03)
plt.xlabel('Time [s]')
plt.ylabel('Strain')
Out[13]:
Text(0, 0.5, 'Strain')
In [14]:
anotherwf = 'SEOBNRv3'

for apx in [anotherwf,'SEOBNRE']:
    if apx == 'SEOBNRE':
        ls = '-.'
    else:
        ls = '-'
    if apx == 'SEOBNRE':
        offset = coa_phase_offset(anotherwf,**kwargs)
        kwargs['coa_phase'] = kwargs['coa_phase'] - offset
    hp, hc = get_td_waveform(approximant=apx,
                                 **kwargs)
    plt.plot(hp.sample_times, hp, ls=ls,label='hplus:'+apx)
    #plt.plot(hc.sample_times, hc, ls=ls,label='hcross:'+apx)
plt.legend()
plt.xlim(-0.05,0.03)
plt.xlabel('Time [s]')
plt.ylabel('Strain')
Out[14]:
Text(0, 0.5, 'Strain')
In [ ]:
anotherwf = 'SEOBNRv4'

for apx in [anotherwf,'SEOBNRE']:
    if apx == 'SEOBNRE':
        ls = '-.'
    else:
        ls = '-'
    if apx == 'SEOBNRE':
        offset = coa_phase_offset(anotherwf,**kwargs)
        kwargs['coa_phase'] = kwargs['coa_phase'] - offset
    hp, hc = get_td_waveform(approximant=apx,
                                 **kwargs)
    plt.plot(hp.sample_times, hp, ls=ls,label='hplus:'+apx)
    #plt.plot(hc.sample_times, hc, ls=ls,label='hcross:'+apx)
plt.legend()
plt.xlim(-0.05,0.03)
plt.xlabel('Time [s]')
plt.ylabel('Strain')
In [ ]: