背景

RXTE衛星のPCAという検出器によって取得された 2-10 keV のフラックスのデータを用いてプロットする方法を紹介する。

RXTEはASMという全天モニターの検出器もあるが、AGNのような暗い天体の場合(mCrabレベル)は、PCAのポインティング観測で得られたデータを使わないとまともなデータ解析ができない。

このような長期のデータは、サンプル時間が非一様で疎らであるが、その場合にえいやっと内挿してパワースペクトルを生成する方法も紹介する。スパースなデータの扱いはこれでは不十分であるが、quickにトレンドや傾向をみたい場合には知っておいてたら得するテクニックではある。

コードを読めばわかる人は、私の google Colab を参考ください。

3C273のSED (http://isdc.unige.ch/3c273/) のプロットをする場合は、 google Colab 3C273 SED のテキストの読み方を参照ください。

3C273とは？

3C273 は、1959年に出版されたケンブリッジ電波源カタログ第3版の273番目に掲載された天体で、赤方偏移 z = 0.158 にあるクエーサーである。中心には太陽質量の8億倍ものブラックホールが存在しており、そこから強力なジェットが噴出していることが観測から確認されている。

ここでは、https://cass.ucsd.edu/~rxteagn/ の中から、
3C273という有名なブレーザーの長期の時間変動解析をgoogle Colab でやる方法を紹介する。

https://cass.ucsd.edu/~rxteagn/3C273/F0210_3C273

にライトカーブのデータはある。

サンプルコード

データの取得とライトカーブのプロット

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)

データの内挿

mjd_uniform=np.linspace(np.amin(mjd), np.amax(mjd), 6000) # 最初から最後までを6000等分する．6000という数字は適当だけど，1dayが一ビンになるくらい．
interpolated_y = np.interp(mjd_uniform,mjd,flux) # スプライン補間をして，等間隔のデータを作る．
plt.errorbar(mjd,y=flux,yerr=err, label="raw data")
plt.errorbar(mjd_uniform,y=interpolated_y, label="uniformly interpolated data")
plt.legend()

numpyのinterpという関数で、スプライン補完して、等時間間隔のデータを生成する。

パワースペクトルの計算

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")

これでかなり強引であるが、低周波数までのパワースペクトルがかける。

移動平均を取る方法

# 30日の移動平均
# 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()

移動平均で高周波数成分を落とすには、numpyのconvolveを使うのがもっとも簡単にできる。

移動平均と取る前後でのパワースペクトルの比較

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")

移動平均を取ることで高い周波数の成分が落ちていることが確認できる。

おまけ

これを利用すると、3C273の多波長の時系列データの相関なども簡単にとることができる。
データは手っ取り早くは、

http://isdc.unige.ch/3c273/

などに公開されている。3C273のSED (http://isdc.unige.ch/3c273/) のプロットをする場合は、 google Colab 3C273 SED のテキストの読み方を参照ください。

ただし、ここ紹介した補完は、台形補完で内挿して等間隔なデータを生成しているので、その近似が妥当かどうかはよく考えて使う必要がある。(とりあえずざっくりと相関をみる、というレベルでは問題ないので、実験データなどデータを手にしたらこれくらいはサッとやってみるのは大事かとは思う。)