;+ ;PROCEDURE: thm_load_peef ;PURPOSE: ; Decommutates the raw THEMIS electron esa full data packets, (peef - apid 457), and puts data structure in a common for access by get_th*_peef.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_peef,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_457,tha_457_ind,tha_457_dat & tha_457_dat=0 & tha_457_ind=-1l endif else if sc eq 'b' then begin common thb_457,thb_457_ind,thb_457_dat & thb_457_dat=0 & thb_457_ind=-1l endif else if sc eq 'c' then begin common thc_457,thc_457_ind,thc_457_dat & thc_457_dat=0 & thc_457_ind=-1l endif else if sc eq 'd' then begin common thd_457,thd_457_ind,thd_457_dat & thd_457_dat=0 & thd_457_ind=-1l endif else if sc eq 'e' then begin common the_457,the_457_ind,the_457_dat & the_457_dat=0 & the_457_ind=-1l endif else if sc eq 'f' then begin common thf_457,thf_457_ind,thf_457_dat & thf_457_dat=0 & thf_457_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_457_'+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 (unused as of 2014-03-26, JWL) mode_decode=intarr(256,256) ; decode array mode_decode[*,*]=-1 mode_decode[0,0]=0 ; 88A x 32E, snapshot, 3 spin mode_decode[1,1]=1 ; 88A x 32E, snapshot, 128 spin mode_decode[2,1]=2 ; 88A x 32E, snapshot, 32 spin mode_decode[1,2]=3 ; 88A x 32E, snapshot, 128 spin mode_decode[2,2]=4 ; 88A x 32E, snapshot, 32 spin mode_decode[3,1]=5 ; 88A x 15E, snapshot, 1 spin ;magnetospheric slow/fast modes with reduced energy map - added ????-??-?? mode_decode[1,4]=6 ; 88A x 32E, snapshot, 128 spin mode_decode[2,4]=7 ; 88A x 32E, snapshot, 32 spin ; Another magnetospheric mode JWL 2013-07-23 mode_decode[3,4]=8 ; 88A x 15E, snapshot, 32 spin ; Revised magnetospheric mode JWL 2014-03-26 mode_decode[1,5]=9 ; 88A x 32 E, snapshot, 128 spin mode_decode[2,5]=10 ; 88A x 32 E, snapshot, 1 spin ; define mode variables for different modes nspins = [3,128,32,128,32,1,128,32,1,128,1] ; # of spins between measurements in mode nenergy = [32,32,32,32,32,15,32,32,15,32,32] ; # of energies in mode nangle = [88,88,88,88,88,88,88,88,88,88,88] ; # of angles in mode dat_len = nenergy*nangle ; size data arrays spin_decode = [1,1,1,1,1,3,1,1,3,1,1] ; # measurements in packet case_decode = [0,0,0,0,0,1,0,0,1,0,0] ; datl[2816,2816,2816,2816,2816,1320]==>size[0,0,0,0,0,1] angle_decode = [0,0,0,0,0,1,0,0,1,0,0] ; angle mode index, last element should be 1, jmm, 2013-10-11 energy_decode = [0,1,1,2,2,3,4,4,5,6,6] ; energy mode index ; initialize arrays ndays=n_elements(file) nmax=40000l*ndays dat_0=bytarr(nmax, 2816) dat_1=bytarr(nmax, 1320) n_case=lonarr(2) & 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, '457 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 ;----------------------------------------------------------------- hdrl=16 ; header length pntr=0l pktl=0l npkt=0l ntot=0l tf_spin=1 last_spin_period=0. & mode=0 & every_spin=0 i=0 ;; Loop over all headers in file to determine spin period ;; This spin period will be used if spin model data is not present. 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 ; print,'Config bytes = ',config mode=mode_decode[config[0],config[1]] ;;get mode ;;increment to next packet and repeat pntr = pntr+pktl+hdrl 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 pktl eq pktl2 then begin if mode eq mode2 and mode ne -1 then begin every_spin=0 ;; dt/(n_spins * n_measurements/packet) spin_period = (time2-time)/(nspins[mode]*spin_decode[mode]) if spin_period lt 1. then begin ; assume ETC is commanded for every spin spin_period = (time2-time)/(spin_decode[mode]) every_spin=1 endif 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, diff from last, ;;or first time through loop endif endif ; print,pntr,' ',pktl,' ',time_string(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 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. time = (1.d*((hdr[6]*256ul+hdr[7])*256ul+hdr[8])*256ul+hdr[9])+(hdr[10]+hdr[11]/256.d)/256. ;;this variable does not appear to be used, unsure of its purpose if mode ne -1 then begin if every_spin then begin last_time = time-spin_period*spin_decode[mode] endif else begin last_time = time-spin_period*nspins[mode]*spin_decode[mode] endelse endif ;;loop over packet file while pntr+1 lt data_size do begin hdr = adat[pntr:pntr+hdrl-1] ;;get header for current packet apid_cntr= hdr[3] spin_num = hdr[15] ;;n spins ;;get time and copy 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 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 ; ;; set spin period ; spin_period=spinmodel.spin_period ;; offset=2.12+(6./32.)*(spin_period-3.) ;; if apid_cntr ne last_apid_cntr then print,'ERROR: peef 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: peef 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 ; ;; use spin_time closest to but not exceeding time ; mintm=min(abs(spinmodel.s_time-time),ind) ; 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] ; ;; set spin period ; spin_period=spinmodel.spin_period[ind] ;; offset=2.12+(6./32.)*(spin_period-3.) ;; if apid_cntr ne last_apid_cntr then print,'ERROR: peef 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: peef 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=last_mode if mode ne -1 then datl=dat_len[mode] last_time=time ;;loop over all samples in this packet to populate data arrays while pktp+datl+4 le pktl do begin config_1 = pkt[pktp] config_2 = pkt[pktp+1] mode = mode_decode[config_1,config_2] ;;determine mode from config word ;;skip if invalid or mode has changed if mode ne -1 and mode eq last_mode then begin datl = dat_len[mode] navg = pkt[pktp+3] ;; # spins per measure 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] ;;set sample start/end times for common block structure ;; spin-model-corrected header time + spin period * number of spins to current sample 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 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 or mode change in packet - skipping packet' ; print,mode,config_1,config_2,pntr,pktl,pktp datl=0 & nsp=1 endelse pktp = pktp + 4 + datl ;;go to next sample 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 pntr = pntr + pktl + hdrl ;;go to next header/packet ; print,pntr,' ',ntot last_time=time endwhile tm2=systime(1) dprint, dlevel=4, '457 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,*] ; 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_full_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)) corr=1. corr1=1. ; lf=4. lf=2./3. scale=(2000./450.) ;scale=1. ;scaling from ion to e- preacceleration ; corr=1.0*(1.+.32*(1.-alog(lf*energy*scale)/8.9)*exp(-lf*energy*scale/350.) + .05*(1.-alog(lf*energy*scale)/12.)) ;due to leakage fields into ESA ; corr1= (1.+.10*((1.-alog(energy*scale/.05)/7.8)>0.)) ;due to ESA exit grid en_eff = en_eff*corr*corr1 ;*********************************************************************************** amode=thm_read_esa_angle_full_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] ; this needs to be checked 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, magf, phi_offset - are to be added later efd_dat = {project_name: 'THEMIS', $ spacecraft: sc, $ data_name: 'EESA 3D Full', $ apid: '457'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} if sc eq 'a' then begin common tha_457,tha_457_ind,tha_457_dat & tha_457_dat=efd_dat & tha_457_ind=0l endif else if sc eq 'b' then begin common thb_457,thb_457_ind,thb_457_dat & thb_457_dat=efd_dat & thb_457_ind=0l endif else if sc eq 'c' then begin common thc_457,thc_457_ind,thc_457_dat & thc_457_dat=efd_dat & thc_457_ind=0l endif else if sc eq 'd' then begin common thd_457,thd_457_ind,thd_457_dat & thd_457_dat=efd_dat & thd_457_ind=0l endif else if sc eq 'e' then begin common the_457,the_457_ind,the_457_dat & the_457_dat=efd_dat & the_457_ind=0l endif else if sc eq 'f' then begin common thf_457,thf_457_ind,thf_457_dat & thf_457_dat=efd_dat & thf_457_ind=0l endif ;store time range & eclipse correction status eclipse = undefined(use_eclipse_corrections) ? 0:use_eclipse_corrections thm_part_trange, sc, 'peef', set={trange:timerange(trange),eclipse:eclipse} ; help,efd_dat,/st dprint, dlevel=4, 'Data loading complete: th'+sc+' peef - apid 457' ;*********************************************************************************** ; Make tplot structures cols=get_colors() en_dat=fltarr(ntot,32) energy=fltarr(ntot,32) for i=0l,ntot-1l do begin ; print,ntot,i,cs_ptr[i],cs_ind[i] case cs_ind[i] of 0: en_dat[i,0:31]=total(reform(decompress[efd_dat.dat0[cs_ptr[i],*]],88,32),1) 1: en_dat[i,0:14]=total(reform(decompress[efd_dat.dat1[cs_ptr[i],*]],88,15),1) endcase energy[i,*]=reform(efd_dat.energy[*,en_ind[i]]) endfor if not keyword_set(gap_time) then gap_time=1500. 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+'_peef_en_counts'+suffix store_data,name1,data={x:(efd_dat.time+efd_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+'_peef_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