Welcome to PINT (Pint is not Tempo3)!¶
The power of pulsar timing comes from the ability to precisely predict when pulses will arrive at our telescope. The "timing model" combines parameters describing how the pulsar spins, how the pulsar moves across the sky, how the pulsar orbits around a companion, and how the interstellar medium affects the intrinsic pulse as it propagates through. We can visualize this as follows:
We see that we must account for the position of the telescope on the Earth and as it goes around the Sun. Similarly, we need to account for where the pulsar is in its binary orbit. The arrows show that the two systems are moving with respect to each other in the Galaxy. And, as there is material in the interstellar medium, a single broadband pulse then spreads out in frequencies due to dispersion, such that longer wavelengths (red) arrive later than shorter wavelengths (blue).
A good model (top) will be able to predict when the pulses arrive exactly whereas a bad model (middle) will not. In that case, our mission is to "tweak" the parameters of the timing model until our predictions match the observed pulses (bottom).
The predictions made by the timing model are the result of very arduous mathematical calculations to account for all physical phenomena affecting the pulses (see Lorimer & Kramer 2005):
$ t_{\mathrm{predicted}} = t_{\mathrm{topo}} + t_{\mathrm{corr}} - \Delta D / f^2 + \Delta_{\mathrm{R \odot}} + \Delta_{\mathrm{S \odot}} + \Delta_{\mathrm{E \odot}} + \Delta_{\mathrm{RB \odot}} + \Delta_{\mathrm{SB \odot}} + \Delta_{\mathrm{EB \odot}} + \Delta_{\mathrm{AB \odot}} $
The difference between the observed time of arrival (TOA) of a pulse and the predicted arrival time is known as a timing residual, and plots of residuals help us visualize whether or not we have a good model or not. If we have a good timing model, able to correctly predict the observed times of arrival, the residuals should be randomly distributed close to zero:
GOOD NEWS: Luckily for us, we don't need to perform all these calculations and make the plots by hand (yay!!!) Instead, there are computer codes that do the math-heavy part of the work for us. However, we need to learn how to use them.
PINT (Pint is not Tempo3) is the state-of-art software for performing pulsar analysis. It is based on Object-oriented programming, a computer programming style that organizes software design around data, or objects, rather than functions and logic. An object can be defined as a data field that has unique attributes and behavior.
The structure, or building blocks, of object-oriented programming include the following:
object.method(). Programmers use methods for reusability or keeping functionality encapsulated inside one object at a time.object.attribute.For example, in the world of Pokemon:
First things first, let's load a set of pulsar observations into PINT. Let's import a few basic packages:
import matplotlib.pyplot as plt
from pint.toa import get_TOAs
from pint.models import get_model
from pint.residuals import Residuals
import pint.fitter
import astropy.units as u
import pint.logging
pint.logging.setup("ERROR")
1
Our observations are contained in the file NGC6440E.tim. Let's import them using the get_TOAs function, by doing something life:
my_toas = get_TOAs(path-to-the-file, ephem='DE440')
where path-to-the-file in this case is ./tims/NGC6440E.tim.
# SOME SPACE FOR YOU TO WORK
my_toas =
my_toas will be a PINT object that contains all the information about the observations.Let's take a look at what the observations look like. TOAs are stored in an instance of the TOAs class, which in turn has a table method:
print(my_toas.table['freq', 'mjd', 'error', 'obs'])
These tables allow us to do fancy sorting, selecting, re-arranging very easily. Imagine we want observations with errors smaller than $20~\mu s$. We can first extract the errors by doing
my_toas.get_errors()
and then we can create a "mask" with the observations that fulfill this condition by doing:
my_toas.get_errors() < 20 * u.us
# SOME SPACE FOR YOU TO WORK
select =
print(my_toas.table['freq', 'mjd', 'error', 'obs'][select])
This is a lot of information! Let's print a quick summary:
my_toas.print_summary()
What do you see?
The parameters defining the timing model are in the file NGC6440E.par. Let's import them using the get_model function:
my_timing_model = get_model(path-to-your-file)
# SOME SPACE FOR YOU TO WORK
my_timing_model =
my_timing_model will be a PINT object that contains all the information about the timing model.Let's take a look at what the timing parameters look like:
print(my_timing_model)
Timing models are composed of “delay” terms and “phase” terms, which are computed by the Components of the model.
# This gives a list of all of the component types (so far there are only delay and phase components)
print(my_timing_model.component_types)
# delay_funcs lists all the delay functions in the model, and the order is important!
print(my_timing_model.delay_funcs)
# And phase_funcs holds a list of all the phase functions
print(my_timing_model.phase_funcs)
Each parameter has attributes like the quantity (which includes units). Try printing
my_timing_model.F0.quantity
and
my_timing_model.F0.description
# SOME SPACE FOR YOU TO WORK
print( )
print( )
In a similar fashion, each parameter has attributes that specify the name and type of the parameter, its units, and the uncertainty.
par.quantity and par.uncertainty are both astropy quantities with units. If you need the bare values, access par.value and par.uncertainty_value, which will be numerical values in the units of par.units. Try doing that!# SOME SPACE FOR YOU TO WORK
print( )
Let’s look at those for each of the types of parameters in this model.
printed = []
for p in my_timing_model.params:
par = getattr(my_timing_model, p)
if type(par) in printed:
continue
print("Name ", par.name)
print("Type ", type(par))
print("Quantity ", par.quantity, type(par.quantity))
print("Value ", par.value)
print("units ", par.units)
print("Uncertainty ", par.uncertainty)
print("Uncertainty_value", par.uncertainty_value)
print("Summary ", par)
print("Parfile Style ", par.as_parfile_line())
print()
What do you see?
Plotting time! Let's calculate and then plot the residuals (in microseconds) as a function of the observing time (in Modified Julian Date).
Residual function, which takes two arguments: a TOAs object and a model object:residuals_object = Residuals(TOAs, model)
# SOME SPACE FOR YOU TO WORK
# Hint: remember that the toas are under the variable `my_toas` and the model is under `my_timing_model` !
rs =
What happens if you try to print this new residuals object? Try it!
print(rs)
Note that the resulting Residuals object contains the residuals, but in order to access them we must use the .phase_resids attribute at the end of the residuals object:
Residuals(toas, model).phase_resids
# SOME SPACE FOR YOU TO WORK
rs =
print(rs)
Now let's plot the residuals. First, let extract the TOAs values and the associated errors:
xt = my_toas.get_mjds()
errors = my_toas.get_errors().to(u.us).value
Now do the plot:
# SOME SPACE FOR YOU TO WORK
plt.figure()
plt.errorbar( , yerr=errors, fmt ='o')
plt.title("%s Pre-Fit Timing Residuals" % my_timing_model.PSR.value)
plt.xlabel("MJD")
plt.ylabel("Residual ($\mu s$)")
plt.grid()
Wow, those error bars are way too big. Let's discard the data points with uncertainties $\gt 30.0 \mu s$ - uncertainty estimation is not always reliable when the signal-to-noise is low.
# SOME SPACE FOR YOU TO WORK
error_ok = # Flag all the observations with error bars bigger than 0.1 us
my_new_toas = my_toas[error_ok] # Keep only the ones that pass the check
my_new_toas.print_summary() # Print a summary of the new data set
Now let's calculate and plot residuals again
# SOME SPACE FOR YOU TO WORK
rs =
new_xt =
errors =
plt.plot( , ,"x")
plt.show()
Huh, that's weird. The residuals are not randomly distributed close to zero. Instead, they show some periodic oscillation. I wonder if we could fix that!
In order to fix our model, we need to modify the timing parameters so that they actually fit the observed TOAs. In order to do that, we need to use a fitting algorithm. There are several fitters in $\texttt{PINT}$, each of which is a subclass of Fitter. The most important are:
WLSFitter - PINT’s workhorse fitter, which does a basic Weighted Least-Squares (WLS) minimization of the residuals.
GLSFitter - A Generalized Least-Squares (GLS) fitter, that can handle noise processes like correlated noise and red noise.
Let's use a WLS Fitter to fix our timing model. The return value of most fitters is the final chi^2 value. The general syntax is
pint.fitter.WLSFitter(toas, model)
# SOME SPACE FOR YOU TO WORK
# Instantiate a fitter
wlsfit = pint.fitter.WLSFitter(toas= , model= )
# A fit is performed by calling fit_toas()
wlsfit.fit_toas(maxiter=5)
A summary of the fit and resulting model parameters can easily be printed. Only free parameters will have values and uncertainties in the Postfit column
wlsfit.print_summary()
Let's plot the residuals and compare
plt.figure()
plt.errorbar(
new_xt.value,
wlsfit.resids.time_resids.to(u.us).value,
my_new_toas.get_errors().to(u.us).value,
fmt="x",
)
plt.title("%s Post-Fit Timing Residuals" % my_timing_model.PSR.value)
plt.xlabel("MJD")
plt.ylabel("Residual (us)")
plt.grid()
There also a convenience function for plotting the residuals with the new timing model.
wlsfit.plot()
Now that looks so much better!