;;; This is a crib seet for PWE/EFD spinfit electric field data;;;

;;; initial settings (PPT: p5)
;;; Settings for ERG data analysis (erg_init)

;.full_reset_session
erg_init
;tplot_options, 'charsize', 1.5 ;;enlarge character size, if needed
;time_stamp, /off ;; remove the time stamp in the plot panel
timespan, '2017-12-13'
set_erg_var_label

;;; load the PWE/EFD data () ;;;
erg_load_pwe_efd
erg_load_pwe_efd, datatype='pot'

tplot_names

stop

;;; setting options, and plot (PPT: p7)

options, ['erg_pwe_efd_l2_Eu_*', 'erg_pwe_efd_l2_Ev_*'], datagap=60., labflag=-1

calc, '"erg_pwe_efd_l2_Eu_residual"=abs("erg_pwe_efd_l2_Eu_residual")'
calc, '"erg_pwe_efd_l2_Ev_residual"=abs("erg_pwe_efd_l2_Ev_residual")'

ylim, ['erg_pwe_efd_l2_Eu_residual', 'erg_pwe_efd_l2_Ev_residual'], 0, 30

;;; plot spinfit E-fields and their residual ;;;
tplot, ['erg_pwe_efd_l2_Eu_dsi', 'erg_pwe_efd_l2_Ev_dsi', $
        'erg_pwe_efd_l2_Eu_residual', 'erg_pwe_efd_l2_Ev_residual', $
        'erg_pwe_efd_l2_Vave']

stop

;;; check the quality flag by erg_pwe_qflag (PPT: p18)

erg_pwe_qflag
erg_pwe_qflag, /full

stop

;;; calculate VxB in DSI coordinate (PPT: p21)
;;; load orbit and magnetic field data 

erg_load_orb
erg_load_mgf

tmagname = 'erg_mgf_l2_mag_8sec_dsi'

;;; rotate velocity vector from GSE to DSI via J2000
spd_cotrans, 'erg_orb_l2_vel_gse', 'vel_j2000', in_coord='gse', out_coord='j2000'
erg_cotrans, 'vel_j2000', 'vel_dsi', in_coord='j2000', out_coord='dsi'


;;; calcualte VxB in DSI coordinate
tinterpol_mxn, 'vel_dsi', tmagname
tcrossp, 'vel_dsi_interp', tmagname, newname='VxB_dsi'
calc, '"VxB_dsi"="VxB_dsi"*1e-3' ;;; km/s*nT=10^3*
options, 'VxB_dsi', ysubtitle='[mV/m]', colors=4, labflag=-1

stop

;;; calculate Ez in DSI coordinate & check the ratio of Ez/Exy (PPT: p22)

tinterpol, tmagname, 'erg_pwe_efd_l2_Eu_dsi', newname=tmagname+'_interp'

get_data, 'erg_pwe_efd_l2_Eu_dsi', data=eu
get_data, 'erg_pwe_efd_l2_Ev_dsi', data=ev
get_data, tmagname+'_interp', data=b

;;; for Eu

ez= -(eu.y[*,0]*b.y[*,0]+eu.y[*,1]*b.y[*,1])/b.y[*,2]

eu3=fltarr(n_elements(eu.x),3)
eu3[*,0:1]=eu.y
eu3[*,2]=ez

;;; for Ev

ez= -(ev.y[*,0]*b.y[*,0]+ev.y[*,1]*b.y[*,1])/b.y[*,2]

ev3=fltarr(n_elements(ev.x),3)
ev3[*,0:1]=ev.y
ev3[*,2]=ez

store_data, 'erg_pwe_efd_l2_Eu3_dsi', data={x:eu.x, y:eu3}
options, 'erg_pwe_efd_l2_Eu3_dsi', colors=0
store_data, 'erg_pwe_efd_l2_Ev3_dsi', data={x:ev.x, y:ev3}
options, 'erg_pwe_efd_l2_Ev3_dsi', colors=2

store_data, 'erg_pwe_efd_l2_Eu3_ratio', data={x:eu.x,y:eu3[*,2]/sqrt(eu3[*,0]^2+eu3[*,1]^2)}
store_data, 'erg_pwe_efd_l2_Ev3_ratio', data={x:ev.x,y:ev3[*,2]/sqrt(ev3[*,0]^2+ev3[*,1]^2)}

stop

;;; Comparison of each component of the electric fields
;;; with VxB field. (PPT: p23)

split_vec, 'erg_pwe_efd_l2_Eu3_dsi'
split_vec, 'erg_pwe_efd_l2_Ev3_dsi'
split_vec, 'VxB_dsi'

store_data, 'Ex_with_VxB', data=['erg_pwe_efd_l2_Eu3','erg_pwe_efd_l2_Ev3','VxB']+'_dsi_x'
store_data, 'Ey_with_VxB', data=['erg_pwe_efd_l2_Eu3','erg_pwe_efd_l2_Ev3','VxB']+'_dsi_y'
store_data, 'Ez_with_VxB', data=['erg_pwe_efd_l2_Eu3','erg_pwe_efd_l2_Ev3','VxB']+'_dsi_z'

options, 'E'+['x','y','z']+'_with_VxB', labflag=-1

tplot, 'E'+['x','y','z']+'_with_VxB'

stop

;;; calculate the difference between the E-field and VxB (PPT: p24)

dif_data, 'erg_pwe_efd_l2_Eu3_dsi', 'VxB_dsi', newname='Eu-VxB_dsi'
dif_data, 'erg_pwe_efd_l2_Ev3_dsi', 'VxB_dsi', newname='Ev-VxB_dsi'

;;; rotate the difference from DSI to SM

erg_cotrans, 'Eu-VxB_dsi', 'Eu-VxB_j2000', in_coord='dsi', out_coord='j2000'
erg_cotrans, 'Ev-VxB_dsi', 'Ev-VxB_j2000', in_coord='dsi', out_coord='j2000'

spd_cotrans, 'Eu-VxB_j2000', 'Eu-VxB_sm', in_coord='j2000', out_coord='sm'
spd_cotrans, 'Ev-VxB_j2000', 'Ev-VxB_sm', in_coord='j2000', out_coord='sm'

options, ['Eu-*','Ev-*'], datagap=60., labels=['Ex','Ey','Ez'], labflag=-1
ylim, ['Eu-*', 'Ev-*'], -10, 10

tplot, ['Eu-VxB_sm', 'Ev-VxB_sm']

stop

;;; Calculate the theoretical co-rotation E-field and compare it with the observations
;;; Angular velocity is We=2pi/24hour ~ 7.27*1e-5 rad/s
;;; The E field is calculated in GEO and rotated into SM coordinate system.
;;; (PPT: p25-p26)

spd_cotrans, 'erg_orb_l2_pos_gse', 'pos_geo', in_coord='gse', out_coord='geo'
get_data, 'pos_geo',data=data

xgeo=data.y[*,0]
ygeo=data.y[*,1]
zgeo=data.y[*,2]

nd=n_elements(data.x)
v_cor_geo = fltarr(nd,3)

v_cor_geo_t = sqrt(xgeo^2+ygeo^2)*2.*!pi/86400. ;; Re/s
v_cor_geo[*,0] = -ygeo*v_cor_geo_t/sqrt(xgeo^2+ygeo^2)
v_cor_geo[*,1] = xgeo*v_cor_geo_t/sqrt(xgeo^2+ygeo^2)
v_cor_geo[*,2] = 0.

store_data, 'v_cor_geo', data={x:data.x, y:v_cor_geo}
tkm2re, 'v_cor_geo', /km, /replace
spd_cotrans, 'v_cor_geo', 'v_cor_sm', in_coord='geo', out_coord='sm'

tinterpol, 'erg_mgf_l2_mag_8sec_sm', 'v_cor_sm', newname='erg_mgf_l2_mag_8sec_sm_interp'
tcrossp, 'erg_mgf_l2_mag_8sec_sm_interp', 'v_cor_sm', newname='E_cor_sm'
calc, '"E_cor_sm"="E_cor_sm"*1e-3'
options, 'E_cor_sm', labels=['Ecor_x', 'Ecor_y', 'Ecor_z'], labflag=-1

options, 'Eu-VxB_sm', colors=0
options, 'Ev-VxB_sm', colors=2
options, 'E_cor_sm', colors=4

split_vec, 'Eu-VxB_sm'
split_vec, 'Ev-VxB_sm'
split_vec, 'E_cor_sm'

store_data, 'Ex_comparison', data=['Eu-VxB_sm', 'Ev-VxB_sm', 'E_cor_sm']+'_x'
store_data, 'Ey_comparison', data=['Eu-VxB_sm', 'Ev-VxB_sm', 'E_cor_sm']+'_y'
store_data, 'Ez_comparison', data=['Eu-VxB_sm', 'Ev-VxB_sm', 'E_cor_sm']+'_z'
ylim, ['Ex','Ey','Ez']+'_comparison', -10, 10
options, ['Ex','Ey','Ez']+'_comparison', labflag=-1

tplot, ['Ex','Ey','Ez']+'_comparison'

stop

;;; Potential & PWE/HFA (UHR) comparison: [Kazama et al., 2018] (PPT: p27-p28)

timespan, '2017-04-12'
erg_load_pwe_hfa
erg_load_pwe_ofa
erg_load_pwe_efd, datatype='pot'
set_erg_var_label

tplot, ['erg_pwe_hfa_l2_lm_spectra_e_mix','erg_pwe_ofa_spec_l2_E_spectra_132', 'erg_pwe_efd_l2_Vave']
tlimit, '2017-04-12/'+['16:50', '17:30']

END