How To

Quick solutions to common tasks (taken from #pint and elsewhere)

Find data files for testing/tutorials

The data files (par and tim) associated with the tutorials and other examples can be located via pint.config.examplefile() (available via the pint.config module):

import pint.config
fullfilename = pint.config.examplefile(filename)

For example, the file NGC6440E.par from the Time a Pulsar notebook can be found via:

import pint
fullfilename = pint.config.examplefile("NGC6440E.par")

Load a par file

from pint.models import get_model
m = get_model(parfile)

Load a tim file

from pint.toa import get_TOAs
t = get_TOAs(timfile)

Load a tim and par file together

from pint.models import get_model_and_toas
m, t = get_model_and_toas(parfile, timfile)

Get the red noise basis functions and the corresponding coefficients out of a PINT fitter object

Select TOAs

You can index by column name into the TOAs object, so you can do toas["observatory"] or whatever the column is called; and that's an array, so you can do toas["observatory"]=="arecibo" to get a Boolean array; and you can index with boolean arrays, so you can do toas[toas["observatory"]=="arecibo"] to get a new TOAs object referencing a subset.

Avoid "KeyError: 'obs_jupiter_pos' error when trying to grab residuals?"

You need to have the TOAs object compute the positions of the planets and add them to the table.

ts.compute_posvels(ephem,planets=True)

This should be done automatically if you load your TOAs with the get_TOAs() or get_model_and_TOAs()

Convert from ELAT/ELONG <-> RA/DEC if I have a timing model

model.as_ICRS(epoch=epoch)

which will give it to you as a model with AstrometryEquatorial components at the requested epoch. Similarly,

model.as_ECL(epoch=epoch)

does the same for AstrometryEcliptic (with an optional specification of the obliquity).

Add a jump programmatically

PINT can handle jumps in the model outside a par file. An example is:

import numpy as np
from astropy import units as u, constants as c
from pint.models import get_model, get_model_and_toas, parameter
from pint import fitter
from pint.models import PhaseJump
import pint.config

m, t = get_model_and_toas(pint.config.examplefile("NGC6440E.par"),
                          pint.config.examplefile("NGC6440E.tim"))

# fit the nominal model
f = fitter.WLSFitter(toas=t, model=m)
f.fit_toas()

# group TOAs: find clusters with gaps of <2h
clusters = t.get_clusters(add_column=True)

# put in the pulse numbers based on the previous fit
t.compute_pulse_numbers(f.model)
# just for a test, add an offset to a set of TOAs
t['delta_pulse_number'][clusters==3]+=3

# now fit without a jump
fnojump = fitter.WLSFitter(toas=t, model=m, track_mode="use_pulse_numbers")
fnojump.fit_toas()


# add the Jump Component to the model
m.add_component(PhaseJump(), validate=False)

# now add the actual jump
# it can be keyed on any parameter that maskParameter will accept
# here we will use a range of MJDs
par = parameter.maskParameter(
    "JUMP",
    key="mjd",
    value=0.0,
    key_value=[t[clusters==3].get_mjds().min().value,
               t[clusters==3].get_mjds().max().value],
    units=u.s,
    frozen=False,
    )
m.components['PhaseJump'].add_param(par, setup=True)

# you can also do it indirectly through the flags as:
# m.components["PhaseJump"].add_jump_and_flags(t.table["flags"][clusters == 3])

# and fit with a jump
fjump = fitter.WLSFitter(toas=t, model=m, track_mode="use_pulse_numbers")
fjump.fit_toas()

print(f"Original chi^2 = {f.resids.calc_chi2():.2f} for {f.resids.dof} DOF")
print(f"After adding 3 rotations to some TOAs, chi^2 = {fnojump.resids.calc_chi2():.2f} for {fnojump.resids.dof} DOF")
print(f"Then after adding a jump to those TOAs, chi^2 = {fjump.resids.calc_chi2():.2f} for {fjump.resids.dof} DOF")
print(f"Best-fit value of the jump is {fjump.model.JUMP1.quantity} +/- {fjump.model.JUMP1.uncertainty} ({(fjump.model.JUMP1.quantity*fjump.model.F0.quantity).decompose():.3f} +/- {(fjump.model.JUMP1.uncertainty*fjump.model.F0.quantity).decompose():.3f} rotations)")

which returns:

Original chi^2 = 59.57 for 56 DOF
After adding 3 rotations to some TOAs, chi^2 = 19136746.30 for 56 DOF
Then after adding a jump to those TOAs, chi^2 = 56.60 for 55 DOF
Best-fit value of the jump is -0.048772786677935796 s +/- 1.114921182802775e-05 s (-2.999 +/- 0.001 rotations)

showing that the offset we applied has been absorbed by the jump (plus a little extra, so chi^2 has actually improved).

See maskParamter documentation on the ways to select the TOAs.

below are several tools i use for new pulsar timing, including cleaning TOAs, adding wraps, and fitting parameters

various packages I use as well as some little convenience functions import numpy as np import matplotlib.pyplot as plt import pint.residuals as res import copy from pint.models import BinaryELL1, BinaryDD, PhaseJump, parameter, get_model from pint.simulation import make_fake_toas_uniform as mft from astropy import units as u, constants as c def dot(l1,l2): return np.array([v1 and v2 for v1,v2 in zip(l1,l2)]) def inv(l): return np.array([not i for i in l])

zapping TOAs on given MJDs, before/after some MJD, or within a window of days

def mask_toas(toas,before=None,after=None,on=None,window=None): cnd=np.array([True for t in toas.get_mjds()]) if before is not None: cnd = dot(cnd,toas.get_mjds().value > before) if after is not None: cnd = dot(cnd,toas.get_mjds().value < after) if on is not None: on=np.array(on) for i,m in enumerate(on): m=mu.day if type(m) is int: cnd = dot(cnd,inv(np.abs(toas.get_mjds()-m).astype(int) == np.abs((toas.get_mjds()-m)).min().astype(int))) else: cnd = dot(cnd,inv(np.abs(toas.get_mjds()-m) == np.abs((toas.get_mjds()-m)).min())) if window is not None: if len(window)!=2: raise ValueError("window must be a 2 element list/array") window = windowu.day lower = window[0] upper = window[1] cnd = dot(cnd,toas.get_mjds() < lower)+dot(cnd,toas.get_mjds() > upper) print(f'{sum(cnd)}/{len(cnd)} TOAs selected') return toas[cnd]

add integer phase wraps between TOAs

def add_wraps(toas,mjd,sign,nwrap=1): cnd = toas.table['mjd'] > Time(mjd,scale='utc',format='mjd') if sign == '-': toas.table['pulse_number'][cnd] -= nwrap elif sign == '+': toas.table['pulse_number'][cnd] += nwrap else: raise TypeError('sign must be "+" or "-"')

plot residuals in phase

def plot_fit(toas,model,track_mode="use_pulse_numbers",title=None,xlim=None,ylim=None): rs=res.Residuals(toas,model,track_mode=track_mode) fig, ax = plt.subplots(figsize=(12,10)) if xlim is not None: ax.set_xlim(xlim) if ylim is not None: ax.set_ylim(ylim) ax.errorbar(rs.toas.get_mjds().value,rs.calc_phase_resids(),
yerr=(rs.toas.get_errors()*model.F0.quantity).decompose().value,fmt='x') ax.tick_params(labelsize=15)

if title is None:
    ax.set_title('%s Residuals, %s toas' %(model.PSR.value,len(toas.get_mjds())),fontsize=18)
else:
    ax.set_title(title,fontsize=18)
ax.set_xlabel('MJD',fontsize=15)
ax.set_ylabel(f'Residuals [phase, P0 = {((1/model.F0.quantity).to(u.ms)).value:2.0f} ms]',fontsize=15)
ax.grid()
return fig, ax

model uncertainties in pulsar phase over a given range of time

def calculate_phase_uncertainties(model, MJDmin, MJDmax, Nmodels=100, params = 'all', error=1*u.us): mjds = np.arange(MJDmin,MJDmax) Nmjd = len(mjds) phases_i = np.zeros((Nmodels,Nmjd)) phases_f = np.zeros((Nmodels, Nmjd)) tnew = mft(MJDmin,MJDmax,Nmjd,model=model, error=error) pnew = {} if params == 'all': params = model.free_params for p in params: pnew[p] = getattr(model,p).quantity + np.random.normal(size=Nmodels) * getattr(model,p).uncertainty for imodel in range(Nmodels): m2 = copy.deepcopy(model) for p in params: getattr(m2,p).quantity=pnew[p][imodel] phase = m2.phase(tnew, abs_phase=True) phases_i[imodel] = phase.int phases_f[imodel] = phase.frac phases = phases_i+ phases_f phases0 = model.phase(tnew, abs_phase = True) dphase = phases - (phases0.int + phases0.frac) return tnew, dphase

plot uncertainties in phase (uses calculate_phase_uncertainties())

def plot_phase_unc(model,start,end,params='all'): if params == 'all': print("calculating phase uncertainty due to all parameters") plab = 'All params' t, dp = calculate_phase_uncertainties(model, start, end) else: if type(params) is list: print("calculating phase uncertainty due to params "+str(params)) plab = str(params) t, dp = calculate_phase_uncertainties(model, start, end, params = params) else: raise TypeError('"params" should be either list or "all"')

plt.gcf().set_size_inches(12,10)
plt.plot(t.get_mjds(),dp.std(axis=0),'.',label=plab)
dt = t.get_mjds() - model.PEPOCH.value*u.d
#plt.plot(t.get_mjds(), np.sqrt((model.F0.uncertainty * dt)**2 + (0.5*model.F1.uncertainty*dt**2)**2).decompose(),label='Analytic')
plt.xlabel('MJD')
plt.ylabel('Phase Uncertainty (cycles)')
plt.legend()

When I'm working on new pulsar solutions, my work flow usually follows something like this:

load pulsar model/TOAs via get_TOAs() and get_model() model = get_model('parfile.par') toas = get_toas('toas.tim')
make a copy of the TOAs that I will edit (for easy resets) newtoas = toas
use plot_fit(newtoas), identify bad toas/missing wraps and mask those data, i.e.

newtoas = mask_toas(newtoas,before=59000,on=59054.321) newtoas.compute_pulse_numbers(model) add_wraps(newtoas,59166,'-') plot_fit(newtoas,model)

once i have the TOAs i want, i fit the data and look at the resids and see how the model changed:

f=WLSFitter(newtoas,model,track_mode='use_pulse_numbers') f.model.free_params=['F0'] # XX TOAs, [lowMJD,highMJD] f.fit_toas() plot_fit(newtoas,f.model) f.print_summary() f.model.compare(model,verbosity='check')

when the model appears to be a good fit, i update the model and add new observations

model=f.model newtoas=mask_toas(toas,before=58000,on=59054.321) plot_fit(newtoas,model)

iterate these last two steps until the fit breaks/i run out of data.

some other, less frequently used snippets

Plot frequency vs residuals

rs=res.Residuals(newtoas,f.model) fig,ax = plt.subplots(figsize=(12,10)) ax.tick_params(labelsize=15) ax.set_ylabel('Frequency [MHz]',fontsize=18) ax.set_xlabel('Phase residuals',fontsize=18)

y = newtoas.get_freqs().to('MHz').value x = rs.calc_phase_resids()

ax.errorbar(x,y,xerr=newtoas.get_errors().to('s').value*f.model.F0.value,elinewidth=2,lw=0,marker='+')

Plot residuals as a function of orbital phase

x = f.model.orbital_phase(newtoas.get_mjds()).value rs=res.Residuals(newtoas,f.model) y = rs.calc_phase_resids()

fig, ax = plt.subplots(figsize=(12,10)) ax.tick_params(labelsize=15) ax.set_xlabel('Orbital Phase',fontsize=18) ax.set_ylabel('Phase Residuals',fontsize=18) ax.grid() for mjd in np.unique(newtoas.get_mjds().astype(int)): cnd = dot(newtoas.get_mjds().astype(int) == mjd,newtoas.get_errors().astype(int) <= 125*u.us) ax.errorbar(x[cnd],y[cnd],yerr=(newtoas.get_errors().to('s')*f.model.F0.quantity).value[cnd],elinewidth=2,lw=0,marker='+',label=mjd.value)

t0 = c.Gc.M_sun/c.c**3 mc = 1.3u.Msun r = t0*mc s = 1

sd = -2rnp.log(1-s*np.sin(x))

ax.plot(x,((sd/c.M_sun)*f.model.F0.quantity).to(''),color='black')
ax.legend(fontsize=15)

adding ELL1 binary component

if 'BinaryELL1' in model.components: model.remove_component('BinaryELL1')

cmp = BinaryELL1() cmp.PB.value = 10 cmp.EPS1.value = 1e-5 cmp.EPS2.value = 1e-5 cmp.TASC.value = 59200 cmp.A1.value = 10

model.add_component(cmp,setup=True)

adding DD binary model

cmp = BinaryDD() cmp.PB.value = 10 cmp.ECC.value = 0.8 cmp.T0.value = 59251. cmp.OM.value = 269. cmp.A1.value = 136.9