In [1]:
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
%load_ext Cython
%load_ext autoreload
%autoreload 2
from __future__ import division
import matplotlib
import matplotlib.pyplot as plt
matplotlib.rcParams['savefig.dpi'] = 1.5 * matplotlib.rcParams['savefig.dpi']
import numpy as np
import sys,os,glob,h5py
import libstempo as T2
import libstempo.plot as LP
import NX01_psr
import NX01_datafile
try:
from IPython.core.display import clear_output
have_ipython = True
except ImportError:
have_ipython = False
Msol = 1.98855*10.0**30.0
In [ ]:
parentpath = '/Users/staylor/Research/NANOGrav/stochastic_11yr_analysis/data/'
parpath = parentpath + 'partim_no_noise'
timpath = parentpath + 'partim_no_noise'
noisepath = parentpath + 'nano_11_noisefiles_30_marg'
In [ ]:
# Find the parfiles and timfiles
parfiles = sorted(glob.glob(parpath+'/*.par'))
timfiles = sorted(glob.glob(timpath+'/*.tim'))
In [ ]:
# Find the noise files
noisefiles = sorted(glob.glob(noisepath+'/*.txt'))
In [ ]:
len(parfiles), len(timfiles), len(noisefiles)
In [ ]:
stripped_pars = list(parfiles)
In [ ]:
for ii in range(len(stripped_pars)):
stripped_pars[ii] = stripped_pars[ii].replace('9yv1.gls.par', '9yv1.gls.strip.par')
stripped_pars[ii] = stripped_pars[ii].replace('9yv1.t2.gls.par', '9yv1.t2.gls.strip.par')
In [ ]:
for ii in range(len(stripped_pars)):
os.system('awk \'($1 !~ /T2EFAC/ && $1 !~ /T2EQUAD/ && $1 !~ /ECORR/ && $1 !~ /RNAMP/ && $1 !~ /RNIDX/ ) {{print $0}} \' {0} > {1}'.format(parfiles[ii],stripped_pars[ii]))
In [ ]:
parfiles = sorted(glob.glob(parpath+'/*.strip.par'))
In [ ]:
#######################################
# PASSING THROUGH TEMPO2 VIA libstempo
#######################################
t2psr = []
for ii in range(len(parfiles)):
t2psr.append( T2.tempopulsar(parfile = parfiles[ii], timfile = timfiles[ii],
maxobs=30000, ephem='DE421') )
#if np.any(np.isfinite(t2psr[ii].residuals())==False)==True:
# t2psr[ii] = T2.tempopulsar(parfile = parfiles[ii], timfile = timfiles[ii])
clear_output()
print '\r', '{0} of {1}'.format(ii+1,len(parfiles))
sys.stdout.flush()
In [ ]:
# Check out some plots if you want
#LP.plotres(t2psr[0])
plt.errorbar(t2psr[0].toas(),t2psr[0].residuals()/1e-6,t2psr[0].toaerrs,fmt='.',alpha=0.2)
plt.xlabel(r'MJD')
plt.ylabel(r'Residuals [$\mu$s]')
plt.show()
In [ ]:
## Cell for time-slicing of data
Nyears = 'tot'
Tmin = np.min([np.min(p.toas()) for p in t2psr])
if Nyears == 'tot':
Tcutoff = np.max([np.max(p.toas()) for p in t2psr])
else:
Tcutoff = Tmin + Nyears*365.25
ind_slice = []
for ii,p in enumerate(t2psr):
mask = np.logical_and(p.toas() >= Tmin, p.toas() <= Tcutoff)
if np.sum(mask) > 0:
Tobs = (p.toas()[mask].max() - p.toas()[mask].min()) / 365.25
if Tobs > 1.0:
ind_slice.append(ii)
In [ ]:
# Pass all tempopulsar objects to NX01 pulsar class
psr = [NX01_psr.PsrObj(t2psr[ii]) for ii in ind_slice]
[p.grab_all_vars(jitterbin=1.0,makeGmat=False,
fastDesign=False,planetssb=True,
startMJD=Tmin, endMJD=Tcutoff) for p in psr] # according to the 9 year paper,
# the jitterbin used is 1s
In [ ]:
# Pass all tempopulsar objects to NX01 pulsar class
#psr = [NX01_psr.PsrObj(p) for p in t2psr]
#[p.grab_all_vars(jitterbin=1.0,makeGmat=False,
# fastDesign=False,planetssb=True,
# startMJD=Tmin, endMJD=Tcutoff) for p in psr] # according to the 9 year paper,
# # the jitterbin used is 1s
In [ ]:
# Fill in the locations of par, tim, and noise files
for ii,jj in enumerate(ind_slice):
psr[ii].parfile = parfiles[jj]
psr[ii].timfile = timfiles[jj]
psr[ii].noisefile = noisefiles[jj]
In [ ]:
# Only need to execute if you want roemer delays
# from other ephemerides.
# Dummy libstempo passes to grab delays.
ephems = ['DE421', 'DE430', 'DE435', 'DE436']
for ii,jj in enumerate(ind_slice):
for eph in ephems:
if eph != psr[ii].ephemname:
dummy_t2psr = T2.tempopulsar(parfile = parfiles[jj], timfile = timfiles[jj],
maxobs=30000, ephem=eph)
psr[ii].roemer[eph] = np.double(dummy_t2psr.roemer)
if psr[ii].tmask is not None:
psr[ii].roemer[eph] = psr[ii].roemer[eph][psr[ii].tmask]
if psr[ii].isort is not None:
psr[ii].roemer[eph] = psr[ii].roemer[eph][psr[ii].isort]
clear_output()
print '\r', '{0} of {1}'.format(ii+1,len(parfiles))
sys.stdout.flush()
In [ ]:
# Make HDF5 file directory
dirname = os.getcwd() + '/11yr_psr_hdf5_files/DE421/tm_svdstabilized/'
if not os.path.exists(dirname):
os.makedirs(dirname)
In [ ]:
# Dump all pulsars into HDF5 files
for ii,p in enumerate(psr):
df = NX01_datafile.DataFile(dirname + '/' + p.name + '.hdf5')
df.addTempoPulsar(p)
clear_output()
print '\r', '{0} of {1}'.format(ii+1,len(psr))
sys.stdout.flush()
In [ ]:
# Create information text file for pulsar hdf5 file locations
fil = open(dirname + '/psrList_nano11yr.txt','w')
print >>fil, "NAME", "HDF5-PATH", "PARFILE-PATH", "TIMFILE-PATH"
print >>fil, "#############################################"
for p in psr:
print >>fil, p.name, dirname+'/'+p.name+'.hdf5'
fil.close()
In [ ]:
rankings = np.genfromtxt(parentpath+'psrlist.txt',dtype=str)
fil = open(dirname + '/psrList_nano11_ranked.txt','w')
print >>fil, "NAME", "HDF5-PATH", "PARFILE-PATH", "TIMFILE-PATH"
print >>fil, "#############################################"
for p in rankings:
print >>fil, p, dirname+'/'+p+'.hdf5'
fil.close()
In [ ]:
rankings = np.genfromtxt(parentpath+'psrlist_Tg3yr.txt',dtype=str)
fil = open(dirname + '/psrList_nano11_ranked_Tg3yr.txt','w')
print >>fil, "NAME", "HDF5-PATH", "PARFILE-PATH", "TIMFILE-PATH"
print >>fil, "#############################################"
for p in rankings:
print >>fil, p, dirname+'/'+p+'.hdf5'
fil.close()
In [ ]:
gmu_uls = np.array([3.87595541e-13, 4.71210227e-13, 5.71803686e-13,
6.95797061e-13, 8.32522097e-13, 1.02163421e-12,
1.26938114e-12, 1.54532333e-12, 1.86832268e-12,
2.25810134e-12, 2.83286775e-12, 3.51803788e-12,
4.52328994e-12, 6.02546201e-12, 7.67416070e-12,
1.04064636e-11, 1.39940969e-11, 2.06332855e-11,
3.08028537e-11, 5.48868339e-11])
In [ ]:
plt.loglog(gmu_uls,10.0**prob)
plt.fill_betweenx(10.0**prob,gmu_uls,1e-9,alpha=0.4)
plt.xlim(1e-13,1e-9)
plt.ylim(1e-3,1.0)
plt.xlabel(r'$G\mu / c^2$')
plt.ylabel(r'$p$')
plt.show()