GalDynPsr 1.0.0

A practise package

GalDynPsr Package

GalDynPsr is a package for calculating dynamical contributions in the period derivatives of pulsars in the Galactic field. It can calculate the fractional contributions or the excess terms, e.g. \dot{P}/P|_excess where P is either the orbital period or the spin period. Various dynamical contributions, including the Shklovskii effect (the contribution due to the proper motion) or due to the acceleration of the pulsar caused by the gravitational potential of the Galaxy, can be calculated. Various methods or models to perform these tasks are available in the package. Using the measured values of the periods and period derivatives, GalDynPsr can even compute the "intrinsic" values of the period derivatives, provided no other extra contribution exist.

Details on various dynamical effects and various models to estimate those are available in Pathak and Bagchi (arXiv:1712.06590).

Brief outline of usage of GalDynPsr is given below.

Questions and queries can be sent to Dhruv Pathak pathakdhruv9786 at gmail dot com

1) After installing numpy, scipy, and galpy successfully, install GalDynPsr

pip3 install GalDynPsr

If wished, one can change the values of Rs (galactocentric cylindrical radius of the sun), Vs (rotational speed of the sun) etc in parameters.in file that can be found inside the GalDynPsr (installed directory).

But remember that galpy also has these values defined in the file '$home/.galpyrc'. One can in principle change these values.

However, the Milky Way potential in galpy was fitted with Rs = 8.0 kpc and vs = 220.0 in galpy

2) Import GalDynPsr in your code as

import GalDynPsr

3) In all of the following modules, the ordering of the arguments of the functions are very important.

4) Calculate excess terms of modelX ** Very important, model names are case sensitive!

GalDynPsr.modelX.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) and GalDynPsr.modelX.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd)

The function Expl() calculates the excess term due to the relative acceleration of the pulsar parallel to the Galactic plane and the function Exz() calculates the excess term due to the relative acceleration of the pulsar perpendicular to the Galactic plane. Total dynamical contribution will be the sum of above two terms. Both of the above functions takes the arguments as the Galactic longitude in degrees (say ldeg), uncertainty in the Galactic longitude in degrees (say sigl), the Galactic latitude in degrees (say bdeg), uncertainty in the Galactic latitude in degrees (say sigb), the distance of the pulsar from the solar system barycenter in kpc (say dkpc), and the uncertainty in the distance in kpc (say sigd); strictly in this order. One needs to assign the values of ldeg, sigl, bdeg, sigb, dkpc, and sigd before calling the above functions.

The models that do not involve galpy can return uncertainities in the excess terms following standard error propagation method, like

modelX.Errpl(ldeg, sigl, bdeg, sigb, dkpc, sigd) and modelX.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd)

where the function Errpl() returns the uncertainty in the excess terms due to the relative acceleration of the pulsar parallel to the Galactic plane and Errz() returns the uncertainty in the excess terms due to the relative acceleration of the pulsar perpendicular to the Galactic plane. The arguments are the same as before, and needs to be assigend values before using the functions.

Remeber to replace modelX by the model you want to use, like modelA, modelLa, etc. Model names are case sensitive!

5) Calculate the Shklovskii term

The Shklovskii term can be calculated as GalDynPsr.Shk.Exshk(dkpc, sigd, mu_alpha, sigmua, mu_delta, sigmud) and the error in the Shklovskii term as GalDynPsr.Shk.errShk(dkpc, sigd, mu_alpha, sigmua, mu_delta, sigmud)

where mu_alpha is the proper motion in the Right Ascesion direction (mas/yr), mu_delta is the proper motion in the Declination direction (mas/yr), sigmua is the uncertainty in mu_alpha, and sigmud is the uncertainty in mu_delta. dkpc is as usual the distance of the pulsar from the solar system barycenter in kpc and sigd is the uncertainty in dkpc. One needs to assign the values of these parameters first.

6) Print the basic parameters of the pulsars

GalDynPsr.read_parameters.Rskpc returns the Galactocentric cylindrical radius of the sun.

GalDynPsr.read_parameters.Vs returns the rotational speed of the sun.

GalDynPsr.read_parameters.Rpkpc(ldeg, sigl, bdeg, sigb, dkpc, sigd) returns the value of Galactocentric cylindrical radius of the pulsar in kpc and GalDynPsr.read_parameters.z(ldeg, sigl, bdeg, sigb, dkpc, sigd) returns the verical height of the pulsar from the Galactic plane. The meaning of the arguments are as usual:

One can get the uncertainty in the Galactocentric cylindrical radius of the pulsar in kpc using GalDynPsr.read_parameters.ErrRpkpc(ldeg, sigl, bdeg, sigb, dkpc, sigd) and the uncertainty in the verical height as GalDynPsr.read_parameters.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd).

Obviously, one first needs to assign the values of ldeg, sigl, bdeg, sigb, dkpc, and sigd.

7) For a pulsar in a globular cluster, one can read the cluster parameters, e.g. the Galactic longitude, latitude and the distance from the Harric catalogue (GC_harris.txt coming with GalDynPsr), and calculate crude values of the dynamical terms

One first needs to import the module to read Harris cataloge as:

import GalDynPsr.ReadGC as readgc

Then one can input the name of the cluster, e.g. readgc.gcinput("47Tuc")

Then the excess terms using modelX can be calculated as

sigl=0.000001 sigb=0.000001 sigd=0.000001

GalDynPsr.modelX.Expl(readgc.ldeg,sigl,readgc.bdeg,sigb,readgc.dkpc,sigd) # the parallel component GalDynPsr.modelX.Exz(readgc.ldeg,sigl,readgc.bdeg,sigb,readgc.dkpc,sig) # the perpendicular component

Note that the Harris catalogue does not provide uncertainties in ldeg, bdeg and dkpc. That is why we need to define them before calculating the excess terms. As the pulsar parameters differ from the cluster parameters as in Harris catalogue, it is meaningless to calculate uncetrainities. So, one SHOULD simply call the above functions as:

GalDynPsr.modelX.Expl(readgc.ldeg,0.0,readgc.bdeg,0.0,readgc.dkpc,0.0) # the parallel component GalDynPsr.modelX.Exz(readgc.ldeg,0.0,readgc.bdeg,0.0,readgc.dkpc,0.0) # the perpendicular component

8) models involving galpy does not return error. However, that can be done using simulation. See Example-2 below.

9) Calculation of instrinsic period. Ordering of the arguments is important.

First assign values of ldeg, sigl, bdeg, sigb, dkpc, sigd (meaning of the parameters are as usual) and calculate the excess terms as usual:

Ex_pl = GalDynPsr.modelA.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) # excess term parallel to the Galactic plane

Ex_z = GalDynPsr.modelA.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) # excess term perpendicular to the Galactic plane

errpl = GalDynPsr.modelA.Errpl(ldeg, sigl, bdeg, sigb, dkpc, sigd) # uncertainty in the parallel term

errz = GalDynPsr.modelA.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd) # uncertainty in the parallel term

Ex_shk = GalDynPsr.Shk.Exshk(dkpc, sigd, mu_alpha, sigmua, mu_delta, sigmud) # excess term due to the Shklovskii effect

errshk = GalDynPsr.Shk.errShk(dkpc, sigd, mu_alpha, sigmua, mu_delta, sigmud) # uncertainties in the Shklovskii term

Now assign the values of the period 'P' in seconds, the uncertainy in period 'sigP' in seconds, the measured value of the period derivative 'Pdot_obs' in seconds/seconds and the uncertainty in the measured value of the period derivative 'sigPdot' in seconds/seconds; and calculate the dynamically caused values of the period derivatives (in seconds/seconds) as:

GalDynPsr.pdotint.PdotExpl(Ex_pl,P) # due to the parallel component of the acceleration

GalDynPsr.pdotint.PdotExz(Ex_z,P) # due to the perpendicular component of the acceleration

GalDynPsr.pdotint.PdotShk(Ex_shk,P) # due to the Shklovskii effect

The total dynamically caused value of Pdot is the addition of the above three terms. The intrinsic value of the period derivative can be calculated by subtracting that sum from the measured value of the period derivative. GalDynPsr can perform this task for us by:

GalDynPsr.pdotint.Pdotint(Ex_pl,Ex_z,Ex_shk,Pdot_obs,P)

The total dynamical contribution in the period derivative, i.e. the sum of GalDynPsr.pdotint.PdotExpl(Ex_pl,P) and GalDynPsr.pdotint.PdotExz(Ex_z,P) can be printed as: GalDynPsr.pdotint.PdotGal(Ex_pl,Ex_z,P)

One can calculate the uncertainity in the above period drivatives as:

GalDynPsr.pdotint.ErrPdotExpl(Ex_pl,errpl,P,sigP) # uncertainty in the planar contribution to the period derivative

GalDynPsr.pdotint.ErrPdotExz(Ex_z,errz,P,sigP) # uncertainty in the perpendicular contribution to the period derivative

GalDynPsr.pdotint.ErrPdotGal(Ex_pl,errpl,Ex_z,errz,P,sigP) # uncertainty in the sum of planar and perpendicular contributions to the period derivative

GalDynPsr.pdotint.ErrPdotShk(Ex_shk,errshk,P,sigP) # uncertainty in the Shklovskii term contribution to the period derivative

GalDynPsr.pdotint.ErrPdotint(Ex_pl,errpl,Ex_z,errz,Ex_shk,errshk,Pdot_obs,sigPdot_obs,P,sigP) # uncertainty in the intrinsic period derivative

• Remember the fact that GalDynPsr reads the value of the period in seconds. Orbital periods are usually measured in days (or sometimes in hours). One needs to convert those values in seconds before calling GalDynPsr modules.

All the above points will be clearer from the following examples

Instructions to use GalDynPsr

Call GalDynPsr as:-

import GalDynPsr

Provide the following values in your code

ldeg = Galactic logitude in degrees, bdeg = Galactic latitude in degrees, dkpc = distance to pulsar in kpc, sigl = error in ldeg, sigb = error in bdeg, sigd = error in dkpc

########### Example-1 ################# ldeg =20.0 sigl = 2.0

bdeg = 45.0 sigb = 0.2

dkpc = 1.2 sigd = 0.3

############# Extract important parameters say values of Rp (in kpc) and z (in kpc) and errors in them ##########

Rpkpc = GalDynPsr.read_parameters.Rpkpc(ldeg, sigl, bdeg, sigb, dkpc, sigd)

zkpc = GalDynPsr.read_parameters.z(ldeg, sigl, bdeg, sigb, dkpc, sigd)

ErrRp = GalDynPsr.read_parameters.ErrRp(ldeg, sigl, bdeg, sigb, dkpc, sigd)

Errz = GalDynPsr.read_parameters.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd)

################# Compute excess Shklovskii term. Here Exshk() calculates the dynamical contribution (or the excess term) due to Shklovskii effect and errShk() calculates the error in that #################################

We need to provide the values of proper motion in right ascesion direction (with the corresponding error) and proper motion in the declination direction (with the corresponding error)########

##mu_alpha = proper motion in Right Ascesion direction (mas/yr), mu_delta = proper motion in Declination direction (mas/yr), sigmua = error in mu_alpha, and sigmud = error in mu_delta

mu_alpha = 22.0 sigmua = 0.2

mu_delta = 1.2 sigmud = 0.3

ExcessSh = GalDynPsr.Shk.Exshk(dkpc, sigd, mu_alpha, sigmua, mu_delta, sigmud)

sigmaExcessSh = GalDynPsr.Shk.errShk(dkpc, sigd, mu_alpha, sigmua, mu_delta, sigmud)

Example-1.1

Model A (can propagate error)######

Apl = GalDynPsr.modelA.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

Az = GalDynPsr.modelA.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

Aplsigma = GalDynPsr.modelA.Errpl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the planar contribution to the excess term

Azsigma = GalDynPsr.modelA.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the perpendicular contribution to the excess term

totalA = Apl + Az

SigmatotalA = math.sqrt(AplsigmaAplsigma+AzsigmaAzsigma) # assuming no correlation between excepp_pl and excess_z

Example-1.2

Model B (can propagate error)######

Bpl = GalDynPsr.modelB.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

Bz = GalDynPsr.modelB.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

Bplsigma = GalDynPsr.modelB.Errpl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the planar contribution to the excess term

Bzsigma = GalDynPsr.modelB.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the perpendicular contribution to the excess term

totalB = Bpl + Bz

SigmatotalB = math.sqrt(BplsigmaBplsigma+BzsigmaBzsigma) # assuming no correlation between excepp_pl and excess_z

Example-1.3

Model C (can propagate error)######

Cpl = GalDynPsr.modelC.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

Cz = GalDynPsr.modelC.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

Cplsigma = GalDynPsr.modelC.Errpl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the planar contribution to the excess term

Czsigma = GalDynPsr.modelC.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the perpendicular contribution to the excess term

totalC = Cpl + Cz

SigmatotalC = math.sqrt(CplsigmaCplsigma+CzsigmaCzsigma) # assuming no correlation between excepp_pl and excess_z

Example-1.4

Model D (can propagate error)######

Dpl = GalDynPsr.modelD.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

Dz = GalDynPsr.modelD.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

Dplsigma = GalDynPsr.modelD.Errpl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the planar contribution to the excess term

Dzsigma = GalDynPsr.modelD.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the perpendicular contribution to the excess term

totalD = Dpl + Dz

SigmatotalD = math.sqrt(DplsigmaDplsigma+DzsigmaDzsigma) # assuming no correlation between excepp_pl and excess_z

Example-1.5

############ Model Ea (error only in parallel component) ###############

Eapl = GalDynPsr.modelEa.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

SigmaEapl = GalDynPsr.modelEa.Errpl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the planar contribution to the excess term

Eaz = GalDynPsr.modelEa.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

totalEa = Eapl+Eaz

Example-1.6

############ Model Eb (error only in parallel component) ###############

Ebpl = GalDynPsr.modelEb.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

SigmaEbpl = GalDynPsr.modelEb.Errpl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the planar contribution to the excess term

Ebz = GalDynPsr.modelEb.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

totalEb = Ebpl+Ebz

Example-1.7

############ Model Fa (error only in parallel component) ###############

Fapl = GalDynPsr.modelFa.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

SigmaFapl = GalDynPsr.modelFa.Errpl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the planar contribution to the excess term

Faz = GalDynPsr.modelFa.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

totalFa = Fapl+Faz

Example-1.8

############ Model Fb (error only in parallel component) ###############

Fbpl = GalDynPsr.modelFb.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

SigmaFbpl = GalDynPsr.modelFb.Errpl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the planar contribution to the excess term

Fbz = GalDynPsr.modelFb.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

totalFb = Fbpl+Fbz

Example-1.9

############ Model Ga (error in perpendicular component) ###############

Gapl = GalDynPsr.modelGa.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

Gaz = GalDynPsr.modelGa.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

SigmaGaz = GalDynPsr.modelGa.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the perpendicular contribution to the excess term

totalGa = Gapl+Gaz

Example-1.10

############ Model Gb (error in perpendicular component) ###############

Gbpl = GalDynPsr.modelGb.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

Gbz = GalDynPsr.modelGb.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

SigmaGbz = GalDynPsr.modelGb.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the perpendicular contribution to the excess term

totalGb = Gbpl+Gbz

Example-1.11

############ Model Ha (error in perpendicular component) ###############

Hapl = GalDynPsr.modelHa.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

Haz = GalDynPsr.modelHa.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

SigmaHaz = GalDynPsr.modelHa.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the perpendicular contribution to the excess term

totalHa = Hapl+Haz

Example-1.12

############ Model Hb (error in perpendicular component) ###############

Hbpl = GalDynPsr.modelHb.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

Hbz = GalDynPsr.modelHb.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

SigmaHbz = GalDynPsr.modelHb.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the perpendicular contribution to the excess term

totalHb = Hbpl+Hbz

Example-1.13

############ Model Ia (error in perpendicular component) ###############

Iapl = GalDynPsr.modelIa.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

Iaz = GalDynPsr.modelIa.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

SigmaIaz = GalDynPsr.modelIa.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the perpendicular contribution to the excess term

totalIa = Iapl+Iaz

Example-1.14

############ Model Ib (error in perpendicular component) ###############

Ibpl = GalDynPsr.modelIb.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

Ibz = GalDynPsr.modelIb.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

SigmaIbz = GalDynPsr.modelIb.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the perpendicular contribution to the excess term

totalIb = Ibpl+Ibz

Example-1.15

############ Model Ja (error in perpendicular component) ###############

Japl = GalDynPsr.modelJa.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

Jaz = GalDynPsr.modelJa.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

SigmaJaz = GalDynPsr.modelJa.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the perpendicular contribution to the excess term

totalJa = Japl+Jaz

Example-1.16

############ Model Jb (error in perpendicular component) ###############

Jbpl = GalDynPsr.modelJb.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

Jbz = GalDynPsr.modelJb.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

SigmaJbz = GalDynPsr.modelJb.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the error in the perpendicular contribution to the excess term

totalJb = Jbpl+Jbz

Example-1.17

############ Model Ka (full galpy, no error) ###############

Kapl = GalDynPsr.modelLa.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

Kaz = GalDynPsr.modelLa.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

totalKa = Kapl+Kaz

Example-1.18

############ Model Kb (full galpy, no error) ###############

Kbpl = GalDynPsr.modelLa.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

Kbz = GalDynPsr.modelLa.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

totalKb = Kbpl+Kbz

Example-1.19

############ Model La (full galpy, no error) ###############

Lapl = GalDynPsr.modelLa.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

Laz = GalDynPsr.modelLa.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

totalLa = Lapl+Laz

Example-1.20

############ Model Lb (full galpy, no error) ###############

Lbpl = GalDynPsr.modelLa.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the planar contribution to the excess term

Lbz = GalDynPsr.modelLa.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd) #calculates the perpendicular contribution to the excess term

totalLb = Lbpl+Lbz

Example-2 error estimation using Model La (full galpy, no error)

import numpy as np import statistics as stat from astropy import units as u from astropy.coordinates import SkyCoord import GalDynPsr

Number=50000 # number of instances to simulate

mul, sigmal = # define mean and standard deviation of l_degree

mub, sigmab = # define mean and standard deviation of b_degree

mud, sigmad = # define mean and standard deviation of d_kpc

now simulate 'Number' instances of l, b and d

lsim = np.random.normal(mul, sigmal, Number)

bsim = np.random.normal(mub, sigmab, Number)

dsim = np.random.normal(mud, sigmad, Number)

declare arrays for excess_parallel, excess_z and excess_total

Lapl = [] Laz = [] totalLa = []

now run GalDynPsr for 'Number' times

count = 0 while (count < Number):

randomly select l, b, d from the simulated values

x=np.random.random_integers(low=1, high=Number-1, size=3)

l=lsim[x[0]] b=bsim[x[1]] d=dsim[x[2]]

for each value of l, b, d, calculate excess_parallel, excess_z and excess_total

Lapl.append(GalDynPsr.modelLa.Expl(l, sigmal, b, sigmab, d, sigmad)) Laz.append(GalDynPsr.modelLa.Exz(l, sigmal, b, sigmab, d, sigmad))

totalLa=np.add(Lapl,Laz)

count = count + 1

now calculate mean and standard deviation of excess_parallel, excess_z and excess_total

print(stat.mean(Lapl), stat.stdev(Lapl), count) print(stat.mean(Laz), stat.stdev(Laz), count) print(stat.mean(totalLa), stat.stdev(totalLa), count)

############################# ######### For a Globular Cluster Pulsar ############################# #Check names of Globular Clusters from GC_Harris.txt and write the name of the cluster. We do not include errors in this case.

Example 3

#Getting Excess_pl for model A for the GC 47Tuc:- #GalDynPsr.ReadGC.gcinput("GC Name") finds the GC in the GC_Harris list and calculates the corresponding galactic coordinates and distance to the GC. Input the arguments to any model as:- (GalDynPsr.ReadGC.ldeg,0.,GalDynPsr.ReadGC.bdeg,0.,GalDynPsr.ReadGC.dkpc,0.)

GalDynPsr.ReadGC.gcinput("47Tuc")

Apl = GalDynPsr.modelA.Expl(GalDynPsr.ReadGC.ldeg,0.,GalDynPsr.ReadGC.bdeg,0.,GalDynPsr.ReadGC.dkpc,0.)

Az = GalDynPsr.modelA.Exz(GalDynPsr.ReadGC.ldeg,0.,GalDynPsr.ReadGC.bdeg,0.,GalDynPsr.ReadGC.dkpc,0.)

Similarly for other models.

#############################

######################### For Intrinsic period derivative calculations ###############################

Parameters required for spin period derivative calcuations:- P = period in seconds, sigP = error in P in second, Pdot_obs = measured period derivative in ss^-1, sigPdot_obs = error in Pdot_obs

#PdotExpl(Excess_pl,P) calculates the planar contribution to the period derivative

#PdotExz(Excess_z,P) calculates the perpendicular contribution to the period derivative

#PdotGal(Excess_pl,Excess_z,P) calculates the sum of planar and perpendicular contributions to the period derivative

#PdotShk(Excess_shk,P) calculates the Shklovskii contribution to the period derivative

#Pdotint(Excess_pl,Excess_z,Excess_shk,Pdot_obs,P) calculates the intrinsic period derivative

#ErrPdotExpl(Excess_pl,error_Excess_pl,P,sigP) calculates the error in the planar contribution to the period derivative

#ErrPdotExz(Excess_z,error_Excess_z,P,sigP) calculates the error in the perpendicular contribution to the period derivative

#ErrPdotGal(Excess_pl,error_Excess_pl,Excess_z,error_Excess_z,P,sigP) calculates the error in the sum of planar and perpendicular contributions to the period derivative

#ErrPdotShk(Excess_shk,error_Excess_shk,P,sigP) calculates the error in the Shklovskii contribution to the spin period derivative

#ErrPdotint(Excess_pl,error_Excess_pl,Excess_z,error_Excess_z,Excess_shk,error_Excess_shk,Pdot_obs,sigPdot_obs,P,sigP) calculates the error in the intrinsic period derivative

Example 4

Say using La

ldeg =20.0 sigl = 2.0

bdeg = 45.0 sigb = 0.2

dkpc = 1.2 sigd = 0.3

mu_alpha = 22.0 sigmua = 0.2

mu_delta = 1.2 sigmud = 0.3

Example spin period and itsderivatives:-

#P = 0.003061844088094608 #sigP = 0.000000000000000015

#Pdot_obs = 0.959013e-20 #sigPdot_obs = 0.000014e-20

#Example orbital period and its derivatives (convert the period and corresponding error in seconds):-

P = 1.19851257518424.03600.0

sigP = 0.00000000001324.03600.0

Pdot_obs = 4.8e-14 sigPdot_obs = 1.1e-14

Ex_pl = GalDynPsr.modelA.Expl(ldeg, sigl, bdeg, sigb, dkpc, sigd)

Ex_z = GalDynPsr.modelA.Exz(ldeg, sigl, bdeg, sigb, dkpc, sigd)

errpl = GalDynPsr.modelA.Errpl(ldeg, sigl, bdeg, sigb, dkpc, sigd)

errz = GalDynPsr.modelA.Errz(ldeg, sigl, bdeg, sigb, dkpc, sigd)

Ex_shk = GalDynPsr.Shk.Exshk(dkpc, sigd, mu_alpha, sigmua, mu_delta, sigmud)

errshk = GalDynPsr.Shk.errShk(dkpc, sigd, mu_alpha, sigmua, mu_delta, sigmud)

PdotExpl = GalDynPsr.pdotint.PdotExpl(Ex_pl,P)

PdotExz = GalDynPsr.pdotint.PdotExz(Ex_z,P)

PdotGal = GalDynPsr.pdotint.PdotGal(Ex_pl,Ex_z,P)

PdotShk = GalDynPsr.pdotint.PdotShk(Ex_shk,P)

Pdotint = GalDynPsr.pdotint.Pdotint(Ex_pl,Ex_z,Ex_shk,Pdot_obs,P)

#Errors in spin period derivatives:-

ErrPdotExpl = GalDynPsr.pdotint.ErrPdotExpl(Ex_pl,errpl,P,sigP)

ErrPdotExz = GalDynPsr.pdotint.ErrPdotExz(Ex_z,errz,P,sigP)

ErrPdotGal = GalDynPsr.pdotint.ErrPdotGal(Ex_pl,errpl,Ex_z,errz,P,sigP)

ErrPdotShk = GalDynPsr.pdotint.ErrPdotShk(Ex_shk,errshk,P,sigP)

ErrPdotint = GalDynPsr.pdotint.ErrPdotint(Ex_pl,errpl,Ex_z,errz,Ex_shk,errshk,Pdot_obs,sigPdot_obs,P,sigP)

#==========================================================================================

Contents of the Package

Datafiles:

parameters.in: Input file contains values of different constants which are subject to change with improvements in measurements. User can change the values of constants if the need be.

GC_Harris.txt: Names and Galactic coordinates of the Globular Clusters.

README.txt: Contents of this README.md file inside package along with code files.

Description of different codes:

read_parameters.py: contains parameters (constants) used in the package

ExcessZ.py: Calculates az/c, both Nice-Taylor (1995) expression as well as our fit of HF2004

Excesspl.py: Calculates apl/c. Can take any value of the slope parameter dv/dR from the parameter file. b0 = (Vs/Rs)*(dV/dR) = 0 is the conventional.

modelA.py: Model A - Calculates apl_old/c (i.e. b0 = 0; DT91) and az_NT95.

modelB.py: Model B - Calculates apl_old/c (i.e. b0 = 0; DT91) and az_HF04.

modelC.py: Model C - Calculates apl_new/c (i.e. dv/dR = -0.2; Reid2014) and az_NT95.

modelD.py: Model D - Calculates apl_new/c (i.e. dv/dR = -0.2; Reid2014) and az_HF04.

modelEa.py: Model Ea - Calculates apl_old/c (i.e. b0 = 0; DT91) and ‘zforce’ in galpy (without BH) modelEb.py: Model Eb - Calculates apl_old/c (i.e. b0 = 0; DT91) and ‘zforce’ in galpy (with BH)

modelFa.py: Model Fa - Calculates apl_new/c (i.e. dv/dR = -0.2; Reid2014) and ‘zforce’ in galpy (without BH) modelFb.py: Model Fb - Calculates apl_new/c (i.e. dv/dR = -0.2; Reid2014) and ‘zforce’ in galpy (with BH)

modelGa.py: Model Ga - Calculates apl/c using Vp'/Vs from galpy (without BH) and az_NT95 modelGb.py: Model Gb - Calculates apl/c using Vp'/Vs from galpy (with BH) and az_NT95

modelHa.py: Model Ha - Calculates ‘Rforce’ in galpy (without BH) and az_NT95 modelHb.py: Model Hb - Calculates ‘Rforce’ in galpy (with BH) and az_NT95

modelIa.py: Model Ia - Calculates apl/c using Vp'/Vs from galpy (without BH) and az_HF04 modelIb.py: Model Ib - Calculates apl/c using Vp'/Vs from galpy (with BH) and az_HF04

modelJa.py: Model Ja - Calculates ‘Rforce’ in galpy (without BH) and az_HF04 modelJb.py: Model Jb - Calculates ‘Rforce’ in galpy (with BH) and az_HF04

modelKa.py: Model Ka - Calculates apl/c using Vp'/Vs from galpy (without BH) and ‘zforce’ in galpy (without BH) modelKb.py: Model Kb - Calculates apl/c using Vp'/Vs from galpy (with BH) and ‘zforce’ in galpy (with BH)

modelLa.py: Model La - Calculates ‘Rforce’ in galpy (without BH) and ‘zforce’ in galpy (without BH) modelLb.py: Model Lb - Calculates ‘Rforce’ in galpy (with BH) and ‘zforce’ in galpy (with BH)

Shk.py: Calculates the Shklovskii term d(mu_T*mu_T)/c. Takes user input of proper motion in right ascension and declination.

galpyMWpl.py: Using vcirc function in galpy to get Vp'(without BH) and eventually apl/c

galpyMWBHpl.py: Using vcirc function in galpy to get Vp'(with BH) and eventually apl/c

galpyMWRfo.py: Using evaluateRforces function in galpy(without BH) to get apl/c

galpyMWBHRfo.py: Using evaluateRforces function in galpy(with BH) to get apl/c

err_excesspl_Damour.py: Calculates error for Excess_pl (DT91)

err_excesspl_Reid.py: Calculates error for Excess_pl (Reid2014)

err_HFhigh.py: Calculates error for Excess_z (our fit for HF04 data)

err_HFlow.py: Calculates error for Excess_z (our fit for HF04 data)

err_NT.py: Calculates error for Excess_z (NT95)

err_Shklovskii.py: Calculates error for Excess_Shk

#==========================================================================================================================================

Though the description of models is provided in the paper, for the sake of completeness we are still providing that here.

Description of Models:-

Model-A: apl_old/c (i.e. b0 = 0; DT91) and Z_NT95 Model-B: apl_old/c (i.e. b0 = 0; DT91) and Z_HF04 Model-C: apl_new/c (i.e. dv/dR = -0.2; Reid2014) and Z_NT95 Model-D: apl_new/c (i.e. dv/dR = -0.2; Reid2014) and Z_HF04

Model-Ea: apl_old/c (i.e. b0 = 0; DT91) and ‘zforce’ in galpy (without BH) Model-Eb: apl_old/c (i.e. b0 = 0; DT91) and ‘zforce’ in galpy (with BH)

Model-Fa: apl_new/c (i.e. dv/dR = -0.2; Reid2014) and ‘zforce’ in galpy (without BH) Model-Fb: apl_new/c (i.e. dv/dR = -0.2; Reid2014) and ‘zforce’ in galpy (with BH)

Model-Ga: Using Vp'/Vs from galpy (without BH) and Z_NT95 Model-Gb: Using Vp'/Vs from galpy (with BH) and Z_NT95

Model-Ha: ‘Rforce’ in galpy (without BH) and Z_NT95 Model-Hb: ‘Rforce’ in galpy (with BH) and Z_NT95

Model-Ia: Using Vp'/Vs from galpy (without BH) and Z_HF04 Model-Ib: Using Vp'/Vs from galpy (with BH) and Z_HF04

Model-Ja: ‘Rforce’ in galpy (without BH) and Z_HF04 Model-Jb: ‘Rforce’ in galpy (with BH) and Z_HF04

Model-Ka: Using Vp'/Vs from galpy (without BH) and ‘zforce’ in galpy (without BH) Model-Kb: Using Vp'/Vs from galpy (with BH) and ‘zforce’ in galpy (with BH)

Model-La: ‘Rforce’ in galpy (without BH) and ‘zforce’ in galpy (without BH) Model-Lb: ‘Rforce’ in galpy (with BH) and ‘zforce’ in galpy (with BH)

Reid2014: Reid M. J., Menten K. M., Brunthaler A., Zheng X. W., Dame T. M., 2014, ApJ, 783, 130.

DT91: Damour T., Taylor J.H., 1991, ApJ, 366, 501.

NT95: Nice D.J., Taylor J.H., 1995, ApJ, 441, 429.

HF04: Holmberg J., Flynn C., 2004, MNRAS, 352, 440.

############################################################