Pulsar_timing_model

#this is the code for plotting different time delays: solar system Roemer delay, Dispersion measure and Period derivative. p.s More model delays will be added in the future.
# the code generates one plot which shows individual delays and delay sum at the bottom.
# Solar system Roemer delay is the classical light travel time across earth orbit. It is the simple vaccum delay between the arrival of pulse at the obervatory and the solar system barycentric. To illustrate, imagine an Earth-Sun-Pulsar system, when earth is at the far orbit( of sun) relative to pulsar, the roemer delay will be the biggest and vice versa. The following plot starts at vernal equinox, where earth is the futherest from the pulsar, which creates the biggest delay. Later, as earth travels around the sun through autumn equinox and solstice points, the delay first gets smaller and then  bigger. (a continous cosine wave)
# Dispersion Measure delay is the delay due to dispersion in the interstella medium, so that the arrival time of pulse depends on the observing frequency. Here the oberserving frequency is set to be 600 Mhz and 800 Mhz. 
# Another DM delay changes as a function of time. It has a base value DM and a changing variable DMX. So, DM as a fucntion of time = DM + DMX
# Period Derivative is the pulsar spin down rate. As pulsar gets older and olders, it gradually slows down which will cause a time delay.
# the y-axis unit is in seconds

# HOW TO RUN:
# to RUN this code, you will need to type in pulsar's par file name in the command line. eg: python pulsar_timing_model.py J1714_3707Nanoxxxx.par

import sys
import matplotlib.pyplot as plt
import math
from pylab import *
from numpy import *


def Cel2Ecl( A, B ):
  betaa = math.asin( math.sin(math.radians(B) )*math.cos(math.radians(23.43) )-math.sin(math.radians(23.43) )*math.sin(math.radians(A))*math.cos(math.radians(B) ))
  lamdaa = math.acos( math.cos(math.radians(A) )*math.cos(math.radians(B))/math.cos(math.radians(betaa)))
  return (lamdaa,betaa)

def remove_space( A ):
    ret = []
    for curr in A:
        if curr != '':
            ret.append( curr )
    return ret





inFile = open( sys.argv[1], "r" )

All = inFile.readlines()

# From here, this code automatically read all the parameters needed for plotting.It should be working for every single par file. 

DMR = []

for line in All:
    X = line.split(' ')
    if X[0][0:4] == "DMX_":
        X = remove_space( X )
        DMR.append( float( X[1] ) )
DMX = DMR + DMR

DMXR1 = []

for line in All:
    X = line.split(' ')
    if X[0][0:5] == "DMXR1":
        X = remove_space( X )
        DMXR1.append( float( X[1] ) )

DMXR2 = []

for line in All:
    X = line.split(' ')
    if X[0][0:5] == "DMXR2":
        X = remove_space( X )
        DMXR2.append( float( X[1] ) )

#DM, FO, F1, PEPOCH, DECJ, RAJ, M2, START, FINISH. 

for line in All:
    X = line.split(' ')
    if X[0][0:4] == "DM":
        X = remove_space( X )
        DM = float( X[1] )



for line in All:
    X = line.split(' ')
    if X[0][0:4] == "F0":
        X = remove_space( X )
        F0=  float( X[1] )


for line in All:
    X = line.split(' ')
    if X[0][0:4] == "F1":
        X = remove_space( X )
        F1 = float( X[1] )


for line in All:
    X = line.split(' ')
    if X[0][0:4] == "DECJ":
        X = remove_space( X )
        DECJ =  X[1].split(':')


for line in All:
    X = line.split(' ')
    if X[0][0:3] == "RAJ":
        X = remove_space( X )
        RAJ = X[1].split(':')


for line in All:
    X = line.split(' ')
    if X[0][0:4] == "M2":
        X = remove_space( X )
        M2 = float( X[1] )


for line in All:
    X = line.split(' ')
    if X[0][0:6] == "PEPOCH":
        X = remove_space( X )
        PEPOCH = float( X[1] )


for line in All:
    X = line.split(' ')
    if X[0][0:5] == "START":
        X = remove_space( X )
        START = float( X[1] )


for line in All:
    X = line.split(' ')
    if X[0][0:6] == "FINISH":
        X = remove_space( X )
        FINISH = float( X[1] )


#sort DMX and DM
for i in range(len(DMX)):
    DMX[i]+=DM

DMXRI = DMXR1 + DMXR2
index = {}
for i in range(len(DMXRI)):
    index[DMXRI[i]]=i

DMXRI.sort()
d=[]
for i in range(len(DMXRI)):
    d.append(DMX[index[DMXRI[i]]])
    
#convert DM to delays for plotting 
b = [x * 4.419 * 10**5 * (600**-2 - 800**-2) * 0.001 for x in d]


#convert RAJ and DECJ to degrees

rah = float( RAJ[0] )
ram = float( RAJ[1] )
ras = float( RAJ[2] )
deg = (rah) + (ram/60) + (ras/3600)

dech = float( DECJ[0] )
decm = float( DECJ[1] )
decs = float( DECJ[2] )
deg2 = (dech) + (decm/60) + (decs/3600)

# Convert RA and DEC to longitude and latitude
lamdaa, betaa = Cel2Ecl( deg, deg2 )




# range depends on whether the par file is a 5 year or 9 year data.
t = arange(START, FINISH, 0.01)

s1 = 500  * cos(betaa) * cos(( t  - START  ) * 2 * pi  /365.25 + lamdaa) #here roemer has a much larger coefficient compare to the other two, so I reduce its value by 10^-22 to make them on the same scale
s2 = - F1 * ((t - START) * 2 * pi /365.25)  /( F0**2)
s3 = 4.149 * 10**5 * DM * t/t * (600**-2 - 800**-2) * 0.001 

ax4 = subplot(414)
plt.plot(t,s1 + s2 - s3,'c',label='delay sum')
setp( ax4.get_xticklabels(), fontsize=10)
plt.legend(bbox_to_anchor=(0, 1), loc=2, borderaxespad=0.,prop={'size':7})


ax1 = subplot(411,sharex=ax4)
plot(t,s1,'r',label='SS \n Roemer \n Delay')
plt.title('Pulsar Timing Models \n (x-axis:date y-axis:seconds ) ')
setp( ax1.get_xticklabels(), visible=False)
plt.legend(bbox_to_anchor=(0, 1), loc=2, borderaxespad=0.,prop={'size':7})


ax2 = subplot(412, sharex=ax4)
plt.plot(t,s2,'b',label='Period Derivatives')
setp( ax2.get_xticklabels(), visible=False)
plt.legend(bbox_to_anchor=(0, 1), loc=2, borderaxespad=0.,prop={'size':10})


ax3 = subplot(413, sharex=ax4)
plt.plot(DMXRI, b, 'k_-',label='DMxth')
plt.plot(DMXRI, b, 'g-')
setp( ax3.get_xticklabels(), visible=False)
plt.legend(bbox_to_anchor=(0, 1), loc=2, borderaxespad=0.,prop={'size':8.5})
ax4.yaxis.label.set_color('g')


plt.show()