try_b1937.py

"""Test of the scintillometry dedispersion routines on PSR B1937+21"""
import numpy as np
from numpy.lib.format import open_memmap
import astropy.units as u
from baseband import guppi
from scintillometry.dispersion import Dedisperse
from scintillometry.functions import Square
from scintillometry.fourier import get_fft_maker
from scintillometry.integration import Stack
import matplotlib.pylab as plt


FFT = get_fft_maker('pyfftw', flags=['FFTW_ESTIMATE'], threads=4)

dm = 71.02227638 * u.pc / u.cm**3
# from polyco
f0 = (641.928224582360 + 2.78557046583437451/60.) * u.cycle / u.s

fh = guppi.open('puppi_58245_B1937+21_0799.0010.raw')
frequency = (fh.header0['OBSFREQ']
             - fh.header0['OBSBW'] / 2.  # Down to bottom end
             + fh.header0['CHAN_BW'] / 2.  # Up to center of channel 0
             + fh.header0['CHAN_BW'] * np.arange(4)) * u.MHz


def phase(t):
    return ((t - square.start_time) * f0).to(u.cycle)


nbin = 64
dedisperse = Dedisperse(fh, dm, frequency=frequency, sideband=1, FFT=FFT)
square = Square(dedisperse)
stack = Stack(square, nbin, phase)

mm = open_memmap('stack.npy', 'w+', dtype='f4',
                 shape=(7000,) + stack.sample_shape)
stack.read(out=mm)
lc = np.roll(mm, -nbin//3, axis=1)

plt.ion()
plt.plot(np.arange(nbin)/nbin, lc.mean((0, 2, 3)))

ip = lc.mean((2, 3))[:, 46:53].sum(1)
mp = lc.mean((2, 3))[:, 12:19].sum(1)