pro training_20180328_advanced

; *****************
; settings for the training
; *****************
erg_init, remote_data_dir='https://ergsc.isee.nagoya-u.ac.jp/data/ergsc_training/nagoya_201803/'
uname = '***************'
pass = '****************'

; *****************
; set time span
; *****************
timespan, '2017-03-28', 1, /day

; *****************
; load orbit CDF & set var label
; *****************
set_erg_var_label

; *****************
; load OFA-MATRIX
; *****************
erg_load_pwe_ofa, datatype='matrix', uname=uname, pass=pass
pr_matrix = 'erg_pwe_ofa_matrix_l2_'

; *****************
; load MGF L2 CDF
; *****************
erg_load_mgf, datatype='64hz', coord='sgi', uname=uname, pass=pass
erg_load_mgf, datatype='8sec', uname=uname, pass=pass

; *****************
; analyze OFA-MATRIX
; *****************
get_data, pr_matrix+'Bx_Bx_132', data=s00
get_data, pr_matrix+'Bx_By_132', data=s01
get_data, pr_matrix+'Bx_Bz_132', data=s02

get_data, pr_matrix+'By_Bx_132', data=s10
get_data, pr_matrix+'By_By_132', data=s11
get_data, pr_matrix+'By_Bz_132', data=s12

get_data, pr_matrix+'Bz_Bx_132', data=s20
get_data, pr_matrix+'Bz_By_132', data=s21
get_data, pr_matrix+'Bz_Bz_132', data=s22

rr=dblarr(3,3,n_elements(s00.x),n_elements(s00.v),2)

rr[0,0,*,*,*]=s00.y & rr[1,0,*,*,*]=s01.y & rr[2,0,*,*,*]=s02.y
rr[0,1,*,*,*]=s10.y & rr[1,1,*,*,*]=s11.y & rr[2,1,*,*,*]=s12.y
rr[0,2,*,*,*]=s20.y & rr[1,2,*,*,*]=s21.y & rr[2,2,*,*,*]=s22.y

; *****************
; analyze MGF data
; *****************
split_vec, 'erg_mgf_l2_mag_64hz_sgi'

tinterpol_mxn, 'erg_mgf_l2_mag_64hz_sgi_x', pr_matrix+'Bx_Bx_132'
tinterpol_mxn, 'erg_mgf_l2_mag_64hz_sgi_y', pr_matrix+'Bx_Bx_132'
tinterpol_mxn, 'erg_mgf_l2_mag_64hz_sgi_z', pr_matrix+'Bx_Bx_132'

get_data, 'erg_mgf_l2_mag_64hz_sgi_x_interp', data=data_x, dlim=dlim_x, lim=lim_x
get_data, 'erg_mgf_l2_mag_64hz_sgi_y_interp', data=data_y, dlim=dlim_y, lim=lim_y
get_data, 'erg_mgf_l2_mag_64hz_sgi_z_interp', data=data_z, dlim=dlim_z, lim=lim_z

rotmat=dblarr(3,3,n_elements(data_x.x))
rotmat_t=dblarr(3,3,n_elements(data_x.x))

; rotation matrix
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

; MGF rotation
data_rot=dblarr(n_elements(data_x.x), 3)
for i=0, n_elements(data_x.x)-1 do begin
  data=[[data_x.y[i]], [data_y.y[i]],[data_z.y[i]]]
  data_rot[i,*] = rotmat[*,*,i] ## data
endfor

store_data, 'erg_mgf_l2_mag_64hz_sgi_rot_x', data={x:data_x.x, y:data_rot[*,0]}, dlim=dlim_x, lim=lim_x
store_data, 'erg_mgf_l2_mag_64hz_sgi_rot_y', data={x:data_y.x, y:data_rot[*,1]}, dlim=dlim_y, lim=lim_y
store_data, 'erg_mgf_l2_mag_64hz_sgi_rot_z', data={x:data_z.x, y:data_rot[*,2]}, dlim=dlim_z, lim=lim_z

ylim, 'erg_mgf_l2_mag_64hz_sgi_rot_?', -4E4, 4E4, 0
tplot, [ 'erg_mgf_l2_mag_64hz_sgi_' + ['x', 'y', 'z'], 'erg_mgf_l2_mag_64hz_sgi_rot_' + ['x', 'y', 'z'] ]
stop

; OFA-MATRIX rotation
for i=0, n_elements(s00.x)-1 do begin
  for j=0, n_elements(s00.v)-1 do begin
    for k=0, 1 do begin
      rr[*,*,i,j,k] = (rotmat[*,*,i] ## rr[*,*,i,j,k]) ## rotmat_t[*,*,i]
    endfor
  endfor
endfor

; ************************************
; auto-spectra
; ************************************
get_data, pr_matrix + 'Bx_Bx_132',data=data, dlim=dlim, lim=lim
store_data, pr_matrix + 'Bx_Bx_re', data={x:data.x,y:data.y[*,*,0],v:data.v}, dlim=dlim, lim=lim
store_data, pr_matrix + 'Bx_Bx_im', data={x:data.x,y:data.y[*,*,1],v:data.v}, dlim=dlim, lim=lim
;
get_data, pr_matrix + 'By_By_132',data=data, dlim=dlim, lim=lim
store_data, pr_matrix + 'By_By_re', data={x:data.x,y:data.y[*,*,0],v:data.v}, dlim=dlim, lim=lim
store_data, pr_matrix + 'By_By_im', data={x:data.x,y:data.y[*,*,1],v:data.v}, dlim=dlim, lim=lim
;
get_data, pr_matrix + 'Bz_Bz_132',data=data, dlim=dlim, lim=lim
store_data, pr_matrix + 'Bz_Bz_re', data={x:data.x,y:data.y[*,*,0],v:data.v}, dlim=dlim, lim=lim
store_data, pr_matrix + 'Bz_Bz_im', data={x:data.x,y:data.y[*,*,1],v:data.v}, dlim=dlim, lim=lim

options, pr_matrix + 'Bx_Bx_re', ytitle='Re(BxBx*)' & options, pr_matrix + 'Bx_Bx_im', ytitle='Im(BxBx*)'
options, pr_matrix + 'By_By_re', ytitle='Re(ByBy*)' & options, pr_matrix + 'By_By_im', ytitle='Im(ByBy*)'
options, pr_matrix + 'Bz_Bz_re', ytitle='Re(BzBz*)' & options, pr_matrix + 'Bz_Bz_im', ytitle='Im(BzBz*)'
options, pr_matrix + 'B?_B?_??', ysubtitle = 'Frequency [kHz]'
options, pr_matrix + 'B?_B?_??', ztitle = '[pT!U2!N/Hz]'
ylim, pr_matrix + 'B?_B?_??', 0.064, 20, 1 ; kHz
zlim, pr_matrix + 'B?_B?_??', 1E-3,1E2, 1  ; nT

tplot, pr_matrix + ['Bx_Bx_re', 'Bx_Bx_im', 'By_By_re', 'By_By_im', 'Bz_Bz_re', 'Bz_Bz_im']
stop

; ************************************
; wave normal angle
; ************************************
wna=dblarr(n_elements(s00.x),n_elements(s00.v))

for i=0, n_elements(s00.x)-1 do begin
  for j=0, n_elements(s00.v)-1 do begin
    ; wave normal
    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
endfor

store_data, 'kvec',data={x:s00.x, y:wna, v:s00.v}
options, 'kvec', ytitle='wave normal angle', ztitle='[degree]', ysubtitle='Frequency [kHz]', spec = 1
ylim, 'kvec', 0.064, 20, 1 ; kHz
zlim, 'kvec', 0., 90, 0    ; degree
tplot,[pr_matrix+'Bx_Bx_re', pr_matrix+'By_By_re', pr_matrix+'Bz_Bz_re', 'kvec']
stop

; ************************************
; polarization
; ************************************
polarization=dblarr(n_elements(s00.x),n_elements(s00.v))

for i=0, n_elements(s00.x)-1 do begin
  for j=0, n_elements(s00.v)-1 do begin
    ; polarization
    polarization[i,j] = (2 * rr[1,0,i,j,1]) / (rr[0,0,i,j,0] + rr[1,1,i,j,0])
  endfor
endfor

store_data, 'polarization',data={x:s00.x, y:polarization, v:s00.v}
options, 'polarization', ytitle='polarization', ysubtitle='Frequency [kHz]', spec = 1
ylim, 'polarization', 0.064, 20, 1 ; kHz
zlim, 'polarization', -1., 1., 0   ;
tplot,[pr_matrix+'Bx_Bx_re', pr_matrix+'By_By_re', pr_matrix+'Bz_Bz_re', 'polarization']
stop

; ************************************
; mask
; ************************************
get_data, pr_matrix+'Btotal_132', data=data_ref

; kvec
get_data, 'kvec', data=data, dlim=dlim, lim=lim
data.y[where(data_ref.y LT 1E-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 1E-2)] = 'NaN'
store_data, 'polarization_mask', data={x:data.x, y:data.y, v:data.v}, dlim=dlim, lim=lim

tplot,[pr_matrix+'Bx_Bx_re', pr_matrix+'By_By_re', pr_matrix+'Bz_Bz_re', 'kvec_mask', 'polarization_mask']
stop

; ************************************
; overplot fce
; ************************************
; calculate fce
get_data, 'erg_mgf_l2_magt_8sec', data=data
fce = data.y/ 10^(9.) * 1.6 * 10^(-19.) / (9.1093D * 10^(-31.)) / 2. / !pi / 1000.
fce_half = fce / 2.

store_data, 'fce', data={x:data.x, y:fce}
store_data, 'fce_half', data={x:data.x, y:fce_half}

options, 'fce', colors=fsc_color('yellow'), thick=2
options, 'fce_half', colors=fsc_color('magenta'), thick=2

store_data, pr_matrix+'Btotal_132_gyro', data=[pr_matrix+'Btotal_132', 'fce', 'fce_half']
store_data, 'kvec_mask_gyro', data=['kvec_mask', 'fce', 'fce_half']
store_data, 'polarization_mask_gyro', data=['polarization_mask', 'fce', 'fce_half']

ylim, '*_gyro', 0.064, 20, 1 ; kHz
zlim, pr_matrix+'Btotal_132_gyro', 1E-2, 1E2, 1 ; pT^2/Hz
options, pr_matrix+'Btotal_132_gyro', ytitle='B total', ysubtitle='Frequency [kHz]', ztitle='[pT!U2!N/Hz]'
tplot, [pr_matrix+'Btotal_132', 'kvec_mask', 'polarization_mask']+'_gyro'
stop

; *****************
; load OFA-COMPLEX
; *****************
erg_load_pwe_ofa, datatype='complex', uname=uname, pass=pass
pr_complex = 'erg_pwe_ofa_complex_l2_'

get_data, pr_complex+'Ex_132', data=data_ex, dlim=dlim, lim=lim
c_ex = dcomplex(data_ex.y[*,*,0]*1E-3, data_ex.y[*,*,1]*1E-3)
get_data, pr_complex+'Ey_132', data=data, dlim=dlim, lim=lim
c_ey = dcomplex(data.y[*,*,0]*1E-3, data.y[*,*,1]*1E-3)
get_data, pr_complex+'Bx_132', data=data_bx, dlim=dlim, lim=lim
c_bx = dcomplex(data_bx.y[*,*,0]*1E-12, data_bx.y[*,*,1]*1E-12)
get_data, pr_complex+'By_132', data=data, dlim=dlim, lim=lim
c_by = dcomplex(data.y[*,*,0]*1E-12, data.y[*,*,1]*1E-12)
get_data, pr_complex+'Bz_132', data=data, dlim=dlim, lim=lim
c_bz = dcomplex(data.y[*,*,0]*1E-12, data.y[*,*,1]*1E-12)

; *****************
; calc Ez
; *****************
c_ez = dcindgen(n_elements(data_ex.x), n_elements(data_ex.v))

idx = nn(data_bx.x, data_ex.x)

for i=0, n_elements(data_ex.x)-1 do begin
  c_ez[i,*] = (-c_ex[i,*]*c_bx[idx[i],*]-c_ey[i,*]*c_by[idx[i],*]) / c_bz[idx[i],*]
endfor

Sx=dindgen(n_elements(data_bx.x),n_elements(data_bx.v))
Sy=dindgen(n_elements(data_bx.x),n_elements(data_bx.v))
Sz=dindgen(n_elements(data_bx.x),n_elements(data_bx.v))

; calc Poynting flux
idx = nn(data_ex.x, data_bx.x)
for i=0, n_elements(data_bx.x)-1 do begin
  for j=0, n_elements(data_bx.v)-1 do begin
    Sx[i,j]= double(c_ey[idx[i],j]*conj(c_bz[i,j])-c_ez[idx[i],j]*conj(c_by[i,j]))
    Sy[i,j]=-double(c_ex[idx[i],j]*conj(c_bz[i,j])-c_ez[idx[i],j]*conj(c_bx[i,j]))
    Sz[i,j]= double(c_ex[idx[i],j]*conj(c_by[i,j])-c_ey[idx[i],j]*conj(c_bx[i,j]))
  endfor
endfor


; *****************
; analyze MGF data
; *****************
split_vec, 'erg_mgf_l2_mag_64hz_sgi'

; interpolation
tinterpol_mxn, 'erg_mgf_l2_mag_64hz_sgi_x',pr_complex+'Bx_132',newname='erg_mgf_l2_mag_64hz_sgi_x_interp'
tinterpol_mxn, 'erg_mgf_l2_mag_64hz_sgi_y',pr_complex+'Bx_132',newname='erg_mgf_l2_mag_64hz_sgi_y_interp'
tinterpol_mxn, 'erg_mgf_l2_mag_64hz_sgi_z',pr_complex+'Bx_132',newname='erg_mgf_l2_mag_64hz_sgi_z_interp'

get_data,'erg_mgf_l2_mag_64hz_sgi_x_interp',data=data_x
get_data,'erg_mgf_l2_mag_64hz_sgi_y_interp',data=data_y
get_data,'erg_mgf_l2_mag_64hz_sgi_z_interp',data=data_z

; rotation matrix
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

; *****************
; Poynting vector calculation
; *****************
S=dindgen(n_elements(data_bx.x),n_elements(data_bx.v),3)

for i=0, n_elements(data_bx.x)-1 do begin
  for j=0, n_elements(data_bx.v)-1 do begin
      S[i,j,*] = rotmat[*,*,i] ## [Sx[i,j],Sy[i,j],Sz[i,j]]
  endfor
endfor

theta=acos(S[*,*,2]/sqrt(S[*,*,0]^2+S[*,*,1]^2+S[*,*,2]^2))/!dtor

store_data, 'S', data={x:data_bx.x, y:theta, v:data_bx.v}, dlim=dlim, lim=lim

ylim, [pr_complex+'Etotal_132', pr_complex+'Btotal_132', 'S'], 0.064, 20, 1
zlim, pr_complex+'Etotal_132', 1E-7, 1E0, 1 ; mV^2/m^2/Hz
zlim, pr_complex+'Btotal_132', 1E-2, 1E2, 1 ; pT^2/Hz
zlim, 'S', 0., 180., 0 ; degree

options, [pr_complex+'Etotal_132', pr_complex+'Btotal_132', 'S'], ysubtitle='Frequency [kHz]' 
options, pr_complex+'Etotal_132', ytitle='E total', ztitle='[mV!U2!N/m!U2!N/Hz]'
options, pr_complex+'Btotal_132', ytitle='B total', ztitle='[pT!U2!N/Hz]'
options, 'S', ytitle='Poynting vector', ztitle='[degree]'

tplot, [pr_complex+'Etotal_132', pr_complex+'Btotal_132', 'S']
stop

; ************************************
; mask
; ************************************
get_data, pr_complex+'Btotal_132', data=data_ref

; S
get_data, 'S', data=data, dlim=dlim, lim=lim
data.y[where(data_ref.y LT 1E-2)] = 'NaN'
store_data, 'S_mask', data={x:data.x, y:data.y, v:data.v}, dlim=dlim, lim=lim

tplot, [pr_complex+'Etotal_132', pr_complex+'Btotal_132', 'S_mask']
stop

; ************************************
; overplot fce
; ************************************
store_data, pr_complex+'Etotal_132_gyro', $
data=[pr_complex+'Etotal_132', 'fce', 'fce_half']
store_data, pr_complex+'Btotal_132_gyro', data=[pr_complex+'Btotal_132', 'fce', 'fce_half']
store_data, 'S_mask_gyro', data=['S_mask', 'fce', 'fce_half']

ylim, '*_gyro', 0.064, 20, 1 ; kHz
zlim, pr_complex+'Etotal_132_gyro', 1E-7, 1E0, 1 ; mV^2/m^2/Hz
zlim, pr_complex+'Btotal_132_gyro', 1E-2, 1E2, 1 ; pT^2/Hz
options, pr_complex+'Etotal_132_gyro', ytitle='E total', ysubtitle='Frequency [kHz]', ztitle='[mV!U2!N/m!U2!N/Hz]'
options, pr_complex+'Btotal_132_gyro', ytitle='B total', ysubtitle='Frequency [kHz]', ztitle='[pT!U2!N/Hz]'
tplot, [pr_complex+'Etotal_132', pr_complex+'Btotal_132', 'kvec_mask', 'polarization_mask', 'S_mask'] + '_gyro'
stop

end