;+
;PROCEDURE: thm_load_peer
;PURPOSE:
; Decommutates the raw THEMIS electron esa reduced data packets, (peer - apid 458), and puts data structure in a common for access by get_th*_peer.pro and generates tplot energy spectrogram
;INPUT:
;
;KEYWORDS:
; file: string, strarr the complete path and filename to the raw pkt files
; if not set, uses timerange to select files
; sc: string themis spacecraft - "a", "b", "c", "d", "e"
; if not set defaults to "a"
; themishome: string path to data dir, where data dir contains the th* dir, where *=a,b,c,d,e
; gap_time real data gaps greater than gap_time have a NAN replacing a data point in tplot structures
; suffix string "suffix" is appended to the default tplot names
;
;CREATED BY: J. McFadden 07/03/18
;VERSION: 1
;LAST MODIFICATION: 08/05/14
;MOD HISTORY:
; suffix keyword added 08/05/14
;
;NOTES:
;
;-
pro thm_load_peer,file=file,sc=sc,themishome=themishome, $
gap_time=gap_time,suffix=suffix,trange=trange, $
use_eclipse_corrections=use_eclipse_corrections
; suffix default
if ~keyword_set(suffix) then suffix = ''
; sc default
if not keyword_set(sc) then begin
dprint, 'S/C number not set, default = a'
sc='a'
endif
; zero data common
if sc eq 'a' then begin
common tha_458,tha_458_ind,tha_458_dat & tha_458_dat=0 & tha_458_ind=-1l
endif else if sc eq 'b' then begin
common thb_458,thb_458_ind,thb_458_dat & thb_458_dat=0 & thb_458_ind=-1l
endif else if sc eq 'c' then begin
common thc_458,thc_458_ind,thc_458_dat & thc_458_dat=0 & thc_458_ind=-1l
endif else if sc eq 'd' then begin
common thd_458,thd_458_ind,thd_458_dat & thd_458_dat=0 & thd_458_ind=-1l
endif else if sc eq 'e' then begin
common the_458,the_458_ind,the_458_dat & the_458_dat=0 & the_458_ind=-1l
endif else if sc eq 'f' then begin
common thf_458,thf_458_ind,thf_458_dat & thf_458_dat=0 & thf_458_ind=-1l
endif
if not keyword_set(themishome) then themishome=!themis.local_data_dir
; get filenames if file keyword not set
if not keyword_set(file) then begin
tt=timerange(trange)
t1=time_double(strmid(time_string(tt(0)),0,10))
t2=time_double(strmid(time_string(tt(1)-1.),0,10))
ndays=1+fix((t2-t1)/(24.*3600.))
; print,ndays
file=strarr(ndays)
i=0
while t1 le t2 do begin
ts=time_string(t1)
yr=strmid(ts,0,4)
mo=strmid(ts,5,2)
da=strmid(ts,8,2)
if i eq 0 then yrmoda1=yr+mo+da
yrmoda2=yr+mo+da
dir=themishome+'th'+sc+'/l0/'+yr+'/'+mo+'/'+da+'/'
name='th'+sc+'_l0_458_'+yr+mo+da+'.pkt'
file[i]=dir+name
i=i+1
t1=t1+24.*3600.
endwhile
endif
; check that files exist
nfiles=n_elements(file)
if nfiles eq 1 then begin
if not file_test(file) then begin
dprint, file+' --- does not exist.'
return
endif
endif else begin
ind=-1
for i=0,nfiles-1 do begin
if file_test(file[i]) then ind=[ind,i] else dprint, file[i]+' --- does not exist.'
endfor
n_ind=n_elements(ind)
if n_ind eq 1 then return else begin
ind=ind[1:n_ind-1]
file=file[ind]
endelse
endelse
; decompression array
decompress=1.*[indgen(16),16+2*indgen(8),32+4*indgen(24),128+8*indgen(16),256+16*indgen(16),$
512+32*indgen(16),1024+64*indgen(16),2048+128*indgen(48),8192+256*indgen(32),$
16384+512*indgen(32),32768+1024*indgen(32)]
; define mode decoder
; n_modes=5 ; number of modes (not used as of 2014-03-26 JWL)
mode_decode=intarr(256,256) ; decode array
mode_decode[*,*]=-1
mode_decode[0,0]=0 ; 1A x 16E, snapshot, 1 spin
mode_decode[1,1]=1 ; 1A x 32E, snapshot, 1 spin
mode_decode[2,1]=2 ; 6A x 32E, snapshot, 1 spin
mode_decode[1,2]=3 ; 1A x 32E, snapshot, 1 spin
mode_decode[2,2]=4 ; 6A x 32E, snapshot, 1 spin
mode_decode[3,1]=5 ; 1A x 32E, snapshot, 1 spin
;magnetospheric slow/fast modes with reduced energy map - added ????-??-??
mode_decode[1,4]=6 ; 1A x 32E, snapshot, 1 spin
mode_decode[2,4]=7 ; 6A x 32E, snapshot, 1 spin
;Another low-energy magnetospheric mode 0x0304 JWL 2013-07-24
mode_decode[3,4]=8 ; 1A x 32E, snapshot, 1 spin
;Revised low-energy magnetospheric modes 0x0105, 0x0205
mode_decode[1,5]=9 ; 1A x 32E, snapshot, 1 spin
mode_decode[2,5]=10 ; 6A x 32E, snapshot, 1 spin
;magnetospheric slow/fast modes, low E - added 2016-03-18
mode_decode[1,6]=11 ; 1A x 32E, snapshot, 1 spin
mode_decode[2,6]=12 ; 6A x 32E, snapshot, 1 spin
; define mode variables for different modes
nspins = [1,1,1,1,1,1,1,1,1,1,1,1,1] ; # of spins between measurements in mode
nenergy = [16,32,32,32,32,32,32,32,32,32,32,32,32] ; # of energies in mode
nangle = [1,1,6,1,6,1,1,6,1,1,1,1,6] ; # of angles in mode
dat_len = nenergy*nangle ; size data arrays
spin_decode = [160,96,16,96,16,96,96,16,96,96,96,96,16]; # measurements in packet
case_decode = [0,1,3,1,3,1,1,3,1,1,1,1,3] ; datl[16,32,96,192]==>size[0,1,2,3]
angle_decode = [0,1,2,1,2,1,1,2,1,1,1,1,2] ; angle mode index
energy_decode = [0,1,1,2,2,1,3,3,3,4,4,5,5] ; energy mode index
; initialize arrays
ndays=n_elements(file)
nmax=40000l*ndays
dat_0=bytarr(nmax,16)
dat_1=bytarr(nmax,32)
dat_2=bytarr(nmax,96)
dat_3=bytarr(nmax,192)
n_case=lonarr(4) & n_case[*]=0 ; # elements in n_case is number of dat_? arrays
config1=bytarr(nmax)
config2=bytarr(nmax)
an_ind = intarr(nmax)
en_ind = intarr(nmax)
md_ind = intarr(nmax)
cs_ind = intarr(nmax)
cs_ptr = lonarr(nmax)
s_time = dblarr(nmax)
e_time = dblarr(nmax)
data_size=0l
; get the files
tm1=systime(1)
dprint, dlevel=2, 'Downloading data'
if n_elements(file) eq 1 then openr,fp,file,/get_lun else openr,fp,file[0],/get_lun
fs = fstat(fp)
; printdat,fs
if fs.size ne 0 then begin
adat = bytarr(fs.size)
data_size=fs.size*1l
readu,fp,adat
free_lun,fp
endif
if n_elements(file) ne 1 then begin
nfile=n_elements(file)
for i=1,nfile-1 do begin
openr,fp,file[i],/get_lun
fs = fstat(fp)
; printdat,fs
if fs.size ne 0 then begin
adat2 = bytarr(fs.size)
readu,fp,adat2
if data_size ne 0 then adat=[adat,adat2] else adat=adat2
data_size=data_size+fs.size
free_lun,fp
endif
endfor
endif
tm2=systime(1)
dprint, dlevel=2, 'Data download complete'
dprint, dlevel=4, '458 Download time= ',tm2-tm1
if data_size eq 0 then return
; Calculate spin period using header times the 3rd and 4th packets,
; assume spin period does not change within data.
;; This spin period will be used if spin model data is not present.
;-----------------------------------------
hdrl=16 ; header length
pntr=0l
pktl=0l
npkt=0l
ntot=0l
tf_spin=1 ;;track if spin period has been successfully determined
last_spin_period=0.
i=0
while (tf_spin and i lt 100 and pntr+hdrl+pktl lt data_size and data_size gt 4000) do begin
i=i+1
; print,tf_spin,i
hdr = adat[pntr:pntr+hdrl-1] ;;header
time = (1.d*((hdr[6]*256ul+hdr[7])*256ul+hdr[8])*256ul+hdr[9])+(hdr[10]+hdr[11]/256.d)/256.
pktl = hdr[4]*256+hdr[5]+7-hdrl ;;packet length
config = [adat[pntr+hdrl],adat[pntr+hdrl+1]] ;;get index into mode table
mode=mode_decode[config[0],config[1]] ;;get mode from table
pntr = pntr+pktl+hdrl ;;increment to next packet
hdr2 = adat[pntr:pntr+hdrl-1]
time2 = (1.d*((hdr2[6]*256ul+hdr2[7])*256ul+hdr2[8])*256ul+hdr2[9])+(hdr2[10]+hdr2[11]/256.d)/256.
pktl2 = hdr2[4]*256+hdr2[5]+7-hdrl
configg = [adat[pntr+hdrl],adat[pntr+hdrl+1]]
mode2 = mode_decode[configg[0],configg[1]]
;;calculate spin period
;;if successful the loop will be exited
if pktl eq pktl2 then begin
if mode eq mode2 and mode ne -1 then begin
;; dt/(n_spins * n_measurements/packet), n_spins always 1
spin_period = (time2-time)/(nspins[mode]*spin_decode[mode])
if spin_period lt 6. and spin_period gt 1.5 $
and abs(spin_period-last_spin_period) lt .005 $
then tf_spin=0 $ ;;spin per in (1.5,6) & close to previous
else last_spin_period=spin_period ;;spin per outside range, different from last, or first run
endif
endif
; print,pntr,' ',pktl,' ',time_to_str(time),' ',adat[pntr+hdrl],' ',adat[pntr+hdrl+1],' spin period= ',spin_period
endwhile
; if tf_spin then print,'ERROR : spin period calculation error, setting spin period to 3. seconds'
;; If spin period was never in (1.5,6) or varied >.005 from previous each time
;; or while loop only ran once then the spin period is set to 3.
if tf_spin then spin_period=3.
dprint, dlevel=4, 'spin_period = ',spin_period
last_spin_period = spin_period
tm1=systime(1)
;; Load spin data
; spinmodel=thm_load_spinmodel(sc=sc,themishome=themishome,available=spindata)
; 2012-10-16
; Load spin model
; Support data should already be loaded at this point
;-------------------------------------
model = spinmodel_get_ptr(sc,use_eclipse_corrections=use_eclipse_corrections)
if ~obj_valid(model) then begin
dprint, dlevel=2, 'Using default spin period: '+strtrim(spin_period,2)+' sec'
endif
; Start decommutation
;-----------------------------------------
; help,adat
dprint, dlevel=4, 'Decommutating data ... Data volume (bytes) =',data_size
; print,'Current data pktl time config bytes '
; initialize variables
t0=time_double('2001-1-1/0')
pntr=0l ;;pointer to current location in file
pktl=0l ;;packet length
npkt=0l ;;number of packets?
ntot=0l ;;index into initialized data
hdr = adat[pntr:pntr+hdrl-1] ;;header
last_apid_cntr= hdr[3]
last_time = (1.d*((hdr[6]*256ul+hdr[7])*256ul+hdr[8])*256ul+hdr[9]) $ ;; ??
+ (hdr[10]+hdr[11]/256.d)/256.
;;loop over packet file
while pntr+1 lt data_size do begin
hdr = adat[pntr:pntr+hdrl-1] ;;get header for this packet
apid_cntr= hdr[3]
spin_num = hdr[15] ;;n spins
;;get header time and copy to new var
time = (1.d*((hdr[6]*256ul+hdr[7])*256ul+hdr[8])*256ul+hdr[9])+(hdr[10]+hdr[11]/256.d)/256.
time2 = time
;;if spin model data is present:
;; -shift the copied time to be at the begining of the most recent spin
;; -set the spin period to that of the model
; if spindata then begin
;
; if n_elements(spinmodel.s_time) eq 1 then begin
;
; ;; rounded # full spins between time and spin_time
; spin_offset = long((time-spinmodel.s_time)/spinmodel.spin_period+.5d)
; ;; spin_time + time taken for rounded # of spins
; time2 = spinmodel.s_time + spin_offset * spinmodel.spin_period
; ;; subtract one spin from new time if original time exceeded
; if time2 gt time then time2 = time2 - spinmodel.spin_period
;
; spin_period = spinmodel.spin_period ;; set spin per
;
;; offset=2.00+(6./32.)*(spin_period-3.)
;; if apid_cntr ne last_apid_cntr then print,'ERROR: peer apid counter discontinuity : '$
;; ,time2-time,-1.*offset,apid_cntr,last_apid_cntr,' th'+sc,' ',time_string(last_time+t0),' '$
;; ,time_string(time+t0),' ',spinmodel.spin_period,' ',pntr
;; if abs(time2-time+offset) gt .15 then print,'ERROR: peer Header time differs from spin model spin-start time by: '$
;; ,time2-time,-1.*offset,apid_cntr,last_apid_cntr,' th'+sc,' ',time_string(last_time+t0),' '$
;; ,time_string(time+t0),' ',spinmodel.spin_period,' ',pntr
;; print,spinmodel.s_time,time2,time,time2-time
;
; endif else begin
;
; ;; min dt between time and spin time
; mintm = min(abs(spinmodel.s_time-time),ind)
; ;; use spin_time closest to but not exceeding time
; if spinmodel.s_time[ind] gt time and ind ne 0 then ind = ind-1
;
; ;; rounded # full spins between time and spin_time
; spin_offset = long((time-spinmodel.s_time[ind])/spinmodel.spin_period[ind]+.5d)
; ;; spin_time + time taken for roudned # of spins
; time2 = spinmodel.s_time[ind] + spin_offset * spinmodel.spin_period[ind]
; ;; subract one spin from new time if original time exceeded
; if time2 gt time then time2 = time2-spinmodel.spin_period[ind]
;
; spin_period = spinmodel.spin_period[ind] ;; set spin per
;
;; offset=2.00+(6./32.)*(spin_period-3.)
;; if apid_cntr ne last_apid_cntr then print,'ERROR: peer apid counter discontinuity : '$
;; ,time2-time,-1.*offset,apid_cntr,last_apid_cntr,' th'+sc,' ',time_string(last_time+t0),' '$
;; ,time_string(time+t0),' ',spinmodel.spin_period(ind),' ',pntr
;; if abs(time2-time+offset) gt .15 then print,'ERROR: peer Header time differs from spin model spin-start time by: '$
;; ,time2-time,-1.*offset,apid_cntr,last_apid_cntr,' th'+sc,' ',time_string(last_time+t0),' '$
;; ,time_string(time+t0),' ',spinmodel.spin_period(ind),' ',pntr
;; print,spinmodel.s_time(ind),time2,time,time2-time
;
; endelse
;
; endif
;2012-10-16
;If spin model is present:
; -get spin period
; -set adjusted time to be at the begining of the most recent spin
if obj_valid(model) then begin
thm_esa_spin_adjust, model=model, time=time, adjusted_time=time2, $
spin_period=spin_period
endif
pktl = hdr[4]*256+hdr[5]+7-hdrl
pktp = 0l
datl = 0
pkt = adat[pntr+hdrl:pntr+hdrl+pktl-1]
sub_spin=0
last_mode=mode_decode[pkt[0],pkt[1]]
last_apid_cntr= byte(apid_cntr+1)
mode_change=0
;;loop over packet looking for valid samples
;; -set start/end time for each measurement
;; -check for mode change and add timing correction if found
;; (this can only occur once through entire loop; all
;; remaining data in the packet will be corrected)
while pktp+4 lt pktl do begin
config_1 = pkt[pktp] ;; get first two bytes (config word)
config_2 = pkt[pktp+1] ;;
mode = mode_decode[config_1,config_2] ;;determine mode from config word
if mode ne -1 then datl = dat_len[mode] ;;get length of data portion
if mode ne -1 and pktp+datl+4 le pktl then begin ;skip if invalid or incomplete
; datl = dat_len[mode]
navg = pkt[pktp+3] ;; # spins per measure from packet data
nsp = nspins[mode] ;; # spins per measure based on mode
config1[ntot] = config_1 ;; config word
config2[ntot] = config_2 ;;
an_ind[ntot] = angle_decode[mode] ;;angle mode index
en_ind[ntot] = energy_decode[mode] ;;energy mode index
md_ind[ntot] = mode ;;mode index
nc = case_decode[mode]
cs_ind[ntot] = nc
cs_ptr[ntot] = n_case[nc]
;;check for mid-packet mode change, specifically check if:
;; -current mode different from previous
;; -next header in current file
;; -first two conditions not previously met in this loop
;; (i.e. code in this block gets hit once, at most, while looping over current packet)
if mode ne last_mode and pntr+pktl+hdrl+hdrl-1 lt data_size and mode_change ne 1 then begin
;;get next packet's header time
hdr = adat[pntr+pktl+hdrl:pntr+pktl+hdrl+hdrl-1]
; print,fix((time-last_time)/spin_period+.5),spin_num,hdr[15]
mode_change=1 ;;set flag
next_time = (1.d*((hdr[6]*256ul+hdr[7])*256ul+hdr[8])*256ul+hdr[9])+(hdr[10]+hdr[11]/256.d)/256.
;;get # of remaining samples in current packet
n_rem = (pktl - pktp)/(datl+4)
;;increment if remainder was truncated
if n_rem*(datl+4) ne (pktl-pktp) then begin
n_rem=n_rem+1
endif
;;get # spins between initial header time and mode change:
;; rounded total spins in packet - whole # of spins to reach next packet from mode change
rem_spin = fix((next_time-time)/spin_period+.5) - n_rem*nsp
;;unneeded, set later
time0 = time_double('1900-1-1')
if next_time lt 3.e9 then time0 = time_double('2001-1-1') ; if 2001 epoch used
; print,'Mode change in packet: ',time_string(next_time+time0),mode,last_spin,rem_spin, n_rem,pktl,pktp $
; ,datl,nsp,(1.*pktl-pktp)/(datl+4),fix((next_time-time)/spin_period+.5)
;;select larger of:
;; # spins encountered up to this point
;; # total spins in packet - spins remaining after mode change (based on calculation above)
sub_spin = rem_spin > sub_spin
; print,'sub_spin = ',sub_spin
endif
;;set sample start/end times for common block structure
;; spin-model-corrected header time + spin period * rounded number of spins to current location
s_time[ntot] = time2 + spin_period * (sub_spin-1)
e_time[ntot] = time2 + spin_period * (sub_spin+navg)
tmp=reform(pkt[pktp+4:pktp+4+datl-1],1,datl)
case nc of
0: dat_0[n_case[nc],*]=tmp
1: dat_1[n_case[nc],*]=tmp
2: dat_2[n_case[nc],*]=tmp
3: dat_3[n_case[nc],*]=tmp
endcase
ntot = ntot + 1 ;;increment data index
n_case[nc] = n_case[nc] + 1
endif else begin
pktp = pktl
; print,'Error: Invalid ESA config header bytes - skipping packet'
; print,mode,config_1,config_2,pntr,pktl,pktp
datl=0
nsp=1
endelse
pktp = pktp + 4 + datl ;;go to next sample
last_spin = sub_spin
sub_spin = sub_spin + nsp
last_mode = mode
; if ntot-1000l*ceil(ntot/1000l) eq 0 then print,pntr,' ',pktl,' ',time_string(time),' ',adat[pntr+hdrl],' ',adat[pntr+hdrl+1]
endwhile ;;end loop over current packet
pntr = pntr + pktl + hdrl ;;increment main pointer
; print,pntr,' ',ntot
last_time=time
endwhile ;;end loop over data file
tm2=systime(1)
dprint, dlevel=4, '458 Decommutation time= ',tm2-tm1
if ntot eq 0 then return
; Get rid of extra elements in arrays
;-----------------------------------------
config1 = config1[0:ntot-1]
config2 = config2[0:ntot-1]
an_ind = an_ind[0:ntot-1]
en_ind = en_ind[0:ntot-1]
md_ind = md_ind[0:ntot-1]
cs_ind = cs_ind[0:ntot-1]
cs_ptr = cs_ptr[0:ntot-1]
s_time = s_time[0:ntot-1]
e_time = e_time[0:ntot-1]
if n_case[0] gt 0 then dat_0=dat_0[0:n_case[0]-1,*] else dat_0=dat_0[0,*]
if n_case[1] gt 0 then dat_1=dat_1[0:n_case[1]-1,*] else dat_1=dat_1[0,*]
if n_case[2] gt 0 then dat_2=dat_2[0:n_case[2]-1,*] else dat_2=dat_2[0,*]
if n_case[3] gt 0 then dat_3=dat_3[0:n_case[3]-1,*] else dat_3=dat_3[0,*]
; Get eclipse delta phi from spin model for each sample time
; If the model is not available the field will be filled with NaNs
if obj_valid(model) then begin
thm_esa_spin_adjust, model=model, time=s_time, eclipse_dphi=eclipse_dphi
endif else begin
eclipse_dphi = replicate(!values.d_nan,size(s_time,/dim))
endelse
; shift the time to unix time
time0 = time_double('1900-1-1')
if s_time[ntot-1 < 5] lt 3.e9 then time0 = time_double('2001-1-1') ; if 2001 epoch used
s_time = s_time + time0
e_time = e_time + time0
;***********************************************************************************
; may want to fix this section
valid=bytarr(ntot) & valid[*]=1
; ind=where(s_time gt e_time),count)
; if count gt 0 then valid[ind]=0
;***********************************************************************************
;Energy mode determination
;-----------------------------------------
case sc of
'a': fm=1
'b': fm=2
'c': fm=5
'd': fm=4
'e': fm=3
'f': fm=6
endcase
emode=thm_read_esa_sweep_reduced_mode(fm)
energy=emode.e_energy
denergy=emode.e_denergy
nenergy=emode.e_nenergy
;***********************************************************************************
; Calibration: geometric factor and deadtime and efficiency
calib=get_thm_esa_cal(sc=sc,ion=0,time=s_time)
geom_factor=calib.geom_factor ; geom_factor of 22.5x22.5 angle bin
an_eff = calib.an_eff
rel_gf = calib.rel_gf
an_en_eff = calib.an_en_eff
deadtime=1.7e-7
; Make electron energy efficiency array
; REF: Relative electron detection efficiency of microchannel plates from 0-3 keV
; R.R. Goruganthu and W. G. Wilson, Rev. Sci. Instrum. Vol. 55, No. 12 Dec 1984
tmax = 2.283
aa = 1.35
delta_max = 1.0
emax = 325. ;eV
en_acc = 450. ;eV for THEMIS
delta = delta_max * ((energy+en_acc)/emax)^(1-aa)*(1-exp(-tmax*((energy+en_acc)/emax)^aa))/(1-exp(-tmax))
k = 2.2
en_eff = (1-exp(-k*delta/delta_max))/(1-exp(-k))
;***********************************************************************************
amode=thm_read_eesa_angle_reduced_mode()
nbins=amode.e_nbins
theta=amode.e_theta
dtheta=amode.e_dtheta
phi=amode.e_phi
dphi=amode.e_dphi
domega=amode.e_domega
nrg_wt=amode.e_nrg_wt
anodes=amode.e_anodes
an_map=amode.e_an_map
;***********************************************************************************
; to get magnetometer data, use thm_load_esa_mag.pro
magf=fltarr(ntot,3) & magf[*]=!values.f_nan
;***********************************************************************************
; gf, dt_arr, dt
dt_arr=(dphi/11.25)*anodes*nrg_wt
gf=dphi/22.5*dtheta/22.5*nrg_wt
dt=e_time-s_time
;***********************************************************************************
; to load background data, use thm_load_pee_bkg.pro
; bkg, bkg_pse, bkg_pee are the time dependent background counts for the entire analyzer in 1/32 of a spin
; bkg is roughly the MCP rate per cm^2-s
; bkg_pse is an eSST determined background
; bkg_pee is an eESA determined background
; bkg is generally the smaller of bkg_pse, bkg_pee -- see keywords in thm_load_pee_bkg.pro
; an_bkg - relative bkg per anode - determined from sc='e' on 08-01-16,
; an_bkg should probably be a function of sc and perhaps of time
; bkg_arr is the relative background per energy-angle bin
; background arrays are constructed as: bkg*bkg_arr
bkg_pse=fltarr(ntot) & bkg_pse[*]=!values.f_nan
bkg_pee=fltarr(ntot) & bkg_pee[*]=!values.f_nan
bkg=fltarr(ntot) & bkg[*]=0.
an_bkg = [1.32653,1.12855,0.759310,0.658778,0.689404,0.764592,1.17875,1.49409]
mapdim = size(an_map)
maxbins = mapdim[2]
nummaps = mapdim[3]
tmparr = fltarr(maxbins,nummaps)
for i=0,nummaps-1 do tmparr[*,i] =total((an_bkg#replicate(1.,maxbins))*an_map[*,*,i],1)
bkg_arr = reform(replicate(1.,32)#reform(tmparr,maxbins*nummaps),32,maxbins,nummaps)
bkg_arr = bkg_arr*gf/128.
;***********************************************************************************
; Make 3D data structures
; Note: sc_pot, b_gse, magf, rot, phi_offset - are to be added later
erd_dat = {project_name: 'THEMIS', $
spacecraft: sc, $
data_name: 'EESA 3D Reduced', $
apid: '458'x, $
units_name: 'compressed', $
units_procedure: 'thm_convert_esa_units', $
valid: valid, $
time: s_time, $
end_time: e_time, $
delta_t: dt, $
integ_t: dt/1024., $
dt_arr: dt_arr, $
cs_ptr: cs_ptr, $
cs_ind: cs_ind, $
config1: config1, $
config2: config2, $
an_ind: an_ind, $
en_ind: en_ind, $
md_ind: md_ind, $
nenergy: nenergy, $
energy: energy, $
denergy: denergy, $
; eff: en_eff, $
nbins: nbins, $
theta: theta, $
dtheta: dtheta, $
phi: phi, $
dphi: dphi, $
phi_offset: replicate(0.,ntot), $
domega: domega, $
gf: gf, $
an_map: an_map, $
eclipse_dphi: eclipse_dphi, $
rel_gf: rel_gf, $
an_eff: an_eff, $
en_eff: en_eff, $
an_en_eff: an_en_eff, $
geom_factor: geom_factor, $
dead: deadtime, $
mass: 5.68566e-06, $
charge: -1., $
sc_pot: replicate(0.,ntot), $
magf: magf, $
bkg_pse: bkg_pse, $
bkg_pei: bkg_pee, $
bkg_pee: bkg_pee, $
bkg: bkg, $
bkg_arr: bkg_arr, $
dat0: dat_0, $
dat1: dat_1, $
dat2: dat_2, $
dat3: dat_3}
if sc eq 'a' then begin
common tha_458,tha_458_ind,tha_458_dat & tha_458_dat=erd_dat & tha_458_ind=0l
endif else if sc eq 'b' then begin
common thb_458,thb_458_ind,thb_458_dat & thb_458_dat=erd_dat & thb_458_ind=0l
endif else if sc eq 'c' then begin
common thc_458,thc_458_ind,thc_458_dat & thc_458_dat=erd_dat & thc_458_ind=0l
endif else if sc eq 'd' then begin
common thd_458,thd_458_ind,thd_458_dat & thd_458_dat=erd_dat & thd_458_ind=0l
endif else if sc eq 'e' then begin
common the_458,the_458_ind,the_458_dat & the_458_dat=erd_dat & the_458_ind=0l
endif else if sc eq 'f' then begin
common thf_458,thf_458_ind,thf_458_dat & thf_458_dat=erd_dat & thf_458_ind=0l
endif
;store time range & eclipse correction status
eclipse = undefined(use_eclipse_corrections) ? 0:use_eclipse_corrections
thm_part_trange, sc, 'peer', set={trange:timerange(trange),eclipse:eclipse}
; help,erd_dat,/st
dprint, dlevel=4, 'Data loading complete: th'+sc+' peer - apid 458'
;***********************************************************************************
; Make tplot structures
cols=get_colors()
en_dat=fltarr(ntot,32)
energy=fltarr(ntot,32)
for i=0l,ntot-1l do begin
; print,i,cs_ptr[i],cs_ind[i]
case cs_ind[i] of
0: en_dat[i,0:15]=reform(decompress[erd_dat.dat0[cs_ptr[i],*]])
1: en_dat[i,0:31]=reform(decompress[erd_dat.dat1[cs_ptr[i],*]])
2: en_dat[i,0:15]=total(reform(decompress[erd_dat.dat2[cs_ptr[i],*]],6,16),1)
3: en_dat[i,0:31]=total(reform(decompress[erd_dat.dat3[cs_ptr[i],*]],6,32),1)
endcase
energy[i,*]=reform(erd_dat.energy[*,en_ind[i]])
endfor
if not keyword_set(gap_time) then gap_time=10.
if ntot gt 1 then begin
ind=where(s_time[1:ntot-1]-s_time[0:ntot-2] gt gap_time,n_ind)
if n_ind gt 0 then begin
en_dat[ind,*]=!values.f_nan
en_dat[ind+1,*]=!values.f_nan
endif
endif
name1='th'+sc+'_peer_en_counts'+suffix
store_data,name1,data={x:(erd_dat.time+erd_dat.end_time)/2.,y:en_dat,v:energy}
zlim,name1,0,0,1,/def
ylim,name1,3.,40000.,1,/def
options,name1,'ztitle','Counts',/def
options,name1,'ytitle','e- th'+sc+'!C!C eV',/def
options,name1,'spec',1,/def
options,name1,'x_no_interp',1,/def
options,name1,'y_no_interp',1,/def
name1='th'+sc+'_peer_mode'+suffix
store_data,name1,data={x:(s_time+e_time)/2.,y:[[config1-.1],[config2+.1],[md_ind]]}
ylim,name1,-1,10,0,/def
options,name1,'ytitle','e- th'+sc+'!C!C mode',/def
options,name1,'colors',[cols.red,cols.green,cols.blue],/def
options,name1,labels=['config1', 'config2', 'mode'],constant=0.,/def
options,name1,'labflag',3,/def
options,name1,'labpos',[0,3,6],/def
;***********************************************************************************
end