We will introduce a method of plotting using 2-10 keV flux data acquired by a detector called PCA of the RXTE satellite.
RXTE also has an all-sky monitor detector called ASM, but in the case of dark celestial bodies such as AGN (mCrab level), proper data analysis cannot be performed without using the data obtained by PCA pointing observations.
Such long-term data have non-uniform and sparse sample times, but we will also show how to finally interpolate and generate a power spectrum in that case. This isn't enough to handle sparse data, but it's a technique you'll get if you know the trends and trends in quick.
If you can read the code, please refer to my google Colab.
3C273 is the 273rd object in the 3rd edition of the Cambridge Radio Source Catalog, published in 1959, with a redshift z = 0.158. It is a quasar in. Observations have confirmed that there is a black hole in the center, which is 800 million times the mass of the sun, and a powerful jet is ejected from it.
Here, from https://cass.ucsd.edu/~rxteagn/ I will introduce how to do long-term time fluctuation analysis of the famous blazer called 3C273 with google Colab.
There is light curve data in.
plot_3C273
import urllib.request
url="https://cass.ucsd.edu/~rxteagn/3C273/F0210_3C273"
urllib.request.urlretrieve(url, '3C273.txt')
import numpy as np
mjd, flux, err, exp = np.loadtxt("3C273.txt", comments='#', unpack=True)
import matplotlib.pyplot as plt
plt.errorbar(mjd,y=flux,yerr=err)
python
mjd_uniform=np.linspace(np.amin(mjd), np.amax(mjd), 6000) #Divide 6000 equal parts from the beginning to the end. The number 6000 is appropriate, but one day becomes one bottle.
interpolated_y = np.interp(mjd_uniform,mjd,flux) #Spline interpolation is performed to create evenly spaced data.
plt.errorbar(mjd,y=flux,yerr=err, label="raw data")
plt.errorbar(mjd_uniform,y=interpolated_y, label="uniformly interpolated data")
plt.legend()
Spline completion is performed with a function called interp of numpy to generate data at equal time intervals.
python
import matplotlib.mlab as mlab
fNum=len(mjd_uniform)
timebin=(mjd_uniform[1] - mjd_uniform[0]) * 24*60*60 # day to sec
print("timebin = ", timebin, " fNum = ", fNum)
psd2_nofil, freqlist_nofil = mlab.psd(interpolated_y,fNum,1./timebin, sides='onesided', scale_by_freq=True)
plt.xscale("log")
plt.yscale("log")
plt.xlabel(r'Frequency (Hz)')
plt.errorbar(freqlist_nofil, psd2_nofil, fmt='c-', label="FFT for uniformly sampled data")
This is pretty aggressive, but it applies a power spectrum down to low frequencies.
python
#30-day moving average
# https://qiita.com/yamadasuzaku/items/fb744a207e7694d1598d
ave = np.convolve(interpolated_y, np.ones(30)/float(30), 'same')
plt.errorbar(mjd_uniform,y=interpolated_y, label="uniformly interpolated data")
plt.errorbar(mjd_uniform,y=ave, label="30-day average")
plt.legend()
The easiest way to get rid of high frequency components on a moving average is to use numpy's convolve.
python
ave_psd, freq = mlab.psd(ave,fNum,1./timebin, sides='onesided', scale_by_freq=True)
plt.xscale("log")
plt.yscale("log")
plt.xlabel(r'Frequency (Hz)')
plt.errorbar(freqlist_nofil, psd2_nofil, fmt='c-', label="FFT for uniformly sampled data")
plt.errorbar(freqlist_nofil, ave_psd, fmt='r-', label="FFT for uniformly sampled data")
By taking the moving average, it can be confirmed that the high frequency components are dropped.
By using this, it is possible to easily correlate 3C 273 multi-wavelength time series data. Data is quick,
It is published in.
However, since the complement introduced here is interpolated by trapezoidal complement to generate data at regular intervals, it is necessary to carefully consider whether the approximation is valid or not. (There is no problem at the level of looking at the correlation roughly for the time being, so I think it is important to try this quickly once you have data such as experimental data.)
Recommended Posts