Vela-Data Analysis

The analysis of Vela pulsar(PSR B0835-45) observed using the Ooty Radio Telescope.

The Ooty Radio Telescope consists of a cylindrical paraboloid reflecting surface which is $530\;m$ long and $30\;m$ wide, placed on a slope of $11.2 ^\circ$ in N-S direction. The signal detectors consists of an array of 1056 half-wave dipoles which produce phased array in the sky.

The data is nyquist sampled data over 16.5 MHz bandwidth, at a center frequency of 326.5 MHz.

The ASCII data files contains voltage signals in two columns, north and south.

- Devansh Shukla

Setup

from google.colab import drive
# drive.mount('/content/drive')

!pip install numpy pandas lmfit scipy
import numpy as np
import pandas as pd
from scipy.signal import find_peaks
from lmfit import Model
from IPython.display import display, Latex

Requirement already satisfied: numpy in /usr/local/lib/python3.7/dist-packages (1.19.5)
Requirement already satisfied: pandas in /usr/local/lib/python3.7/dist-packages (1.1.5)
Requirement already satisfied: lmfit in /usr/local/lib/python3.7/dist-packages (1.0.3)
Requirement already satisfied: scipy in /usr/local/lib/python3.7/dist-packages (1.4.1)
Requirement already satisfied: python-dateutil>=2.7.3 in /usr/local/lib/python3.7/dist-packages (from pandas) (2.8.2)
Requirement already satisfied: pytz>=2017.2 in /usr/local/lib/python3.7/dist-packages (from pandas) (2018.9)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.7/dist-packages (from python-dateutil>=2.7.3->pandas) (1.15.0)
Requirement already satisfied: uncertainties>=3.0.1 in /usr/local/lib/python3.7/dist-packages (from lmfit) (3.1.6)
Requirement already satisfied: asteval>=0.9.22 in /usr/local/lib/python3.7/dist-packages (from lmfit) (0.9.25)
Requirement already satisfied: future in /usr/local/lib/python3.7/dist-packages (from uncertainties>=3.0.1->lmfit) (0.16.0)

loc = "drive/MyDrive/VelaProcessing/"
filename = "ch00_B0833-45_20150612_191438_011_1"

def read_full(filename):
    f = open(filename, "rb")
    f.seek(0, 2)
    fsize = f.tell()
    f.seek(0, 0)

    north = []
    south = []
    i = 0
    while True:
        n_value, s_value = f.readline().decode().split(" ")
        north.append(n_value)
        south.append(s_value)
        if i%1000000==0: 
          print(i)
        i += 1
        if f.tell() >= fsize:
            print(f"Stopping... {fsize} {f.tell()}")
            return np.array(north, dtype=np.int), np.array(south, dtype=np.int)

# Plotting
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

custom_rcparams = {
    "axes.labelsize": 10,
    "axes.titlesize": 10,
    "axes.grid": True,
    'axes.edgecolor': "white",
    "text.color" : "white",
    "xtick.color" : "white",
    "ytick.color" : "white",
    "axes.facecolor": "black",
    "axes.labelcolor" : "white",
    # Figure
    "figure.facecolor": "black",
    "figure.autolayout": True,
    "figure.titlesize": 9,
    "figure.dpi": 150,
    "savefig.format": "pdf",
    "lines.linewidth": 1,
    # Legend
    "legend.fontsize": 8,
    "legend.frameon": True,
    # Ticks
    "xtick.minor.visible": False,
    "xtick.direction": "in",
    "ytick.direction": "in",
    "ytick.minor.visible": True,
    # TeX
    # "pgf.texsystem": "lualatex",
}
mpl.rcParams.update(custom_rcparams)

def plot(x, limit=1e6, cmap="twilight_r", xoffset=None, ts=None, title=None, xlabel="Time(bins)", ylabel="Freq. channels", grid=True):
    fig, ax = plt.subplots()
    ax.grid(grid)
    plt.imshow(x.real.transpose(), cmap=cmap, aspect="auto", origin="lower")
    plt.colorbar()
    plt.clim(-limit, limit)
    ax.set_ylabel(ylabel)
    ax.set_xlabel(xlabel)
    if title:
      ax.set_title(title)
    if xoffset is not None:
      xticks = ax.get_xticks()
      print(xticks)
      print([str(round((xoffset + x), 2)) for x in xticks])
      xticklabels = [str(round((xoffset + x)*ts, 2)) for x in xticks]
      print(xticklabels)
      ax.set_xticklabels(xticklabels)
      ax.set_xlabel("Time(s)")
    return

def compute_fft(north, south, channels=512, avg=60):
  # computing fft
  channelsd2 = int(channels/2)
  size = int(north.size/(channels*avg))
  ssize = int(north.size/channels)

  dyn_n = np.zeros((size, channelsd2), dtype=complex)
  dyn_s = np.zeros((size, channelsd2), dtype=complex)
  m=0
  for i in range(0, size):
      sum_n = np.zeros(channelsd2, dtype=complex)
      sum_s = np.zeros(channelsd2, dtype=complex)
      for j in range(0, avg):
          north_chunk = north[channels*m : channels*(m+1)]
          south_chunk = south[channels*m : channels*(m+1)]
          sum_n += np.abs(np.fft.fft(north_chunk, channels)[channelsd2:])**2
          sum_s += np.abs(np.fft.fft(south_chunk, channels)[channelsd2:])**2
          m += 1
      dyn_n[i] = sum_n/avg
      dyn_s[i] = sum_s/avg
  return dyn_n, dyn_s

def params(channels, avg):
    global cf, bw, fs
    ts = channels*avg/fs
    cw = bw*2/channels
    start_freq = 326.5e6 - cw*channels/4
    stop_freq = 326.5e6 + cw*channels/4
    return start_freq, stop_freq, cw, ts

Initial parameters

Center frequency, $f_c = 326.5\;$MHz

Bandwidth, $bw = 16.5\;$MHz

Sampling freq, $f_s = 33\;$MHz

cf = 326.5e6  # 326.5 MHz
bw = 16.5e6  # 16.5 MHz
fs = 33e6  # 33 MHz

Reading the data

north, south = read_full(loc + filename)

0
1000000
2000000
3000000
4000000
5000000
6000000
7000000
8000000
9000000
10000000
11000000
12000000
13000000
14000000
15000000
16000000
17000000
18000000
19000000
20000000
21000000
22000000
23000000
24000000
25000000
26000000
27000000
28000000
29000000
30000000
31000000
32000000
Stopping... 211985964 211985964

Time-Domain Properties | Signal Statistics

Voltage signal characterstics

We expect the signal to have gaussian distribution with mean $\mu$ and standard deviation $\sigma$

The mean will not be zero as one would expect, instead telescope-back-end induces a small bias, thus adding an offset to all voltages.

north_mean = np.mean(north)
north_std = np.std(north)
south_mean = np.mean(south)
south_std = np.std(south)

print("- North")
display(Latex(rf"$\mu_{{\text{{north}}}} = {np.round(north_mean, 2)}$"))
display(Latex(rf"$\sigma_{{\text{{north}}}} = {np.round(north_std, 2)}$"))
print("- South")
display(Latex(rf"$\mu_{{\text{{south}}}} = {np.round(south_mean, 2)}$"))
display(Latex(rf"$\sigma_{{\text{{south}}}} = {np.round(south_std, 2)}$"))

- North

$\mu_{\text{north}} = 3.48$

$\sigma_{\text{north}} = 28.03$

- South

$\mu_{\text{south}} = 0.78$

$\sigma_{\text{south}} = 29.92$

Now, for verifing the gaussian nature of the pulsar signal, we select 100,000 randomly selected data points and observe them in a histogram plot.

random_samples_north = np.random.choice(north, size=100000)
random_samples_south = np.random.choice(south, size=100000)
fig = plt.figure(figsize=(8,4))
gs = gridspec.GridSpec(1, 2)
ax = fig.add_subplot(gs[0, 0])
_ = plt.hist(random_samples_north, bins=50)
ax.set_xlabel("bins")
ax.set_ylabel("No of occurances")
ax.set_title("North")
ax = fig.add_subplot(gs[0, 1])
_ = plt.hist(random_samples_south, bins=50)
ax.set_xlabel("bins")
ax.set_ylabel("No of occurances")
ax.set_title("South")

Text(0.5, 1.0, 'South')

As expected we observe the gaussian nature in the data.

Power signal characterstics

We now investigate the properties of the signal intensity, i.e. the square of the voltage data.

$$P \propto V^2$$

We expect the signal intensity to have an exponential distribution with equal mean and standard deviation.

signal_intensity_north = north**2
signal_intensity_south = south**2

random_samples_intensity_north = np.random.choice(signal_intensity_north, size=100000)
random_samples_intensity_south = np.random.choice(signal_intensity_south, size=100000)
fig = plt.figure(figsize=(8,4))
gs = gridspec.GridSpec(1, 2)
ax = fig.add_subplot(gs[0, 0])
_ = plt.hist(random_samples_intensity_north, bins=50)
ax.set_xlabel("bins")
ax.set_ylabel("No of occurances")
ax.set_title("North")
ax = fig.add_subplot(gs[0, 1])
_ = plt.hist(random_samples_intensity_south, bins=50)
ax.set_xlabel("bins")
ax.set_ylabel("No of occurances")
ax.set_title("South")
plt.suptitle("Signal Intensity", fontsize=12)

Text(0.5, 0.98, 'Signal Intensity')

Frequency domain properties

Power spectrum

channels = 512
avg = 1
north_freq, south_freq = compute_fft(north, south, channels=channels, avg=avg)

fig = plt.figure(figsize=(9,4))
ax = fig.add_subplot(gs[0, 0])
_ = ax.plot(np.mean(north_freq, axis=0).real)
ax.set_xlabel("Frequency Bins")
ax.set_ylabel("Intensity")
ax.set_title("North")
ax = fig.add_subplot(gs[0, 1])
_ = ax.plot(np.mean(south_freq, axis=0).real)
ax.set_xlabel("Frequency Bins")
ax.set_ylabel("Intensity")
ax.set_title("South")
plt.suptitle("Power spectrum", fontsize=12)

Text(0.5, 0.98, 'Power spectrum')

The sharp peaks in the power spectrum represents the Radio-Frequency interferences.

Note: The band of the filter used in the telescope back-end is retained in the data.

Dyanmic Spectrum

Frequency vs Time spectrum

Initially we choose 512 channels and average 60 data points to improve the SNR.

Averaging 60 data points (avg=60) gives we an integration time of about 1ms

$\text{Integration time} = \dfrac{channels \times avg}{fs}$

channels = 512
avg = 60
dyn_n, dyn_s = compute_fft(north, south, channels=channels, avg=avg)

The data contains some interesting signal which repeats after an interval; it first appears in higher frequencies and gradually appears later in lower freqs.

(figure shown below)

plot(dyn_n, title="North")

plot(dyn_s, title="South")

Dispersion Measure

It is integrated column density of free electrons between the observer and the pulsar.

It manifests itself as the broading of an otherwise sharp peak when observed under finite bandwidth.

In presence of charged particles, the electrostatic interaction b/t the light and charged particles causes a delay in the propogation of light. This delay is inversely proportional to the frequency; so low frequencies experiences a larger delay compared to high frequencies.

$$delay \propto \dfrac{1}{frequency}$$$$\Delta t = 4.149 \times 10^{-3} DM \left(\dfrac{1}{\nu_1^2(\text{GHz})} - \dfrac{1}{\nu_2^2(\text{GHz})}\right)$$

This DM can be approximated by observing the position of peaks in various frequency channels and performing a linear-fit.

$$ t = 4.149 \times 10^{-3} \;DM \dfrac{1}{\nu^2}$$

But, it requires a much-better SNR; single pulses have poor SNR so will have to increase the channel width while making sure the remaining dispersion doesn't spoil the signal.

Here we use 64 channels i.e. channel width of 515.625 kHz. Since we are unaware of the pulsar period, we cannot improve the SNR further; now we choose some strong pulses in the data and fit a gaussian curve; further, we use the mean of the fitted-gaussian to obtain the arrival-time of the individual pulses.

Here we are assuming that the signal is gaussian in nature.

Now, we fit the arrival-time vs $1/\nu^2$ with a linear-fit to obtain the slope and using the above equation, we have an approximate estimate of the dispersion measure in that direction.

channels = 64
avg = 60
dyn_n, dyn_s = compute_fft(north, south, channels=channels, avg=avg)

plot(dyn_s, title="South", limit=1e5)

def plot_pulses(a, b, chs, arr):
    for ch in chs:
      x = [*range(a, b)]
      ydata = arr[a:b, ch].real
      plt.figure(figsize=(4,3))
      plt.plot(x, ydata)
      plt.title(f"ch={ch}")
    return
plot_pulses(3400, 4200, [15,16,17,18], dyn_s)

def fit_gaussian(a, b, chs, arr, center=0):
    def gauss(x, A, B, C, D):
        # B is mean; C is std
        return A*np.exp(-((x-B)**2/(2*C**2))) + D
    pulse_gauss = []
    for ch in chs:
      x = np.array([*range(a, b)])
      ydata = arr[a:b, ch].real
      gmodel = Model(gauss)
      result = gmodel.fit(ydata, x=x, A=1, B=center, C=50, D=0.0)
      # print(result.fit_report())
      if result:
          pulse_gauss.append([ch, result.best_values.get("B"), np.sqrt(result.covar[1][1])])  # err=sqrt of diagonal of cov(1, 1)

      plt.figure(figsize=(6,4))
      plt.plot(x, ydata)
      plt.plot(x, result.best_fit, '-', label='best fit')
      plt.errorbar(x, result.best_fit, yerr=result.eval_uncertainty())
      plt.title(f"Pulse (ch={ch})")
      plt.xlabel("Time(bins)")
      plt.ylabel("Signal Intensity")
    return pulse_gauss

pulse_gauss = np.array(fit_gaussian(3400, 4200, chs=[15, 16, 17, 18], arr=dyn_s, center=3800))

start_freq, stop_freq, cw, ts = params(64, 60)

pulse_freq = start_freq + pulse_gauss[:, 0] * cw
pulse_gauss_time = pulse_gauss[:, 1] * ts
pulse_gauss_time_err = pulse_gauss[:, 2] * ts

xdata = (pulse_freq/1e9)**(-2)
ydata = pulse_gauss_time
yerr = pulse_gauss_time_err

def straight_line(x, m, c):
  return m*x + c

gmodel = Model(straight_line)
result = gmodel.fit(ydata, x=xdata, m=1, c=0, weights=1/yerr)
print(result.fit_report())

plt.errorbar(xdata, ydata, yerr=yerr, marker="o", markersize=3.5, linestyle="", ecolor="red")
plt.plot(xdata, result.best_fit, "-", color="C1", label="best_fit")
plt.title("DM Estimation")
plt.xlabel(r"$1/{\nu}^2 ({GHz}^{-2}$)")
plt.ylabel(r"time(s)")

[[Model]]
    Model(straight_line)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 7
    # data points      = 4
    # variables        = 2
    chi-square         = 2.46817733
    reduced chi-square = 1.23408866
    Akaike info crit   = 2.06874238
    Bayesian info crit = 0.84133110
[[Variables]]
    m:  0.27349646 +/- 0.00859489 (3.14%) (init = 1)
    c: -2.11640825 +/- 0.08049609 (3.80%) (init = 0)
[[Correlations]] (unreported correlations are < 0.100)
    C(m, c) = -1.000

Text(0, 0.5, 'time(s)')

slope = result.best_values.get("m")
slope_err = np.sqrt(result.covar[0][0])
DM = slope * 1e3/4.149
DM_err = slope_err * 1e3/4.149

display(Latex(f"\\text{{Dispersion Measure = }}{np.round(DM, 4)} \pm {np.round(DM_err, 4)}\; pc/cc"))

\text{Dispersion Measure = }65.9186 \pm 2.0716\; pc/cc

Distance computation

The distance to the pulsar can be computed by using the defination of the dispersion measure.

$$DM = \int_0^L n_e dl $$$$L = \dfrac{DM}{|n_e|}$$

ne = 0.023  # cc-1 (THE DISTANCE TO THE VELA PULSAR GAUGED WITH HST PARALLAX OBSERVATIONS, 10.1086/323377)
distance = DM/ne
distance_err = DM_err/ne

display(Latex(f"\\text{{Distance to the pulsar }}= {np.round(distance, 4)} \pm {np.round(distance_err, 4)}\; pc"))

\text{Distance to the pulsar }= 2866.0281 \pm 90.0677\; pc

Radio Frequency Interference mitigation

RFI mitigation plays a very important role in extracting the meaning-full signal among the noises.

channels = 512
avg = 60
dyn_n, dyn_s = compute_fft(north, south, channels=channels, avg=avg)

def reject_RFI(data_n, data_s, avg, channels, bound):
    efficiency_matrix_n = (np.mean(data_n, axis=0)/(np.std(data_n, axis=0) * np.sqrt(avg))).real
    efficiency_matrix_s = (np.mean(data_s, axis=0)/(np.std(data_s, axis=0) * np.sqrt(avg))).real

    bounds_n = [np.mean(efficiency_matrix_n) - bound*np.std(efficiency_matrix_n), np.mean(efficiency_matrix_n) + bound*np.std(efficiency_matrix_n)]
    bounds_s = [np.mean(efficiency_matrix_s) - bound*np.std(efficiency_matrix_s), np.mean(efficiency_matrix_s) + bound*np.std(efficiency_matrix_s)]

    rejected_channels = [
                         np.where((efficiency_matrix_n < bounds_n[0]) | (efficiency_matrix_n > bounds_n[1]))[0],
                         np.where((efficiency_matrix_s < bounds_n[0]) | (efficiency_matrix_s > bounds_n[1]))[0]
    ]

    plt.figure(figsize=(5,3.5))
    plt.plot(efficiency_matrix_n, color="C0", label="north")
    plt.plot(efficiency_matrix_s, color="C1", label="south")
    plt.hlines(bounds_n[0], 0, int(channels/2) + 5, linestyle="--", color="red", label=rf"${bound} \sigma\; north$")
    plt.hlines(bounds_n[1], 0, int(channels/2) + 5, linestyle="--", color="red")
    plt.hlines(bounds_s[0], 0, int(channels/2) + 5, linestyle="--", color="purple", label=rf"${bound} \sigma\; south$")
    plt.hlines(bounds_s[1], 0, int(channels/2) + 5, linestyle="--", color="purple")
    plt.xlim(0, int(channels/2) + 5)
    plt.ylim(0.8, 1.1)
    plt.xlabel("Frequency Channels")
    plt.ylabel("Efficiency")
    plt.title(f"Efficiency")
    plt.legend()

    mask_n = np.zeros(int(channels/2), dtype=np.int)
    mask_n[rejected_channels[0]] = 1
    mask_s = np.zeros(int(channels/2), dtype=np.int)
    mask_s[rejected_channels[1]] = 1

    new_efficiency_n = np.ma.array(efficiency_matrix_n, mask=mask_n)
    new_efficiency_s = np.ma.array(efficiency_matrix_s, mask=mask_s)

    plt.figure(figsize=(5,3.5))
    plt.plot(new_efficiency_n, color="C0", label="north")
    plt.plot(new_efficiency_s, color="C1", label="south")
    plt.hlines(bounds_n[0], 0, int(channels/2) + 5, linestyle="--", color="red", label=rf"${bound} \sigma\; north$")
    plt.hlines(bounds_n[1], 0, int(channels/2) + 5, linestyle="--", color="red")
    plt.hlines(bounds_s[0], 0, int(channels/2) + 5, linestyle="--", color="purple", label=rf"${bound} \sigma\; south$")
    plt.hlines(bounds_s[1], 0, int(channels/2) + 5, linestyle="--", color="purple")
    plt.xlim(0, int(channels/2) + 5)
    plt.ylim(0.8, 1.1)
    plt.xlabel("Frequency Channels")
    plt.ylabel("Efficiency")
    plt.title(f"Efficiency after RFI-mitigation")
    plt.legend()

    processed_data_n = np.copy(data_n)
    processed_data_n[:, rejected_channels[0]] = np.nan
    processed_data_s = np.copy(data_s)
    processed_data_s[:, rejected_channels[1]] = np.nan
  
    return rejected_channels, processed_data_n, processed_data_s

rejected_channels, dyn_n, dyn_s = reject_RFI(dyn_n, dyn_s, avg, channels, bound=1.25)

plot(dyn_n, title="North", grid=False)
plot(dyn_s, title="South", grid=False)

Dedispersion

# DM pc/cc
start_freq, stop_freq, cw, ts = params(channels, avg)

f1 = start_freq/1e9
channel_freq = np.array([start_freq + cw*i for i in range(0, int(channels/2))])/1e9
time_delay = 4.149*10**(-3) * DM * (f1**(-2) - channel_freq**(-2))  # s

pack_delay = []
for i in range(0, len(time_delay)):
  pack_delay.append(time_delay[i]/ts)
packs = np.array(pack_delay, dtype=np.int)

transposed_dyn = dyn_n.transpose()

def shift_fill(x, shifts, channels, fill=0):
    to_pad = max(shifts + shifts[::-1])
    arr = []
    for i in range(0, int(channels/2)):
      arr.append(np.pad(x[i], (shifts[i], to_pad - shifts[i]), "constant", constant_values=fill))
    return np.array(arr)
dedispersed_dyn = shift_fill(transposed_dyn, packs, channels)

plot(dedispersed_dyn.transpose(), limit=1e6, title="Dedispersed Dynamic Spectrum")

pulses_to_fold = np.nanmean(dedispersed_dyn.real.transpose(), 1)
fig = plt.figure()
plt.plot(pulses_to_fold)
plt.xlim(0, 1200)
plt.ylim(bottom=np.nanmean(pulses_to_fold))
plt.title("Dedipsersed Pulses")
plt.xlabel("Time(bins)")
plt.ylabel("Signal Intensity")

Text(0, 0.5, 'Signal Intensity')

gs = gridspec.GridSpec(2, 1, height_ratios=[2, 1])
fig = plt.figure(figsize=(6, 6))

ax1 = fig.add_subplot(gs[0, 0])
ax1.set_aspect(2)
im = ax1.imshow(dedispersed_dyn.real, cmap="twilight_r", origin="lower", aspect="auto")
plt.colorbar(im, aspect=50, orientation='horizontal')
im.set_clim(-1e6, 1e6)
ax1.set_title("Dedispersed dynamic spectrum")

ax2 = fig.add_subplot(gs[1, 0])
ax2.get_shared_x_axes().join(ax1, ax2)
ax2.plot(pulses_to_fold)
ax2.set_xlim(left=0)
ax2.set_ylim(bottom=np.nanmean(pulses_to_fold))
ax2.set_title("Dedispersed pulses")

Text(0.5, 1.0, 'Dedispersed pulses')

The Pulsar time-period

def identify_pulse_peaks(a, b, arr, std, width=5):
    peakfinder = pulses_to_fold[a: b]
    peaks, _ = find_peaks(peakfinder, height=peakfinder.mean() + std*peakfinder.std(), width=width)
    fig = plt.figure()
    plt.plot(peakfinder)
    plt.plot(peaks, peakfinder[peaks], "o", markersize="1.5", color="red")
    plt.xlim(0, 1000)
    plt.ylim(bottom=peakfinder.mean() - std*peakfinder.std())
    plt.title("Dedipsersed Pulses")
    plt.xlabel("Time(bins)")
    plt.ylabel("Signal Intensity")
    diff_peaks = [peaks[i] - peaks[i-1] for i in range(1, len(peaks))]
    return peaks, np.mean(diff_peaks), len(diff_peaks)

pulse_peaks, time_period, no_of_pulses = identify_pulse_peaks(200, 1100, pulses_to_fold, std=1.5)
time_period = time_period*ts

print(f"Approximate time period = {np.round(time_period*1000, 4)} ms; with {no_of_pulses} peaks")

Approximate time period = 89.5003 ms; with 7 peaks

def straight_line(x, m, c):
    return m*x + c

gmodel = Model(straight_line)
result = gmodel.fit(data=pulse_peaks, x=[*range(0, len(pulse_peaks))], m=1, c=0)
print(result.fit_report())

plt.plot(pulse_peaks, "o", markersize=4, color="C0", label="peaks")
plt.plot(result.best_fit, "--", color="C1", label="fit")
plt.xlim(left=0)
plt.ylim(bottom=0)
plt.xlabel("Pulse peaks")
plt.ylabel("Time(bins)")
plt.legend(loc="upper left")

[[Model]]
    Model(straight_line)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 6
    # data points      = 8
    # variables        = 2
    chi-square         = 2.82142857
    reduced chi-square = 0.47023810
    Akaike info crit   = -4.33758560
    Bayesian info crit = -4.17870251
[[Variables]]
    m:  96.0357143 +/- 0.10581184 (0.11%) (init = 1)
    c:  137.000000 +/- 0.44264268 (0.32%) (init = 0)
[[Correlations]] (unreported correlations are < 0.100)
    C(m, c) = -0.837

<matplotlib.legend.Legend at 0x7fe3431d4790>

slope_time_period = result.best_values.get("m")
slope_time_period_err = np.sqrt(result.covar[0][0])
time_period = slope_time_period*ts
time_period_err = slope_time_period_err*ts

display(Latex(f"\\text{{Time Period }} = {np.round(time_period*1000, 4)} \pm {np.round(time_period_err*1000, 4)}\;ms"))
display(Latex(f"\\text{{\% error }} = {np.round(np.round(time_period_err*1000, 4)*100/np.round(time_period*1000, 4), 4)} \;\%"))

\text{Time Period } = 89.4005 \pm 0.0985\;ms

\text{\% error } = 0.1102 \;\%

Finally, using the approximate pulsar period, we can fold the pulses to obtain an integrated profile for the pulsar.

def fold_prepare(arr, pulses, cut_intervals):
    plt.plot(arr)
    plt.xlim(80, 1200)
    plt.ylim(arr.mean(), arr.mean() + arr.std())
    plt.vlines(cut_intervals, 0, arr.max(), linestyles="--", color="red")
    plt.title("Dedispersed Pulses")
    plt.xlabel("Time(bins)")
    plt.ylabel("Signal Intensity")
    return

def after_fold(slices, shift, offset):
    fig = plt.figure()
    plt.xlim(0, len(slices[0]))
    plt.ylim(offset, offset + (len(slices)+2)*shift)
    for i in range(0, len(slices)):
      plt.plot(slices[i] + shift*i)
    plt.title("Stack plot")
    fig = plt.figure()
    plt.stackplot([*range(0, len(slices[0]))], slices)
    plt.xlim(left=0)
    plt.title("Stack plot")
    return

def fold(arr, pulses, start, interval, confirm=False):
    fig, ax = plt.subplots()
    cut_intervals = [start + interval*i for i in range(0, pulses)]
    fold_prepare(arr, pulses, cut_intervals)
    if confirm:
      slices = [arr[cut_intervals[i-1]: cut_intervals[i]] for i in range(1, pulses)]
      slices = np.array(slices)
      return slices
    return []

print(f"interval={round(time_period/ts, 0)}")
slices = fold(pulses_to_fold, pulses=9, start=300, interval=96, confirm=True)

interval=96.0

after_fold(slices, shift=4e4, offset=4e5)

folded_pulse = np.mean(slices, axis=0)
plt.plot(folded_pulse)
plt.title("Folded profile")
plt.xlabel("Time(bins)")
plt.ylabel("Signal Intensity")

Text(0, 0.5, 'Signal Intensity')

Note that the average profile resembles a gaussian function, hence our initial assumption was fairly accurate.

print("Results: \n")
display(Latex(f"\\text{{Dispersion measure = }}{np.round(DM, 4)} \pm {np.round(DM_err, 4)}\; pc/cc"))
display(Latex(f"\\text{{Distance to the pulsar }}= {np.round(distance, 4)} \pm {np.round(distance_err, 4)}\; pc"))
display(Latex(f"\\text{{Time period }} = {np.round(time_period*1000, 4)} \pm {np.round(time_period_err*1000, 4)}\;ms"))

Results: 

\text{Dispersion measure = }65.9186 \pm 2.0716\; pc/cc

\text{Distance to the pulsar }= 2866.0281 \pm 90.0677\; pc

\text{Time period } = 89.4005 \pm 0.0985\;ms

Recommended reads

• Kishalay's Vela analysis report
• Pulsar data analysis with PSRCHIVE (DOI: 1205.6276)
• D.R. Lorimer & M.Kramer, Handbook of Pulsar Astronomy, Vol.4, Cambridge University Press (2004)
• Pulsars: https://www.cv.nrao.edu/~sransom/web/Ch6.html