; +
; Joint ESA + SST moments calculations
; Version 01
; The code is partly written by Pat Cruce.
; Level 1 SST and ESA data are used
; ==================================================================

date='2008-03-01'
startdate = '2008-03-01/00:00'

timespan,startdate, 4.0, /hour

Trange=['08-03-01/00:00','08-03-01/04:00']
Tzoom=['08-03-01/01:40','08-03-01/02:40']

;... select exact time interval to calculate join ESA/SST moments
tbeg = time_double(date+'/00:00')
tend = time_double(date+'/04:00')

;select a probe
sc='b'

cols=get_colors()

;...load all data over the requested interval

thm_load_state,probe=sc,coord='gsm',/get_support
thm_load_fit, level=1, probe=sc, datatype=['efs', 'fgs'],/verbose
thm_cotrans,strjoin('th'+sc+'_fgs'),out_suf='_gsm', in_c='dsl', out_c='gsm'

;
; SST now
thm_load_sst,probe=sc,lev=1
thm_part_moments, probe = sc, instrument = ['ps?f'], $
                  moments = ['density', 'velocity', 't3'], $
                  mag_suffix='_peir_magt3',          $
                  scpot_suffix='_peir_sc_pot';,/median
; work in gsm
thm_cotrans,'th'+sc+'_ps?f_velocity',in_coord='dsl',out_coord='gsm',out_suffix='_gsm'
;
; ESA now
thm_load_esa,probe=sc
; or create from scratch
;thm_load_esa_pkt,probe=sc
;thm_part_moments, probe = sc, instrument = ['pe?r'], $
;                  moments = ['density', 'velocity', 't3'], $
;                 mag_suffix='_peir_magt3',          $
;                  scpot_suffix='_peir_sc_pot'
; work in gsm
;thm_cotrans,'th'+sc+'_pe?r_velocity',in_coord='dsl',out_coord='gsm',out_suffix='_gsm'
;


fix_dat,'*'

  ;... degap all data
  ;this puts nans in the right spots so that it is easier to interpolate
  ;otherwise interpolation might generate data inside gaps
  ;density
   tdegap,'th'+sc+'_pe?r_density',/overwrite,margin=1.5
   tdegap,'th'+sc+'_ps?f_density',/overwrite,margin=150
  ;temperature
   tdegap,'th'+sc+'_pe?r_t3',/overwrite,margin=1.5
   tdegap,'th'+sc+'_ps?f_t3',/overwrite,margin=150
  ;velocity
   tdegap,'th'+sc+'_pe?r_velocity_gsm',/overwrite,margin=1.5
   tdegap,'th'+sc+'_ps?f_velocity_gsm',/overwrite,margin=150
  ;magnetic field
   tdegap,'th'+sc+'_fgs*',/overwrite


; ...now interpolate
; ...matching to esa cadence


;magnetic field
 tinterpol_mxn,'th'+sc+'_fgs_gsm','th'+sc+'_peir_density',/overwrite,/nan_extrapolate
;density
 tinterpol_mxn,'th'+sc+'_peer_density','th'+sc+'_peir_density',/overwrite,/nan_extrapolate
 tinterpol_mxn,'th'+sc+'_ps?f_density','th'+sc+'_peir_density',/overwrite,/nan_extrapolate
;temperature
 tinterpol_mxn,'th'+sc+'_pe?r_t3','th'+sc+'_peir_density',/overwrite,/nan_extrapolate
 tinterpol_mxn,'th'+sc+'_ps?f_t3','th'+sc+'_peir_density',/overwrite,/nan_extrapolate
;velocity
 tinterpol_mxn,'th'+sc+'_pe?r_velocity_gsm','th'+sc+'_peir_density',/overwrite,/nan_extrapolate
 tinterpol_mxn,'th'+sc+'_ps?f_velocity_gsm','th'+sc+'_peir_density',/overwrite,/nan_extrapolate
;
; --------------------------------------------------------------------------------
; ... Joint moments calculation
;

get_data, 'th'+sc+'_psif_density', data=sst_i_n ; sst density
get_data, 'th'+sc+'_psif_velocity_gsm', data=sst_i_v ; ion velocity (GSM)
get_data, 'th'+sc+'_psif_t3', data=sst_i_t3 ; temperature

get_data, 'th'+sc+'_peir_density', data=esa_i_n ; esa ion density
get_data, 'th'+sc+'_peir_velocity_gsm', data=esa_i_V ; ion velocity (GSM)
get_data, 'th'+sc+'_peir_t3', data=esa_i_T

get_data, 'th'+sc+'_peer_density', data=esa_e_n ; esa electron density
get_data, 'th'+sc+'_peer_t3', data=esa_e_T ; esa electron temperature

;sst_i_range = where((sst_i_n.x GT tbeg) AND (sst_i_n.x LE tend))
;esa_i_range = where((esa_i_n.x GT tbeg) AND (esa_i_n.x LE tend))
;esa_e_range = where((esa_e_n.x GT tbeg) AND (esa_e_n.x LE tend))

; ...sst Flux
sstFi = sst_i_v.y*0.
sstFi[*,0] = sst_i_n.y*sst_i_v.y[*,0]
sstFi[*,1] = sst_i_n.y*sst_i_v.y[*,1]
sstFi[*,2] = sst_i_n.y*sst_i_v.y[*,2]

; ...esa Flux
esaFi = esa_i_v.y*0.
esaFi[*,0] = esa_i_n.y*esa_i_v.y[*,0]
esaFi[*,1] = esa_i_n.y*esa_i_v.y[*,1]
esaFi[*,2] = esa_i_n.y*esa_i_v.y[*,2]

; ...total ion density
totNi = sst_i_n.y + esa_i_n.y

store_data, 'th'+sc+'_Ni', data={x:esa_i_n.x, y:totNi}
options, 'th'+sc+'_Ni', 'ytitle', 'Ni !C!C1/cm!U3'
ylim, 'th'+sc+'_Ni', 0.01, 1., 1

; ...total ion velocity (GSM)
totVi = esa_i_v.y*0.
totVi[*,0] = (sstFi[*,0]+esaFi[*,0])/totNi
totVi[*,1] = (sstFi[*,1]+esaFi[*,1])/totNi
totVi[*,2] = (sstFi[*,2]+esaFi[*,2])/totNi

store_data, 'th'+sc+'_Vi_gsm', data={x:esa_i_n.x, y:totVi}
options, 'th'+sc+'_Vi_gsm', 'ytitle', 'Vi !C!Ckm/s'

; ...pressure
; ...SST: perpendicular temperature is utilized only
sst_Tperp = .5*(sst_i_t3.y[*,0]+sst_i_t3.y[*,1])
sst_i_p_nPa = 0.16*.001*sst_i_n.y * sst_Tperp ; perp. pressure in nPa
store_data, 'th'+sc+'_psif_p_perp_nPa', data={x:sst_i_n.x, y:sst_i_p_nPa}
options, 'th'+sc+'_psif_p_perp_nPa', 'ytitle', 'sst Pi !C!CnPa'

; ...ESA: scalar temperature
esa_Ti = total(esa_i_T.y,2)/3.
store_data,'Ti_th'+sc+'_peir', data={x:esa_i_n.x, y:esa_Ti}
; ...ESA ion pressure:
esa_i_p_nPa = 0.16 *.001 * esa_i_n.y*esa_Ti ; scalar pressure in nPa
store_data, 'th'+sc+'_peir_p_nPa', data={x:esa_i_n.x, y:esa_i_p_nPa}
options, 'th'+sc+'_peir_p_nPa', 'ytitle', 'esa Pi !C!CnPa'

; ...Total ion pressure
totPi = sst_i_p_nPa + esa_i_p_nPa
store_data, 'th'+sc+'_i_p_nPa', data={x:esa_i_n.x, y:totPi}
options, 'th'+sc+'_i_p_nPa', 'ytitle', 'Pi !C!CnPa'

; ...ESA electron pressure
esa_Te = esa_e_T.y
store_data,'Te_th'+sc+'_peer', data={x:esa_e_n.x, y:esa_Te}
; ...ESA electron pressure
esa_e_p_nPa = 0.16 *.001 * esa_e_n.y*esa_Te ; scalar pressure in nPa
store_data, 'th'+sc+'_peer_p_nPa', data={x:esa_e_n.x, y:esa_e_p_nPa}
options, 'th'+sc+'_peer_p_nPa', 'ytitle', 'esa Pe !C!CnPa'

; ...Total particle pressure (SST electron pressure is ignored)
totPp = totPi + esa_e_p_nPa
store_data, 'th'+sc+'_Pp_nPa', data={x:esa_e_n.x, y:totPp}
options, 'th'+sc+'_Pp_nPa', 'ytitle', 'Pp !C!CnPa'


; ...Magnetic pressure
get_data, 'th'+sc+'_fgs_gsm', data=B
Btot = sqrt(B.y[*,0]^2+B.y[*,1]^2+B.y[*,2]^2)
Pm= 3.977e-4 * Btot^2 ; nPa
store_data, 'th'+sc+'_pmag_nPa', data={x:B.x, y:Pm}
options, 'th'+sc+'_pmag_nPa', 'ytitle', 'Pm !C!CnPa'
Ptot_nPa = Pm + totPp
store_data, 'th'+sc+'_ptot_nPa', data={x:esa_e_n.x, y:Ptot_nPa}
options, 'th'+sc+'_ptot_nPa', 'ytitle', 'Pt !C!CnPa'


; ... SC position info

get_data,strjoin('th'+sc+'_state_pos'),data=tmp
store_data,strjoin('th'+sc+'_state_pos_gsm_x'),data={x:tmp.x,y:tmp.y(*,0)/6370.}
	options,strjoin('th'+sc+'_state_pos_gsm_x'),'ytitle','th'+sc+'_X-GSM'
store_data,strjoin('th'+sc+'_state_pos_gsm_y'),data={x:tmp.x,y:tmp.y(*,1)/6370.}
	options,strjoin('th'+sc+'_state_pos_gsm_y'),'ytitle','th'+sc+'_Y-GSM'
store_data,strjoin('th'+sc+'_state_pos_gsm_z'),data={x:tmp.x,y:tmp.y(*,2)/6370.}
	options,strjoin('th'+sc+'_state_pos_gsm_z'),'ytitle','th'+sc+'_Z-GSM'


; ...plot the data. To avoid plotting unneccessary variables, pls comment corresponding lines)
;
; Densities first
;
store_data,'th'+sc+'_Ni_all',data='th'+sc+'_peir_density th'+sc+'_psif_density th'+sc+'_Ni'
tvectot,'th'+sc+'_fgs_gsm',newname='th'+sc+'_fgs_gsmt'
tplot, ['th'+sc+'_fgs_gsmt', $
        'th'+sc+'_psif_en_eflux', $
        'th'+sc+'_pe?r_en_eflux', $ 
        'th'+sc+'_Ni_all'], $
Title='THEMIS '+ sc, $
var_label=['th'+sc+'_state_pos_gsm_z', $
           'th'+sc+'_state_pos_gsm_y', $
           'th'+sc+'_state_pos_gsm_x']
tlimit,Trange
;
; Velocities next
;
options,'th'+sc+'_Vi_gsm',colors=['b','g','r']
tplot, ['th'+sc+'_fgs_gsmt', $
        'th'+sc+'_psif_en_eflux', $
        'th'+sc+'_peir_en_eflux', $ 
        'th'+sc+'_peir_velocity_gsm', $
        'th'+sc+'_psif_velocity_gsm', $
        'th'+sc+'_Vi_gsm'], $
Title='THEMIS '+ sc, $
var_label=['th'+sc+'_state_pos_gsm_z', $
           'th'+sc+'_state_pos_gsm_y', $
           'th'+sc+'_state_pos_gsm_x']
;
; Pressures next
;
tplot, ['th'+sc+'_Ni', $
        'th'+sc+'_psif_en_eflux', $
        'th'+sc+'_peir_en_eflux', $ 
        'th'+sc+'_psif_p_perp_nPa', $
        'th'+sc+'_peir_p_nPa', $
        'th'+sc+'_i_p_nPa', $
        'th'+sc+'_peer_en_eflux', $ 
        'th'+sc+'_peer_p_nPa', $
        'th'+sc+'_Pp_nPa',$
        'th'+sc+'_pmag_nPa', $
        'th'+sc+'_ptot_nPa'], $
Title='THEMIS '+ sc, $
var_label=['th'+sc+'_state_pos_gsm_z', $
           'th'+sc+'_state_pos_gsm_y', $
           'th'+sc+'_state_pos_gsm_x']

; ...asii output, uncomment if necessary
; ... Extract ASCII:
;tplot_ascii, 'th'+sc+'_state_pos', dir='E:\work\themis\ascii\'+date
;tplot_ascii, 'th'+sc+'_fgs_gsm', Dir='E:\work\themis\ascii\'+date+'\'
;tplot_ascii, 'th'+sc+'_Ni', Dir='E:\work\themis\ascii\'+date+'\'
;tplot_ascii, 'th'+sc+'_Vi_gsm', Dir='E:\work\themis\ascii\'+date+'\'
;tplot_ascii, 'th'+sc+'_Pp_nPa', Dir='E:\work\themis\ascii\'+date+'\'
;tplot_ascii, 'th'+sc+'_ptot_nPa', Dir='E:\work\themis\ascii\'+date+'\'
;     the text files with the data will appear in the subdirectory
;      specified by Dir variable (specify your PATH please!)

;

end