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/NANOGrav_9y/'

parpath = parentpath + 'par'
timpath = parentpath + 'tim'
noisepath = parentpath + 'noisefiles'


    

In [ ]:

    

# Find the parfiles and timfiles
parfiles = sorted(glob.glob(parpath+'/*.gls.par'))
timfiles = sorted(glob.glob(timpath+'/*.tim'))


    

In [ ]:

    

parfiles = [x for x in parfiles 
            if x!=parpath+'/J1713+0747_NANOGrav_9yv1.gls.par' ]


    

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

    

# 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() + '/9yr_psr_hdf5_files/DE421/tm_svdstabilized/'
if not os.path.exists(dirname):
    os.makedirs(dirname)


    

In [ ]:

    

# Dumpy 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_nano9yr.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 [ ]:

    

# Create ranked information text file
rankings = np.genfromtxt('./PsrListings_GWB.txt',dtype=str)[1:,0]
fil = open(dirname + '/psrList_nano9yr_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()