Bonus: Wideband timing with PulsePortraiture

This is a stretch goal, mostly for people who already know PSRCHIVE well. It also needs PulsePortraiture's own software environment. If the import cell below fails, that environment isn't set up on this machine — that's expected, and not something you've done wrong.

In Part 1 you made a single template (a 1-D profile) and used pat to get one TOA per observation. But a pulsar's profile changes shape with frequency. PulsePortraiture (Pennucci et al. 2014, 2016) models that frequency evolution as a smooth 2-D portrait and then, for each wideband observation, measures a TOA and a dispersion measure simultaneously.

This notebook runs the standard PulsePortraiture pipeline on the real J1903-7051 data, using the 16-channel archives (data_16ch/) so there is a frequency axis for the model to describe:

1. ppalign — align and average the observations into one master portrait (profile vs frequency).
2. ppspline — fit a smooth spline model of how the profile evolves with frequency (via PCA / eigenprofiles).
3. pptoas — measure wideband TOAs + DMs against that model.

It mirrors the official example_make_model_and_TOAs notebook, but points at real data instead of simulated data.

Setup

PulsePortraiture is a set of importable modules (pplib, ppalign, ppspline, pptoas). They must be importable in this kernel. If they live in a directory that isn't on the path, add it with sys.path (an example location is shown, commented out).

import sys
# If the import below fails with ModuleNotFoundError, uncomment and adjust:
# sys.path.insert(0, "/path/to/PulsePortraiture")

import numpy as np
import matplotlib.pyplot as plt
import psrchive as pr

import pplib
import ppalign as ppa
import ppspline as pps
import pptoas as ppt
from pplib import DataPortrait, make_constant_portrait, write_TOAs
print("PulsePortraiture imported OK")

Inputs

We use the 16-channel J1903-7051 archives and the pulsar ephemeris that shipped with the tutorial. The work products (metafile, portrait, model, TOAs) are written into the current directory.

import glob

REAL      = "real_data/J1903-7051"
ephemeris = f"{REAL}/J1903-7015.par"          # ships with the tutorial
datafiles = sorted(glob.glob(f"{REAL}/data_16ch/*.ar"))
print(f"{len(datafiles)} epochs found")

# PulsePortraiture reads a 'metafile': a plain text list of archive paths.
metafile = "j1903.meta"
with open(metafile, "w") as f:
    f.write("\n".join(datafiles) + "\n")
print("wrote", metafile)

(1) Build an average portrait with ppalign

align_archives aligns every observation (correcting a phase offset and a per-epoch DM) and averages them into a single high-S/N portrait — intensity as a function of frequency and pulse phase. We seed the alignment with the average profile of the first observation.

This is the scripted equivalent of the command line:

ppalign.py -M j1903.meta -I <first archive> -T -C 15.0 -o j1903.port --niter 3

# Make an initial-guess portrait from the first archive's average profile.
dp0 = DataPortrait(datafiles[0], quiet=True)
make_constant_portrait(datafiles[0], "j1903-init.fits",
                       profile=dp0.prof, DM=0.0, dmc=False,
                       weights=None, quiet=True)

portfile = "j1903.port"
ppa.align_archives(metafile=metafile, initial_guess="j1903-init.fits",
                   fit_dm=True, tscrunch=True, pscrunch=True,
                   SNR_cutoff=15.0, outfile=portfile, niter=3, quiet=True)
print("wrote", portfile)

(2) Model the profile evolution with ppspline

ppspline decomposes the portrait into a mean profile plus a few significant eigenprofiles (PCA), and fits a B-spline to how their amplitudes change across the band. The result is a smooth model of profile shape vs frequency.

Command-line equivalent:

ppspline.py -d j1903.port -o j1903.spl -N prof -s

dp = pps.DataPortrait(portfile)
dp.normalize_portrait("prof")          # normalise by the mean profile
dp.show_data_portrait()                 # look at the data we're modelling

# Fit the spline model. Defaults are sensible; a couple shown for clarity.
# (max_ncomp=None lets it choose the number of eigenprofiles automatically;
#  with only 16 channels you'll typically get just 1-2 significant ones.)
dp.make_spline_model(max_ncomp=None, smooth=True, model_name=None, quiet=False)

modelfile = "j1903.spl"
dp.write_model(modelfile, quiet=False)
print("wrote", modelfile)

# Inspect the model: the mean + eigenprofiles, and the spline curves.
dp.show_eigenprofiles()
dp.show_spline_curve_projections()

(3) Measure wideband TOAs with pptoas

pptoas fits the spline model to every observation, returning, per epoch, a TOA at a reference frequency (chosen to be uncorrelated with DM) and a DM measurement. This uses a Fourier-domain phase-gradient fit (Pennucci, Demorest & Ransom 2014).

Command-line equivalent:

pptoas.py -d j1903.meta -m j1903.spl -o j1903_wb.tim

gt = ppt.GetTOAs(metafile, modelfile)
gt.get_TOAs()
print("measured", len(gt.TOA_list), "wideband TOAs")

# Look at how the model fit one observation (data, model, residuals).
gt.show_fit(datafile=gt.datafiles[0], isub=0)

timfile = "j1903_wb.tim"
write_TOAs(gt.TOA_list, SNR_cutoff=0.0, outfile=timfile, append=False)
print("wrote", timfile)

(4) Look at the measurements

Each wideband TOA carries a -pp_dm / -pp_dme flag (the per-epoch DM and its error). We can pull the TOA values straight out of the TOA_list without any timing package.

mjds   = np.array([toa.MJD.in_days() for toa in gt.TOA_list], dtype=float)
errs   = np.array([toa.TOA_error for toa in gt.TOA_list], dtype=float)   # us
dms    = np.array([toa.DM  for toa in gt.TOA_list], dtype=float)
dm_err = np.array([toa.DM_error for toa in gt.TOA_list], dtype=float)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 7), sharex=True)
ax1.errorbar(mjds, dms, dm_err, fmt='k.', ms=4)
ax1.set_ylabel(r"DM [pc cm$^{-3}$]"); ax1.set_title("Wideband DM per epoch")
ax2.plot(mjds, errs, 'b.', ms=4)
ax2.set_xlabel("MJD"); ax2.set_ylabel("TOA uncertainty [us]")
plt.tight_layout(); plt.show()
print("median wideband TOA uncertainty: %.3f us" % np.median(errs))

Where this goes next

j1903_wb.tim is a set of wideband TOAs with simultaneous DM measurements. A timing package (tempo2 with -us/wideband support, or tempo with DMDATA 1 + GLS) can use the DM measurements directly, which is the whole point of wideband timing: you separate the dispersive delay from the rest of the timing model instead of fitting it from narrowband TOAs alone.

Notes / things to try

• With only 16 channels the model is modest; PulsePortraiture really shines on many-channel wideband receivers. Try the same pipeline on a higher nchan dataset if you have one.
• ppgauss is an alternative that models the profile with Gaussian components and can also estimate an unscattered portrait, which lets pptoas --fit_scat fit a scattering timescale per epoch.
• See the official examples in the PulsePortraiture repository for the fake-data generator and a tempo-based residual plot.