SPEDAS training session (06 Mar. 2019)
Advanced course [ 2 of 2 ]Arase/ MGF & PWE-EFD
ELF wave analysis
S. Matsuda, S. Kurita, M. Teramoto, M. Shoji, Y. Miyoshi, T.
Hori, T.-F. Chang, Y. Tsugawa, N. Umemura, T. Segawa(ISAS/JAXA & ERG Science Center, ISEE, Nagoya Univ.)

これまでの講習会資料を下記にアップロードしています．ご参考ください．
https://ergsc.isee.nagoya-u.ac.jp/documents/others/data/SPEDAS_training_session/index.shtml.ja

2

PWE/EFD, MGF Dynamic Spectrum Quick Look
https://ergsc.isee.nagoya-u.ac.jp/cef/erg_emic_2hours.cgi

このプロットは，衛星スピン周期よりも速い磁場変動を伴う現象(EMICなど)を 確認することを主な用途とするもので，あくまで Quick Lookプロットとして ご利用いた
だくことを想定しています．

本QLに関してのお問い合わせは
各infoアドレスにお願いいたします．
コンタクト：
erg_pwe_info [at] isee.nagoya-u.ac.jp
erg_mgf_info [at] isee.nagoya-u.ac.jp

3

40. crib sheetを作成する準備
0.1 コンソールからの入力では，複雑な手順を記述することが困難です．
Crib sheet(手順書・スクリプト)を作成しましょう．
2. 新規ファイルを作成
4. プロシージャ名と同じ名前で保存

0.2
以下のコマンドで，cribsheetの途中で「一時停止」「再開」ができます(Crib sheet); *****************; test; *****************
print, 'test1'

training_20190306_adv2

3. Crib sheetの骨格を準備
1. 作業ディレクトリに移動;stopコマンドで一時停止
stop
print, 'test2'(Console)
test1% Stop encountered: *****> .cont
test2

5

1. 必要なファイルの読み込み
1.1 プロットの初期設定と，解析対象の日付を指定します; *****************; initialize; *****************
erg_init
uname = '***************'
pass = '***************'

プロジェクトメンバー限定でファイルを公開しています．
ユーザー名/パスワードの取り扱いには十分ご注意ください!!!

; *****************; set time span
; *****************
timespan, '2017-12-11T13:00:00', 2, /hour

1.2 軌道L2 CDFファイルを取得して，プロット下部に表示するように設定します; *****************; load orbit CDF & set var label; *****************
set_erg_var_label

6

1. 必要なファイルの読み込み
1.3 MGF 64 Hzデータ(DSI座標系) のCDFを取得します'8sec''64hz''256hz': 8秒平均データ: 64 Hz データ: 256 Hzデータ; *****************; load MGF 64 Hz CDF; *****************
erg_load_mgf, datatype='64hz', coord='dsi', uname=uname, pass=passtplot, 'erg_mgf_l2_mag_64hz_dsi'> tplot_names･･･
27 erg_mgf_l2_date_time_64hz
28 erg_mgf_l2_mag_64hz_dsi
29 erg_mgf_l2_dyn_rng_64hz
30 erg_mgf_l2_quality_64hz
31 erg_mgf_l2_quality_64hz_gc
32 erg_mgf_l2_spin_phase_64hz
33 erg_mgf_l2_ti_64hz'sgi': Spinning satellite Geometry Inertia coordinate'dsi': Despun Sun sector Inertia coordinate

7

1. パワースペクトル解析
1.4 tdpwrspcコマンドを用いて，各軸のパワースペクトルを計算します; *****************; FFT analyses (MGF); *****************
nbox = 2048
nshift = 256
split_vec, 'erg_mgf_l2_mag_64hz_dsi'

ERG> get_data,'erg_mgf_l2_mag_64hz_dsi_x_dpwrspc', data=data
ERG> help, data** Structure <2479b5f0>, 3 tags, length=14710800, data
length=14710800, refs=1:
X
DOUBLE Array[1794]
(時刻)
Y
FLOAT
Array[1794, 1024] (パワースペクトル)V
FLOAT
Array[1794, 1024] (周波数)tdpwrspc, 'erg_mgf_l2_mag_64hz_dsi_x', nbox=nbox, nshift=nshift, bin=1
tdpwrspc, 'erg_mgf_l2_mag_64hz_dsi_y', nbox=nbox, nshift=nshift, bin=1
tdpwrspc, 'erg_mgf_l2_mag_64hz_dsi_z', nbox=nbox, nshift=nshift, bin=1
Tips: 時系列波形のShort-time Fourier transformを計算します -> (元のtplot変数名)_dpwrspc
tdpwrspc, 時系列波形のtplot変数名, nbox=FFT点数, nshift=overlap点数, bin=周波数方向の平均化点数

1.5 軸やラベルの設定を行います;y-title
options, 'erg_mgf_l2_mag_64hz_dsi_x_dpwrspc', ytitle='Bx spectra'
options, 'erg_mgf_l2_mag_64hz_dsi_y_dpwrspc', ytitle='By spectra'
options, 'erg_mgf_l2_mag_64hz_dsi_z_dpwrspc', ytitle='Bz spectra';y-subtitle
options, 'erg_mgf_l2_mag_64hz_dsi_?_dpwrspc', ysubtitle='Frequency [Hz]';z-title
options, 'erg_mgf_l2_mag_64hz_dsi_?_dpwrspc', ztitle='nT!U2!N/Hz';
ylim, 'erg_mgf_l2_mag_64hz_dsi_?_dpwrspc', 0, 2, 0zlim, 'erg_mgf_l2_mag_64hz_dsi_?_dpwrspc', 1E-4, 1E+2, 1;
tplot, 'erg_mgf_l2_mag_64hz_dsi_' + ['x','y','z'] + '_dpwrspc'

8

1. パワースペクトル解析
1.6 計算したパワースペクトルを描画します
tplot, 'erg_mgf_l2_mag_64hz_dsi_' + ['x','y','z'] + '_dpwrspc'

9

1. パワースペクトル解析
1.7

MGF 8secデータからイオンサイクロトロン周波数を計算します'8sec''64hz''256hz': 8秒平均データ: 64 Hz データ: 256 Hzデータ; *****************; load MGF L2 CDF; *****************
erg_load_mgf, datatype='8sec', uname=uname, pass=pass; ************************************; overplot fce (EFD & MGF)
; ************************************
get_data, 'erg_mgf_l2_magt_8sec', data=data
fce = data.y / 10^(9.) * 1.6 * 10^(-19.) $/ (9.1093D * 10^(-31.)) / 2. / !pi ; Hz
fch = fce / 1837.
;
store_data, 'fch', data={x:data.x, y:fch}
store_data, 'fche', data={x:data.x, y:fch/4.}
store_data, 'fco', data={x:data.x, y:fch/16.};
options, 'fch', colors=fsc_color('yellow'), thick=2
options, 'fche', colors=fsc_color('magenta'), thick=2
options, 'fco', colors=fsc_color('red'), thick=2;
tplot, ['fch', 'fche', 'fco']

10

1. パワースペクトル解析
1.8

イオンサイクロトロン周波数を重ねてプロットします

Tips: 複数のtplot変数をマージしたtplot変数を作る (store_data)store_data, tplot変数1, data=[tplot変数2, tplot変数3]
tplot変数2とtplot変数3をマージして，tplot変数1に格納する．
といっても，tplot変数1の実体は 'tplot変数2, tplot変数3' という単なる文字列
store_data, 'erg_mgf_l2_mag_64hz_dsi_x_dpwrspc_gyro',data=['erg_mgf_l2_mag_64hz_dsi_x_dpwrspc', 'fch', 'fche', 'fco'];
store_data, 'erg_mgf_l2_mag_64hz_dsi_y_dpwrspc_gyro',data=['erg_mgf_l2_mag_64hz_dsi_y_dpwrspc', 'fch', 'fche', 'fco'];
store_data, 'erg_mgf_l2_mag_64hz_dsi_z_dpwrspc_gyro',data=['erg_mgf_l2_mag_64hz_dsi_z_dpwrspc', 'fch', 'fche', 'fco'];ylim, 'erg_mgf_l2_mag_64hz_dsi_?_dpwrspc_gyro', 0, 2, 0;
tplot, 'erg_mgf_l2_mag_64hz_dsi_' + ['x','y','z'] + '_dpwrspc_gyro'(DSI)

1. パワースペクトル解析
Tips:
スピン周波数近傍の波の解析は，非スピン座標系 (あらせ衛星では DSI座標系 )で行うこと
𝑉𝑠 cos 𝜔signal 𝑡 cos 𝜔spin 𝑡
𝑉𝑠=cos(𝜔spin 𝑡 + 𝜔signal 𝑡) + cos(𝜔spin 𝑡 − 𝜔signal 𝑡)2; *****************; load MGF 64 Hz CDF (SGI); *****************
erg_load_mgf, datatype='64hz', coord='sgi', uname=uname, pass=passtplot, 'erg_mgf_l2_mag_64hz_sgi'; *****************; FFT analyses (MGF, SGI)
; *****************
nbox = 2048
nshift = 256
split_vec, 'erg_mgf_l2_mag_64hz_sgi'
tdpwrspc, 'erg_mgf_l2_mag_64hz_sgi_x', nbox=nbox, nshift=nshift, bin=1
tdpwrspc, 'erg_mgf_l2_mag_64hz_sgi_y', nbox=nbox, nshift=nshift, bin=1
tdpwrspc, 'erg_mgf_l2_mag_64hz_sgi_z', nbox=nbox, nshift=nshift, bin=1;
options, 'erg_mgf_l2_mag_64hz_sgi_x_dpwrspc', ytitle='Bx spectra'
options, 'erg_mgf_l2_mag_64hz_sgi_y_dpwrspc', ytitle='By spectra'
options, 'erg_mgf_l2_mag_64hz_sgi_z_dpwrspc', ytitle='Bz spectra'
options, 'erg_mgf_l2_mag_64hz_sgi_?_dpwrspc', ysubtitle='Frequency [Hz]'options, 'erg_mgf_l2_mag_64hz_sgi_?_dpwrspc', ztitle='nT!U2!N/Hz';
ylim, 'erg_mgf_l2_mag_64hz_sgi_?_dpwrspc', 0, 2, 0zlim, 'erg_mgf_l2_mag_64hz_sgi_?_dpwrspc', 1E-4, 1E+2, 1
tplot, 'erg_mgf_l2_mag_64hz_sgi_' + ['x','y','z'] + '_dpwrspc'

(SGI)

11

12

2. 時間波形からスペクトルマトリクスを作る
2.1 MGF 64Hzデータ(DSI座標系)から，スペクトルマトリクスの要素を計算します; *****************; calculate spectral matrix elements; *****************
get_data, 'erg_mgf_l2_mag_64hz_dsi', data=mgf_64hz
ndata = n_elements(mgf_64hz.x)
ERG> help, 'mgf_64hz_fft'
MGF_64HZ_FFT DCOMPLEX = Array[1794, 2048, 3]

mgf_64hz_fft = dcomplexarr((ndata-nbox)/nshift+1, nbox, 3)
mgf_64hz_t = dblarr((ndata-nbox)/nshift+1)
w = hanning(nbox, /double)
fs = 1.0 / (mgf_64hz.x[1] - mgf_64hz.x[0])ERG> help,w
W
DOUBLE = Array[2048]
ERG> plot, indgen(n_elements(w)), w
ERG> help,fs
FS
DOUBLE =64.103683

13

2. 時間波形からスペクトルマトリクスを作る
2.1 MGF 64Hzデータ(DSI座標系)から，スペクトルマトリクスの要素を計算します(続き)
ERG> help,mgf_64hz_fft[0, 0, 0]<Expression> DCOMPLEX = (43.351752,

0.00000000)

ERG> help,mgf_64hz_fft[0, 1, 0]<Expression> DCOMPLEX = ( -21.687504,

0.015383164)

i = 0Lfor j=0L, ndata-nbox-1, nshift do begin
for k=0, 2 do begin
mgf_64hz_fft[i,*,k] = fft(mgf_64hz.y[j:j+nbox-1,k] * w)
endformgf_64hz_t[i++] = mgf_64hz.x[j]endfornshiftnbox->FFT[j][ｊ+nbox-1]nbox->FFT

i=inshifti = i+1[j]nshifti = i+2[ｊ+nbox-1]nbox->FFT14

2. 時間波形からスペクトルマトリクスを作る
2.1 MGF 64Hzデータ(DSI座標系)から，スペクトルマトリクスの要素を計算します(続き)𝑁

𝑆𝐵 𝑓 =

𝑛|𝐵u |2

1
෍ 𝐵v(𝑛) 𝐵u𝑛 ∗𝑁(𝑛) 𝑛 ∗
𝑛=0
𝐵w 𝐵u

ERG> help, t_s
T_S
DOUBLE = Array[1794]

(𝑛)

𝑛 ∗

(𝑛)

𝑛 ∗

(𝑛)

𝑛 ∗

𝐵u 𝐵v

𝐵u 𝐵w

𝑛

𝐵v 𝐵w

|𝐵v |2
(𝑛)

𝑛 ∗

𝐵w 𝐵v

𝑛

|𝐵w |2npt = (ndata-nbox)/nshift + 1
t_s = mgf_64hz_t[0:n_elements(mgf_64hz_t)-1]

ERG> help, freq
FREQ
DOUBLE = Array[1024]

freq = findgen(nbox/2) * fs / nbox
w_ave = average(w^2)

ERG> help, w_ave<Expression>
DOUBLE = 0.37500000

ERG> help, bubu1
BUBU1
DCOMPLEX = Array[1794, 2048]

bubu1 = mgf_64hz_fft[*,*,0] * conj(mgf_64hz_fft[*,*,0]) / w_ave
bubv1 = mgf_64hz_fft[*,*,0] * conj(mgf_64hz_fft[*,*,1]) / w_ave
bubw1 = mgf_64hz_fft[*,*,0] * conj(mgf_64hz_fft[*,*,2]) / w_ave
bvbu1 = mgf_64hz_fft[*,*,1] * conj(mgf_64hz_fft[*,*,0]) / w_ave
bvbv1 = mgf_64hz_fft[*,*,1] * conj(mgf_64hz_fft[*,*,1]) / w_ave
bvbw1 = mgf_64hz_fft[*,*,1] * conj(mgf_64hz_fft[*,*,2]) / w_ave
bwbu1 = mgf_64hz_fft[*,*,2] * conj(mgf_64hz_fft[*,*,0]) / w_ave
bwbv1 = mgf_64hz_fft[*,*,2] * conj(mgf_64hz_fft[*,*,1]) / w_ave
bwbw1 = mgf_64hz_fft[*,*,2] * conj(mgf_64hz_fft[*,*,2]) / w_ave

15

2. 時間波形からスペクトルマトリクスを作る
2.2 時間平均されたスペクトルマトリクスの要素を計算します; *****************; time average of spectral matrix elements; *****************
bubu2 = dcomplexarr(npt, n_elements(freq))bubv2 = dcomplexarr(npt, n_elements(freq))bubw2 = dcomplexarr(npt, n_elements(freq))
bvbu2 = dcomplexarr(npt, n_elements(freq))bvbv2 = dcomplexarr(npt, n_elements(freq))
bvbw2 = dcomplexarr(npt, n_elements(freq))
bwbu2 = dcomplexarr(npt, n_elements(freq))
bwbv2 = dcomplexarr(npt, n_elements(freq))
bwbw2 = dcomplexarr(npt, n_elements(freq))

𝑁

𝑆𝐵 𝑓 =

𝑛|𝐵u |2

1
෍ 𝐵v(𝑛) 𝐵u𝑛 ∗𝑁(𝑛) 𝑛 ∗
𝑛=0
𝐵w 𝐵u

(𝑛)

𝐵u 𝐵v

𝑛 ∗

𝑛

|𝐵v |2
(𝑛)

𝐵w 𝐵v

𝑛 ∗

for i=0, npt-1 do begin
indx = i + indgen(3)if( (where(indx LT 0))[0] NE -1 ) then breakif( (where(indx GT npt-1))[0] NE -1 ) then breakfor j=0, n_elements(freq)-1 do begin
bubu2[i, j] = average(bubu1[indx, j])bubv2[i, j] = average(bubv1[indx, j])bubw2[i, j] = average(bubw1[indx, j])bvbu2[i, j] = average(bvbu1[indx, j])bvbv2[i, j] = average(bvbv1[indx, j])bvbw2[i, j] = average(bvbw1[indx, j])bwbu2[i, j] = average(bwbu1[indx, j])bwbv2[i, j] = average(bwbv1[indx, j])bwbw2[i, j] = average(bwbw1[indx, j])endfor
endfor

ERG> help, bubu2
BUBU2
DCOMPLEX = Array[1794, 1024]

(𝑛)

𝑛 ∗

(𝑛)

𝑛 ∗

𝐵u 𝐵w
𝐵v 𝐵w
𝑛

|𝐵w |216

2. 時間波形からスペクトルマトリクスを作る
2.3 各要素から，スペクトルマトリクスを構成します𝑁

𝑆𝐵 𝑓 =

𝑛|𝐵u |2

1
෍ 𝐵v(𝑛) 𝐵u𝑛 ∗𝑁(𝑛) 𝑛 ∗
𝑛=0
𝐵w 𝐵u

(𝑛)

𝐵u 𝐵v

𝑛 ∗

𝑛

|𝐵v |2
(𝑛)

𝐵w 𝐵v

𝑛 ∗

(𝑛)

𝑛 ∗

(𝑛)

𝑛 ∗

𝐵u 𝐵w
𝐵v 𝐵w
𝑛

|𝐵w |2

ERG> help, rr
RR
DOUBLE = Array[3, 3, 1794, 1024, 2]
3:
3:
1794:
1024:
2:

列方向のindex
行方向のindex
時刻のindex
周波数のindex
実部/虚部のindex; *****************; make spectral matrix; *****************
rr=dblarr(3, 3, npt, n_elements(freq), 2)
rr[0,0,*,*,0]=real_part(bubu2)
rr[1,0,*,*,0]=real_part(bubv2)
rr[2,0,*,*,0]=real_part(bubw2)
rr[0,1,*,*,0]=real_part(bvbu2)
rr[1,1,*,*,0]=real_part(bvbv2)
rr[2,1,*,*,0]=real_part(bvbw2)
rr[0,2,*,*,0]=real_part(bwbu2)
rr[1,2,*,*,0]=real_part(bwbv2)
rr[2,2,*,*,0]=real_part(bwbw2)
rr[0,0,*,*,1]=imaginary(bubu2)
rr[1,0,*,*,1]=imaginary(bubv2)
rr[2,0,*,*,1]=imaginary(bubw2)
rr[0,1,*,*,1]=imaginary(bvbu2)
rr[1,1,*,*,1]=imaginary(bvbv2)
rr[2,1,*,*,1]=imaginary(bvbw2)
rr[0,2,*,*,1]=imaginary(bwbu2)
rr[1,2,*,*,1]=imaginary(bwbv2)
rr[2,2,*,*,1]=imaginary(bwbw2)

17

3. 磁場座標系への回転行列
3.1 変動成分を取り除くために，MGF 64Hzデータ(DSI座標系)の移動平均を計算します
Tips: データの移動平均を計算する (tsmooth_in_time)
tsmooth_in_time, tplot変数1, 平均幅, newname=tplot変数2
tplot変数1の時系列データに対してrunning averageを計算し、結果をtplot変数2に格納する．
; *****************
Running averageの結果は、平均する前のtplot変数と同じ時刻タグを付けて返される; MGF matrix rotation; *****************
tsmooth_in_time, 'erg_mgf_l2_mag_64hz_dsi', 60., newname='erg_mgf_l2_mag_64hz_dsi_smooth'
tplot, ['erg_mgf_l2_mag_64hz_dsi', 'erg_mgf_l2_mag_64hz_dsi_smooth']18

3. 磁場座標系への回転行列
3.2 スムージングされたMGF 64Hzデータ(DSI座標系)を，スペクトルマトリクスの時刻タグに合わせて内挿します
split_vec, 'erg_mgf_l2_mag_64hz_dsi_smooth'tinterpol_mxn, 'erg_mgf_l2_mag_64hz_dsi_smooth_x', t_s, $newname='erg_mgf_l2_mag_64hz_dsi_smooth_x_interp'
tinterpol_mxn, 'erg_mgf_l2_mag_64hz_dsi_smooth_y', t_s, $newname='erg_mgf_l2_mag_64hz_dsi_smooth_y_interp'
tinterpol_mxn, 'erg_mgf_l2_mag_64hz_dsi_smooth_z', t_s, $newname='erg_mgf_l2_mag_64hz_dsi_smooth_z_interp'
Tips: データを内挿する (tinterpol_mxn)
tinterpol_mxn, tplot変数1, tplot変数2, newname=tplot変数3
tplot変数2の時刻タグに合わせて，tplot変数1のデータを内挿し，tplot変数3に格納する
newnameを省略した場合は，'元のtplot変数名_interp'という名前の変数に格納される
get_data, 'erg_mgf_l2_mag_64hz_dsi_smooth_x_interp', data = data_x
get_data, 'erg_mgf_l2_mag_64hz_dsi_smooth_y_interp', data = data_y
get_data, 'erg_mgf_l2_mag_64hz_dsi_smooth_z_interp', data = data_z193. 磁場座標系への回転行列
3.3 内挿したスムージング済みMGF 64Hzデータ（DSI座標系）から，
DSI座標系 → 背景磁場座標系(B0 // z) の変換行列を求めます
rotmat=dblarr(3,3,n_elements(data_x.x))
rotmat_t=dblarr(3,3,n_elements(data_x.x))
for i=0, n_elements(data_x.x)-1 do begin
bvec=[data_x.y[i],data_y.y[i],data_z.y[i]]zz=[0.,0.,1.]
yhat=crossp(zz,bvec)
xhat=crossp(yhat,bvec)
zhat=bvec
xhat=xhat/sqrt(xhat[0]^2+xhat[1]^2+xhat[2]^2)
yhat=yhat/sqrt(yhat[0]^2+yhat[1]^2+yhat[2]^2)
zhat=zhat/sqrt(zhat[0]^2+zhat[1]^2+zhat[2]^2)
rotmat[*,*,i]=([[xhat],[yhat],[zhat]])rotmat_t[*,*,i]=transpose([[xhat],[yhat],[zhat]])endfor20

3. 磁場座標系への回転
3.4 求めた回転行列を用いて，スペクトルマトリクス (DSI座標系)を背景磁場座標系(B0 // z)に変換します; matrix rotation
for i=0, n_elements(t_s)-1 do begin

ERG> help, rotmat[0,*,*]<Expression> DOUBLE = Array[1, 3, 3]

for j=0, n_elements(freq)-1 do begin
for k=0, 1 do begin

ERG> help, reform(rotmat[0,*,*])<Expression> DOUBLE = Array[3, 3]

rr[*,*,i,j,k] = $(reform(rotmat[i,*,*]) ## rr[*,*,i,j,k]) ## reform(rotmat_t[i,*,*])
endfor
endfor
endforTips: 行列の積 (A # B, A ## B)

分散共分散行列の定義
𝑪 = 𝑬 𝑿𝑿𝑇
回転行列を𝑹𝑿′ = 𝑹𝑿 とする
𝑇

C = A # B → C = BA
C = A ## B → C = AB

𝑪 ′ = 𝑬 𝑿′ 𝑿′= 𝑬 𝑹𝑿𝑿𝑇 𝑹𝑇= 𝑹𝑬 𝑿𝑿𝑇 𝑹𝑇= 𝑹𝑪𝑹𝑇

21

4. パワースペクトルをプロット
4.1 計算したスペクトルマトリクスから，パワースペクトルをプロットしますスペクトルマトリクスの対角成分は，各軸のパワースペクトルを表します（対角成分のimaginary partは常にゼロ）< 𝐵𝑥 𝐵𝑥 ∗ >∗S = < 𝐵𝑦 𝐵𝑥 >< 𝐵𝑧 𝐵𝑥 ∗ >< 𝐵𝑥 𝐵𝑦 ∗ >< 𝐵𝑦 𝐵𝑦 ∗ >< 𝐵𝑧 𝐵𝑦 ∗ >|𝐵𝑥 |2= |𝐵𝑦 ||𝐵𝑥 |𝑒 −𝑖 𝜙𝑥 −𝜙𝑦|𝐵𝑧 ||𝐵𝑥 |𝑒 −𝑖 𝜙𝑥 −𝜙𝑧< 𝐵𝑥 𝐵𝑧 ∗ >< 𝐵𝑦 𝐵𝑧 ∗ >< 𝐵𝑧 𝐵𝑧 ∗ >|𝐵𝑥 ||𝐵𝑦 |𝑒 𝑖 𝜙𝑥 −𝜙𝑦|𝐵𝑥 ||𝐵𝑧 |𝑒 𝑖 𝜙𝑥 −𝜙𝑧|𝐵𝑦 |2|𝐵𝑦 ||𝐵𝑧 |𝑒 𝑖 𝜙𝑦 −𝜙𝑧|𝐵𝑧 ||𝐵𝑦 |𝑒 −𝑖 𝜙𝑦 −𝜙𝑧|𝐵𝑧 |2

ERG> help, rr[0, 0, *, *, 0]<Expression> DOUBLE = Array[1, 1, 1794, 1024]
ERG> help, reform( rr[0, 0, *, *, 0] )<Expression> DOUBLE = Array[1794, 1024]; ************************************; auto-spectra; ************************************
store_data, 'Bx_Bx_re', data={x:t_s, y:reform(rr[0,0,*,*,0])/(fs/nbox), v:freq}
store_data, 'By_By_re', data={x:t_s, y:reform(rr[1,1,*,*,0])/(fs/nbox), v:freq}
store_data, 'Bz_Bz_re', data={x:t_s, y:reform(rr[2,2,*,*,0])/(fs/nbox), v:freq}
calc, '"Btotal"="Bx_Bx_re"+"By_By_re"+"Bz_Bz_re"'options, ['Bx_Bx_re','By_By_re','Bz_Bz_re','Btotal'], ysubtitle = 'Frequency [Hz]'options, ['Bx_Bx_re','By_By_re','Bz_Bz_re','Btotal'], ztitle = '[nT!U2!N/Hz]'
options, ['Bx_Bx_re','By_By_re','Bz_Bz_re','Btotal'], spec=1
ylim, ['Bx_Bx_re','By_By_re','Bz_Bz_re','Btotal'], 0, 2, 0 ; Hz
zlim, ['Bx_Bx_re','By_By_re','Bz_Bz_re','Btotal'], 1E-4, 1E2, 1 ; nT^2/Hz;
tplot, ['Bx_Bx_re','By_By_re','Bz_Bz_re','Btotal']

22

4. パワースペクトルをプロット
4.1 計算したスペクトルマトリクスから，パワースペクトルをプロットします23

5. Means’ methodによるwave normal angle推定
5.1 スペクトルマトリクス (背景磁場座標系)にMeans’ methodを適用します; ************************************; wave normal angle; ************************************
wna = dblarr(n_elements(t_s), n_elements(freq))
for i=0, n_elements(t_s)-1 do begin
for j=0, n_elements(freq)-1 do begin
wna_x = rr[2,1,i,j,1] / sqrt(rr[2,1,i,j,1]^2 + rr[0,2,i,j,1]^2 + rr[1,0,i,j,1]^2)
wna_y = rr[0,2,i,j,1] / sqrt(rr[2,1,i,j,1]^2 + rr[0,2,i,j,1]^2 + rr[1,0,i,j,1]^2)
wna_z = rr[1,0,i,j,1] / sqrt(rr[2,1,i,j,1]^2 + rr[0,2,i,j,1]^2 + rr[1,0,i,j,1]^2)
wna[i,j] = abs(atan(sqrt(wna_x^2+wna_y^2)/wna_z)/!dtor)
endfor
endforMeans et al. (1972)𝑘 ∥Im(𝑆YZ )Im(𝑆ZX )Im(𝑆𝑋𝑌 )store_data, 'kvec',data={x:t_s, y:wna, v:freq}
options, 'kvec', ytitle='wave normal angle', ztitle='[degree]', ysubtitle='Frequency[Hz]', spec = 1
ylim, 'kvec', 0, 2, 0 ; Hz
zlim, 'kvec', 0, 90, 0 ; degree

24

6. Polarizationを計算
6.1

磁力線に対して垂直な面における信号を，次のようにおく
𝐻𝑥 𝑡 = 𝑎𝑥 𝑡 𝑒 𝑖 𝑤𝑡+𝜙𝑥 𝑡
𝐻𝑦 𝑡 = 𝑎𝑦 (𝑡)𝑒 𝑖(𝑤𝑡+𝜙𝑦 𝑡 )スペクトルマトリクス(背景磁場座標系)から， スペクトルマトリクスは< 𝐻𝑥 𝐻𝑥 ∗ > < 𝐻𝑥 𝐻𝑦 ∗ >偏波(polarization)を計算しますS=∗∗< 𝐻𝑦 𝐻𝑥 > < 𝐻𝑦 𝐻𝑦 >; ************************************; polarization; ************************************
polarization = $dblarr(n_elements(t_s), n_elements(freq))=< 𝑎𝑥 2 >< 𝑎𝑥 𝑎𝑦 𝑒 𝑖 𝜙𝑥 −𝜙𝑦 >< 𝑎𝑥 𝑎𝑦 𝑒 −𝑖 𝜙𝑥 −𝜙𝑦< 𝑎𝑦 2 >𝑎𝑥 𝑎𝑦 𝑒 𝑖 𝜙𝑥−𝜙𝑦 = 𝑎𝑥 𝑎𝑦 cos 𝜙𝑥 − 𝜙𝑦 + 𝑖𝑎𝑥 𝑎𝑦 sin 𝜙𝑥 − 𝜙𝑦
偏波率ηは
< 2𝑎𝑥 𝑎𝑦 sin 𝜙𝑥 − 𝜙𝑦 > 2Im(𝑠10 )𝜂==< 𝑎𝑥 2 > + < 𝑎𝑦 2 >𝑠00 + 𝑠11

for i=0, n_elements(t_s)-1 do begin
for j=0, n_elements(freq)-1 do begin
polarization[i,j] = (2 * rr[1,0,i,j,1]) / (rr[0,0,i,j,0] + rr[1,1,i,j,0])
endfor
endforstore_data, 'polarization',data={x:t_s, y:polarization, v:freq}
options, 'polarization', ytitle='polarization', $ysubtitle='Frequency [Hz]', spec = 1
ylim, 'polarization', 0, 2, 0 ; Hz
zlim, 'polarization', -1, 1, 0 ;25

7. 解析結果をマスク
7.1

パワースペクトル(Btotal)を使って，解析結果をマスクします; ************************************; mask; ************************************
get_data, 'Btotal', data=data_ref

Tips: 条件を満たす配列要素の番号を返す (where)
ret = where( 配列 <演算子> 条件 )配列の中から，引数として与えられる条件式を満たす要素の番号を，配列として返す; kvecget_data, 'kvec', data=data, dlim=dlim, lim=lim
data.y[where(data_ref.y LT 5E-2)] = 'NaN'
store_data, 'kvec_mask', data={x:data.x, y:data.y, v:data.v}, dlim=dlim, lim=lim; polarization
get_data, 'polarization', data=data, dlim=dlim, lim=lim
data.y[where(data_ref.y LT 5E-2)] = 'NaN'
store_data, 'polarization_mask', data={x:data.x, y:data.y, v:data.v}, dlim=dlim, lim=lim
tplot, ['Btotal', 'kvec_mask', 'polarization_mask']

26

8. まとめてプロット
8.1

MGF 8秒平均値から電子サイクロトロン周波数を計算して，重ねてプロットします; ************************************; overplot fce; ************************************
store_data, 'Btotal_gyro', data=['Btotal', 'fch', 'fche', 'fco']
store_data, 'kvec_mask_gyro', data=['kvec_mask', 'fch', 'fche', 'fco']
store_data, 'polarization_mask_gyro', $data=['polarization_mask', 'fch', 'fche', 'fco'];ylim, 'Btotal_gyro', 0, 2, 0ylim, 'kvec_mask_gyro', 0, 2, 0ylim, 'polarization_mask_gyro', 0, 2, 0
tplot, ['Btotal', 'kvec_mask', 'polarization_mask'] + '_gyro'27

Memo