Pulsar datareductie

Tammo Jan Dijkema, CAMRAS, 23 oktober 2017

Over dit notebook

Dit is een IPython notebook. Je kunt hem zelf uitvoeren door de bronversie en de gegevens te downloaden:

De gegevens zijn opgenomen met de Dwingeloo Radiotelescoop op 20 oktober 2017 door een groep scholieren die een profielwerkstuk willen schrijven over o.a. de Krabpulsar. Deel van de opdracht was ook een pulsarwaarneming voor Dwingeloo Live, het profielwerkstukprogramma van CAMRAS en ASTRON.

Om de notebook op je eigen computer te draaien heb je nodig:

  • Python
  • Jupyter notebook
  • Matplotlib en Numpy

Dit kun je installeren met bijvoorbeeld Anaconda. Je kunt ook gewoon de resultaten uit dit notebook gebruiken; dan kun je niet zelf met de data spelen.

Het helpt als je een beetje python kan. Als je dat niet kan, wees dan niet bang om her en der een paar waarden aan te passen en te kijken wat er gebeurt.

Je kunt een commando uitvoeren door op Shift-Enter te drukken. Als je op Alt-Enter drukt, wordt de cel geƫvalueerd en wordt een nieuwe cel ingevoegd.

Extra opmerking: voor de grote plaatjes hieronder kan jupyter notebook klagen over IOPub data rate. Om dat op te lossen kun je jupyter starten via jupyter notebook --NotebookApp.iopub_data_rate_limit=10000000000.

Om te beginnen moeten er wat pakketten worden geladen:

In [1]:
import matplotlib.pyplot as plt
from matplotlib import cm
%matplotlib inline
import numpy as np

We slaan vast de kanaalfrequenties (in MHz) alvast op:

In [2]:
frequencies = 441.668359375 - 35.0*np.arange(256)/256.0

B0329+54 (de 'huispulsar')

In [3]:
data = np.load("B0329-clip.npy")

Deze data is al gecorrigeerd voor dispersie, en in losse pulsperiodes ingedeeld. Hiervoor zijn de volgende parameters gebruikt:

In [4]:
dm = 26.8644275665
In [5]:
period = 0.714486089554
In [6]:
date = "Fri Oct 20 09:09:29 CEST 2017"

Met deze periode en DM kunnen we de data over verschillende assen middelen.

In [7]:
data.shape
Out[7]:
(687, 1, 256, 1024)

Deze dataset bevat 687 pulsperiodes, 1 polarisatie, 256 kanalen, met 1024 meetpunten per pulsperiode.

De totale integratietijd is ruim 8 minuten:

In [8]:
687 * period / 60.
Out[8]:
8.1808657253933

We gaan nu deze dataset op verschillende manieren bekijken. Bijvoorbeeld de data voor het derde kanaal van de eerste pulsperiode. Tellen begint voor computers altijd bij 0, : betekent 'alle waarden':

In [9]:
data[0,0,140,:]
Out[9]:
array([58, 61, 61, ..., 55, 65, 70], dtype=uint8)

Deze data kunenn we als volgt in een grafiek zetten:

In [10]:
times = np.linspace(0, period, data.shape[-1])
In [11]:
data_to_plot = data[0,0,140,:]

fig, ax = plt.subplots(1, figsize=(8,6));
ax.plot(times, data_to_plot);
ax.set_title("Eerste pulsperiode, kanaal 140");
ax.set_xlabel("Tijd (s.)");
ax.set_ylabel("Intensiteit (willekeurige schaal)");

We nemen nu het gemiddelde over alle frequenties, voor de eerste pulsperiode. Uit data[0,0,:,:] komt een array met twee assen, namelijk kanaal en meetnummer. Met np.mean(..., axis=(0) nemen we het gemiddelde over de eerste as (kanaal).

In [12]:
data_to_plot = np.mean(data[0,0,:,:], axis=0)

fig, ax = plt.subplots(1, figsize=(8,6));
ax.plot(times, data_to_plot);
ax.set_title("Eerste pulsperiode, gemiddelde alle kanalen");
ax.set_xlabel("Tijd (s.)");
ax.set_ylabel("Intensiteit (willekeurige schaal)");

We kunnen ook het gemiddelde over alle pulsperiodes nemen:

In [13]:
data_to_plot = np.mean(data[:,0,:,:], axis=(0,1))

fig, ax = plt.subplots(1, figsize=(8,6));
ax.plot(times, data_to_plot);
ax.set_title("Gemiddelde alle pulsperiodes, gemiddelde alle kanalen");
ax.set_xlabel("Tijd (s.)");
ax.set_ylabel("Intensiteit (willekeurige schaal)");

We kunnen ook de intensiteit weergeven als kleur, waardoor we de intensiteit per kanaal kunnen zien. We nemen het gemiddelde over alle pulsperiodes en alle polarisaties. Met vmin en vmax kun je het minimum en maximum van de kleurschaal instellen.

In [14]:
data_to_plot = np.mean(data, axis=(0,1));

fig, ax = plt.subplots(1, figsize=(9,6));
spectrumplot = ax.imshow(data_to_plot,
                         #cmap=cm.coolwarm, 
                         origin='upper', interpolation='nearest',
                         extent=([0, period, frequencies[-1], frequencies[0]]),
                         vmin=55, vmax=80
                         );
ax.set_aspect('auto')
ax.set_ylabel('Kanaalfrequentie (MHz)')
ax.set_xlabel('Tijd (s.)')
ax.set_title('Gemiddelde alle pulsperiodes')
cbar = fig.colorbar(spectrumplot,fraction=0.042, pad=0.03);
cbar.ax.set_ylabel('intensiteit (willekeurige eenheid)');

Je ziet hier de puls van de pulsar niet scheeft lopen over frequentie (zoals in Dwingeloo Live), omdat er al voor dispersie is gecorrigeerd.

Pulsar B0329 is zo helder, dat we zelfs in een individuele pulsperiode de puls kunnen zien, bijvoorbeeld de eerste.

In [15]:
data_to_plot = data[0,0,:,:]

fig, ax = plt.subplots(1, figsize=(9,6));
spectrumplot = ax.imshow(data_to_plot,
                         #cmap=cm.coolwarm, 
                         origin='upper', interpolation='nearest',
                         extent=([0, period, frequencies[-1], frequencies[0]]),
                         vmin=55, vmax=80
                         );
ax.set_aspect('auto')
ax.set_ylabel('Kanaalfrequentie (MHz)')
ax.set_xlabel('Tijd (s.)')
ax.set_title('Eerste pulsperiode')
cbar = fig.colorbar(spectrumplot,fraction=0.042, pad=0.03);
cbar.ax.set_ylabel('intensiteit (willekeurige eenheid)');

We kunnen ook de pulsen over tijd laten zien:

In [16]:
data_to_plot = np.mean(data,axis=(1,2))

fig, ax = plt.subplots(1, figsize=(12,100));
timeplot = ax.imshow(data_to_plot, #cmap=cm.PuRd, 
                         origin='lower', interpolation='nearest',
                         extent=([0, period, 0, data_to_plot.shape[0]])
                         );

ax.set_aspect('auto')
ax.set_ylabel('Pulsnummer')
ax.set_xlabel('Tijd (s.)')
ax.set_title('Gemiddelde over alle kanalen')
ax.set_yticks(np.arange(0, data_to_plot.shape[0], 10))
ax.grid()
cbar = fig.colorbar(timeplot, fraction=0.042, pad=0.03);
cbar.ax.set_ylabel('intensiteit (willekeurige eenheid)');