thm_init

timespan,'2012-01-22',10,/day

erg_load_gmag_nipr, site=['syo']
erg_load_gmag_mm210,site=['ktb']  
iug_load_gmag_serc,site=['asb']  
kyoto_dst_load

tplot_names
tplot_names, 'magdas_mag_asb',/verbose

window, 0, xsize=600, ysize=700
tplot,['nipr_mag_syo_1sec','magdas_mag_asb','mm210_mag_ktb_1min_hdz','kyoto_dst']

;芦別地磁気を各成分に分解。H成分から平均値を差し引く
split_vec,'magdas_mag_asb'
tsub_average, 'magdas_mag_asb_0'

;コトタバン地磁気を各成分に分解
split_vec,'mm210_mag_ktb_1min_hdz'

;昭和地磁気を各成分に分解。H成分を1分平均値化
split_vec,'nipr_mag_syo_1sec'
avg_data,'nipr_mag_syo_1sec_x',60

tplot,['nipr_mag_syo_1sec_x_avg','magdas_mag_asb_0-d','mm210_mag_ktb_1min_hdz_x','kyoto_dst']

tlimit, '2012-01-22/00:00','2012-01-24/00:00'

;正規表現によるtplot。;芦別地磁気各成分（絶対値）プロット
tplot,'magdas_mag_asb_?'

;芦別地磁気各成分（平均値を引いたもの）を作成
tsub_average, 'magdas_mag_asb_1'
tsub_average, 'magdas_mag_asb_2'
tsub_average, 'magdas_mag_asb_3'

;y軸の幅を-50〜50に変更後プロット
ylim, 'magdas_mag_asb_?-d', -50, 50
tplot,'magdas_mag_asb_?-d'

;時間微分は deriv_data
deriv_data, 'magdas_mag_asb_0'
deriv_data, 'magdas_mag_asb_1'
deriv_data, 'magdas_mag_asb_2'
tplot, 'magdas_mag_asb_?_ddt'

;大量データのウェーブレットは時間がかかるので、必要なデータを切り出す
newname = 'nipr_mag_syo_1sec_x' + '_clip'
t1=time_double('2012-01-22/00:00')
t2=time_double('2012-01-22/12:00')
trange_clip, 'nipr_mag_syo_1sec_x', t1, t2, newname=newname

;昭和地磁気H成分をウェーブレット解析
wav_data,'nipr_mag_syo_1sec_x_clip',/kol ,maxpoints=24l*3600*2

;カラープロットの幅指定
zlim,'*pow', .0001,.1,1
;昭和地磁気H成分とウェーブレット結果の表示
tplot, ['nipr_mag_syo_1sec_x_clip', 'nipr_mag_syo_1sec_x_clip_wv_pow', 'kyoto_dst']

;時間軸を拡大
tlimit, '2012-01-22/00:00','2012-01-22/12:00'

;昭和地磁気H成分にPi 2帯（40s-150s）でバンドパスフィルタ
;tplot変数から値を取り出す
get_data, 'nipr_mag_syo_1sec_x_clip', data=x

;thm_lsp_filterを逐次実行
flow = 1d/150
fhigh =1d/40
dt=1
db=120.0
nyquist = 0.5d/dt
fhigh = double(fhigh/nyquist) < 1.d
flow  = double( flow/nyquist) > 0.d
fmin = min([flow, fhigh])
npts = long(!pi/fmin) > 1
npts = npts < n_elements(x.y)

;デジタルフィルタの定義
cofs = digital_filter(flow,fhigh,db,npts, /double)

;畳み込みとtplot変数再格納
x.y =  convol(x.y,cofs,/edge_t,/nan)
store_data, 'nipr_mag_syo_1sec_x_clip_filt', data=x

;SYOのH成分とバンドパスとウェーブレット表示
ylim, 'nipr_mag_syo_1sec_x_clip_filt', -10, 10
tplot,['nipr_mag_syo_1sec_x_clip','nipr_mag_syo_1sec_x_clip_filt', 'nipr_mag_syo_1sec_x_clip_wv_pow']

;フーリエ変換によるパワースペクトル表示
newname2 = 'nipr_mag_syo_1sec_x_clip' + '_psd'
tdpwrspc, 'nipr_mag_syo_1sec_x_clip', newname=newname2, nboxpoint=512
tplot,['nipr_mag_syo_1sec_x_clip','nipr_mag_syo_1sec_x_clip_filt', 'nipr_mag_syo_1sec_x_clip_wv_pow','nipr_mag_syo_1sec_x_clip_psd']