Dynamical Power Spectra (on real data)

[1]:

%matplotlib inline

[2]:

# load auxiliary libraries
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits

# import stingray
import stingray

plt.style.use('seaborn-v0_8-talk')

All starts with a lightcurve..

Open the event file with astropy.io.fits

[3]:

f = fits.open('emr_cleaned.fits')

The time resolution is stored in the header of the first extension under the Keyword TIMEDEL

[4]:

dt = f[1].header['TIMEDEL']

The collumn TIME of the first extension stores the time of each event

[5]:

toa = f[1].data['Time']

Let’s create a Lightcurve from the Events time of arrival witha a given time resolution

[6]:

lc = stingray.Lightcurve.make_lightcurve(toa=toa, dt=dt)

[7]:

lc.plot()
../../_images/notebooks_DynamicalPowerspectrum_DynamicalPowerspectrum_tutorial_%5Breal_data%5D_12_0.png

DynamicPowerspectrum

Let’s create a dynamic powerspectrum with the a segment size of 16s and the powers with a “leahy” normalization

[8]:

dynspec = stingray.DynamicalPowerspectrum(lc=lc, segment_size=16, norm='leahy')

The dyn_ps attribute stores the power matrix, each column corresponds to the powerspectrum of each segment of the light curve

[9]:

dynspec.dyn_ps

[9]:

array([[  2.01901704e+00,   2.32485459e+00,   5.14704363e+00, ...,
          9.76872866e-01,   9.49269045e-01,   4.60522187e+02],
       [  2.93960257e+00,   2.48892516e+00,   3.39280288e+00, ...,
          6.23511732e+00,   4.27550837e+00,   1.06261843e+02],
       [  3.64619904e+00,   1.58266627e+00,   3.42614944e-01, ...,
          1.16952148e+00,   3.54994270e+00,   4.56956463e+01],
       ...,
       [  1.69311108e+00,   5.18784072e-01,   1.57151667e+00, ...,
          1.09923562e+00,   3.40274378e-01,   2.53108287e+00],
       [  2.95675687e-01,   2.47939959e+00,   2.84930818e+00, ...,
          2.99674579e-01,   1.48585951e+00,   7.49068264e+00],
       [  8.84156884e-01,   1.65514790e+00,   4.17385519e-01, ...,
          7.54942692e+00,   9.99801389e-01,   2.03835451e-01]])

To plot the DynamicalPowerspectrum matrix, we use the attributes time and freq to set the extend of the image axis. have a look at the documentation of matplotlib’s imshow().

[10]:

extent = min(dynspec.time), max(dynspec.time), max(dynspec.freq), min(dynspec.freq)

plt.imshow(dynspec.dyn_ps, origin="lower left", aspect="auto", vmin=1.98, vmax=3.0,
           interpolation="none", extent=extent)
plt.colorbar()
plt.ylim(700,850)

[10]:

(700, 850)
../../_images/notebooks_DynamicalPowerspectrum_DynamicalPowerspectrum_tutorial_%5Breal_data%5D_20_1.png 
[11]:

print("The dynamical powerspectrun has {} frequency bins and {} time bins".format(len(dynspec.freq), len(dynspec.time)))

The dynamical powerspectrun has 65535 frequency bins and 104 time bins

# Rebinning in Frequency

[12]:

print("The current frequency resolution is {}".format(dynspec.df))

The current frequency resolution is 0.0625

Let’s rebin to a frequency resolution of 2 Hz and using the average of the power

[13]:

dynspec.rebin_frequency(df_new=2.0, method="average")

[14]:

print("The new frequency resolution is {}".format(dynspec.df))

The new frequency resolution is 2.0

Let’s see how the Dynamical Powerspectrum looks now

[15]:

extent = min(dynspec.time), max(dynspec.time), min(dynspec.freq), max(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower", aspect="auto", vmin=1.98, vmax=3.0,
           interpolation="none", extent=extent)
plt.colorbar()
plt.ylim(500, 1000)

[15]:

(500, 1000)
../../_images/notebooks_DynamicalPowerspectrum_DynamicalPowerspectrum_tutorial_%5Breal_data%5D_29_1.png 
[16]:

extent = min(dynspec.time), max(dynspec.time), min(dynspec.freq), max(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower", aspect="auto", vmin=2.0, vmax=3.0,
           interpolation="none", extent=extent)
plt.colorbar()
plt.ylim(700,850)

[16]:

(700, 850)
../../_images/notebooks_DynamicalPowerspectrum_DynamicalPowerspectrum_tutorial_%5Breal_data%5D_30_1.png

Rebin time

Let’s try to improve the visualization by rebinnin our matrix in the time axis

[17]:

print("The current time resolution is {}".format(dynspec.dt))

The current time resolution is 16.0

Let’s rebin to a time resolution of 64 s

[18]:

dynspec.rebin_time(dt_new=64.0, method="average")

[19]:

print("The new time resolution is {}".format(dynspec.dt))

The new time resolution is 64.0

[20]:

extent = min(dynspec.time), max(dynspec.time), min(dynspec.freq), max(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower", aspect="auto", vmin=2.0, vmax=3.0,
           interpolation="none", extent=extent)
plt.colorbar()
plt.ylim(700,850)

[20]:

(700, 850)
../../_images/notebooks_DynamicalPowerspectrum_DynamicalPowerspectrum_tutorial_%5Breal_data%5D_37_1.png

Trace maximun

Let’s use the method trace_maximum() to find the index of the maximum on each powerspectrum in a certain frequency range. For example, between 755 and 782Hz)

[21]:

tracing = dynspec.trace_maximum(min_freq=755, max_freq=782)

This is how the trace function looks like

[22]:

plt.plot(dynspec.time, dynspec.freq[tracing], color='red', alpha=1)
plt.show()
../../_images/notebooks_DynamicalPowerspectrum_DynamicalPowerspectrum_tutorial_%5Breal_data%5D_42_0.png

Let’s plot it on top of the dynamic spectrum

[23]:

extent = min(dynspec.time), max(dynspec.time), min(dynspec.freq), max(dynspec.freq)
plt.imshow(dynspec.dyn_ps, origin="lower", aspect="auto", vmin=2.0, vmax=3.0,
           interpolation="none", extent=extent, alpha=0.7)
plt.colorbar()
plt.ylim(740,800)
plt.plot(dynspec.time, dynspec.freq[tracing], color='red', lw=3, alpha=1)
plt.show()
../../_images/notebooks_DynamicalPowerspectrum_DynamicalPowerspectrum_tutorial_%5Breal_data%5D_44_0.png

The spike at 400 Hz is probably a statistical fluctutations, tracing by the maximum power can be dangerous!

We will implement better methods in the future, stay tunned ;)