;<<<SOTデータ読み込み ＋ 黒点tracking メモ>>>
;前半はSOT(FG)のデータ読み込みに最低限必要な作業です。

; 一例として2007/05/02 の FG のUT 02:00 台のデータファイル名一覧を取得
file=find_file('/solardb/hinode/sot/level0/2007/05/02/FG/H0200/*.fits')

; データ読み込み
read_sot, file,index,data
;;;
;;;index にヘッダ情報、 data にデータキューブが入る。

; キャリブレーション
fg_prep,index,data,index_out,data_out
;;;
;;;注意１：1ファイルにつき数秒かかります。
;;;注意２：NFIのデータを使うときは /no_flatをつけてください。


; G band のデータ番号をとりだす。
g=where(index.wave eq 'G band 4305')

; G band のデータだけをとりだした配列を作成。
index_G=index_out[g]
data_G=data_out[*,*,g]

; Ca II H のデータ番号をとりだす。
ca=where(index.wave  eq 'Ca II H line')

; Ca II H のデータだけをとりだした配列を作成。
index_Ca=index_out[ca]
data_Ca=data_out[*,*,ca]


;;;;以下は黒点trackingのための作業の一例です。

; 黒点を追いかけたいとき：correlation trackingを行う
; (通常はgranuleパターンをもとに追尾されている)
Goffset=get_correl_offsets(data_G)
; Goffsetは「一枚目に対する」ずれの配列。
; Goffset[0,*] x座標
; Goffset[1,*] y座標

; これをこのまま使うとガタガタするので、
; polynomial fittingするとよい。
polyoffx=poly_fit(indgen(60),Goffset[0,*],2)
polyoffy=poly_fit(indgen(60),Goffset[0,*],2)
; 最後の2は次数。この場合２次曲線fitting。

;fitting後の配列を用意。
p_offset=Goffset

; fittingしたずれの座標の関数を p_offset[0,*]にx座標、[1,*]にy座標で格納。
for i=0, 59 do begin
 p_offset[0,i]=float(i)^2*polyoffx[0,2]+float(i)*polyoffx[0,1]+polyoffx[0,0]
 p_offset[1,i]=float(i)^2*polyoffy[0,2]+float(i)*polyoffy[0,1]+polyoffy[0,0]
endfor

;このfittingしたずれを使って、データ自体をshiftさせて、
; tracking済みのデータキューブを p_Gdataに格納。
p_Gdata=shift_img(data_G,p_offset)


;Ca II H のデータも同様に。

以上