Back

RHESSI データの解析方法 II

解析方法 II(マニアック版へ)

このマニュアルは、 RHESSI解析マニュアルに準じて書きました。
Object Syntaxで作られているので、そのように書くようにしました。


(0) 観測データの取得

http://hesperia.gsfc.nasa.gov/hessidata/
ftp://hercules.ethz.ch/pub/hessi/data/
のどちらかから取得します。

ディレクトリー構造は YYYY -> MM -> DD となっており、 欲しい日付までたどり着いたら、hsi_yyyymmdd_hhmmss_vvv.fits というfileがあるので、欲しい時間帯のfileをダウンロードする。
*観測fileの形式はFITSです。
*将来的にはフレアリストに自動でアクセスして、データ取得の手間は なくなる予定です。


(1) ライトカーブを描く

現在はフレアリストが立ち上がっていないので、像合成する時間などは、 ライトカーブなどから決めてやる必要があります。
GOESなどから時間を指定しても良いですが、硬X線(RHESSI)のデータから ライトカーブを描く方法です。

IDL> obj_ltc=hsi_lightcurve()
obj_ltcというobjectを定義。

IDL> obj_ltc->plot, obs_time_interval=['2002/04/04 15:00:00','2002/04/04 16:00:00']
ライトカーブを描かせたい時間帯を指定。更に、エネルギー帯を指定したりするには、
IDL> obj_ltc->plotman, obs_time_interval=['2002/04/04 15:00:00','2002/04/04 16:00:00'], $
ltc_energy_band=[25,50,100,150],ltc_time_resolution=1.0

ここでplotではなく、plotmanにしたのはplot managerを立ち上げるため (色別で描いてくれたりする)。
ltc_time_resolutionはplotする時間を1秒間隔にする。


(2) 像合成

IDL> obj_im = hsi_image(obs_time_interval=['2002/04/04 15:36:10', '2002/04/04 15:37:00'])
IDL> pmtras_analysis,'hsi_20020404_153600_002.fit',/no_strid

IDL> obj_im -> set,as_roll_solution='PMT',ras_time_extension=[-500,500]
IDL> obj_im -> set,energy_band=[25.,50.],image_dim=[64,64], $
IDL> pixel_size=[32,32],time_range=[0,21.6645],xyoffset=[0,0], $
IDL> det_index_mask=[0,0,0,0,1,1,1,1,0],as_spin_period=4.3329
IDL> im = obj_im -> getdata()
IDL> obj_im -> plotman

IDL> obj_im -> set,image_dim=[64,64],pixel_size=[4,4], $
IDL> xyoffset=[-680,-680],det_index_mask=[0,0,1,1,1,1,1,0,0]
IDL> im = obj_im -> getdata()
IDL> obj_im -> plotman

(refining)
IDL> obj_im -> set,pixel_size=[2,2],xyoffset=[-680,-650], $
IDL> det_index_mask=[0,0,1,1,1,1,1,1,1],image_algorithm='clean',clean_niter=100
IDL> im = obj_im -> getdata()
IDL> obj_im -> plotman

(3) 像合成のプログラム

(違う計算方法で像合成する;image reconstruction)
  1. back projection ; まずはざっと放射源をチェックするとき
    (imageのpeakの場所をpoint sourceとして全photon数をばらまく)
    IDL> hessi_image,pars,im,/plot
    するとこんな感じ

  2. clean (back projection のcleaning)
    IDL> hessi_image,im,pars,/plot,image_algorithm='clean'
    するとこんな感じ

    使用するコリメータを変え、細かい構造まで見てみるには、(時間かかる)
    IDL> hessi_image,im,pars,/plot,image_algorithm='clean', det_index_mask=[0,0,1,1,1,1,1,1,1]
    するとこんな感じ

  3. MEM SATO (佐藤さん作成のMEM) ; sharp だけど時間かかります
    なんだかまだうまくいっていないみたい
    (Sato et al. 1999, PASJ 参)
    MEMって? … Maximum Entropy Method
    IDL> hessi_image,im,pars,/plot,image_algorithm='mem sato'
    するとこんな感じ

    (ここで補足)
    小さいピッチのコリメータを使った方が小さい構造まで追えます。 しかし、MEM SATOを計算するにはとっても時間がかかるので、 デフォルトでは番号にして4から8のコリメータ (0〜8の、全部で9つのコリメータがあり、番号が小さい方がピッチが小さい) を使うようにしてあります。
    使用するコリメータを指定するときは、
    IDL> hessi_image,im,pars,/plot,image_algorithm='mem sato', det_index_mask=[0,0,0,0,1,1,1,1,1]
    (この例では4,5,6,7,8のコリメータを使用する)

  4. MEM VIS (MEM SATOより時間かかる)
    まだうまくいっていない
    (light curveではなく、visibility dataでMEMを使用している)
    IDL> hessi_image,im,pars,/plot,image_algorithm='mem vis'
    するとこんな感じ
    MEM SATOのようにdet_index_maskも使えます。

  5. PIXON (更に時間かかる ; マニュアルによると750MHZ Pentium IIIのPCで 3時間半かかったらしい、私はもっとかかったみたい)
    まだうまくいっていない
    (Puetter and Pina, Experimental Astron., 3, 293; Pina and Puetter, 1993, PASP, Metcalf et al., 1996, ApJ 参)
    IDL> hessi_image,im,pars,/plot,image_algorithm='pixon'
    するとこんな感じ

  6. Forward Fitting (結構早くてロバスト)
    たくさんの2D Gaussianの重ね合わせ
    IDL> hessi_image,im,pars,/plot,image_algorithm='forwardfit'
    これはGaussian 1つで計算。
    IDL> hessi_image,im,pars,/plot,image_algorithm='forwardfit', ff_n_gaussians=2
    でGaussianの数を指定する(この場合は2つ)。
    するとこんな感じ

(RHESSI画像の制御キーワード達)

ここには一例を挙げました。 詳しくはhttp://hesperia.gsfc.nasa.gov/ssw/hessi/doc/hsi_params_image_standard.htm を。

(補足)

(4) スペクトル解析

fileの読み込み。
IDL> sp = hsi_spectrum(filename='hsi_20020404_153600_002.fits')
ちなみにパラメタを見るには、
IDL> sp -> print
IDL> sp -> print,/control_only

などとすれば良い。

(指定する必要のあるキーワード達)

以上のことをobjectに放り込みます。
IDL> spec = sp -> getdata()
これでスペクトルを描く準備が出来ました。

(さらに補助的だけど大事なキーワード)

いよいよスペクトルを描きます。
IDL> e = sp -> getdata(/xaxis)
横軸(つまりエネルギー軸)を取り出す。
IDL> plot_oo,e,spec.flux(*,0),ps=3,yr=[1.e-4,1.e+4]
IDL> errplot,e,spec.flux(*,0)-spec.eflux(*,0),spec.flux(*,0)-spec.eflux(*,0), $
width=0.0001

エラーバーを(秘かに)描いているようだ。

*GUIやSPEXでスペクトル解析をするなら http://hesperia.gsfc.nasa.gov/~dennis/spectroscopy/first_steps.htm を参照してください。


Top
Back