thm_load_peir.pro (Documentation for /home/davin/idl/socware/)

;+
;PROCEDURE:	thm_load_peir
;PURPOSE:	
;	Decommutates the raw THEMIS ion esa reduced data packets, (peir - apid 455), and puts data structure in a common for access by get_th*_peir.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_peir,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_455,tha_455_ind,tha_455_dat & tha_455_dat=0 & tha_455_ind=-1l
	endif else if sc eq 'b' then begin
		common thb_455,thb_455_ind,thb_455_dat & thb_455_dat=0 & thb_455_ind=-1l
	endif else if sc eq 'c' then begin
		common thc_455,thc_455_ind,thc_455_dat & thc_455_dat=0 & thc_455_ind=-1l
	endif else if sc eq 'd' then begin
		common thd_455,thd_455_ind,thd_455_dat & thd_455_dat=0 & thd_455_ind=-1l
	endif else if sc eq 'e' then begin
		common the_455,the_455_ind,the_455_dat & the_455_dat=0 & the_455_ind=-1l
	endif else if sc eq 'f' then begin
		common thf_455,thf_455_ind,thf_455_dat & thf_455_dat=0 & thf_455_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_455_'+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=9				; number of modes
	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			; 50A x 24E, snapshot, 1 spin
	mode_decode[3,2]=3			;  1A x 16E, snapshot, 1 spin
	mode_decode[4,2]=4			; 72A x 16E, snapshot, 1 spin
	mode_decode[1,3]=5			;  1A x 32E, snapshot, 1 spin
	mode_decode[5,3]=6			; 50A x 24E, snapshot, 1 spin
	mode_decode[6,1]=7			;  6A x 32E, snapshot, 1 spin
	mode_decode[7,2]=8			;  6A x 16E, snapshot, 1 spin

; define mode variables for different modes
	nspins = [1,1,1,1,1,1,1,1,1]					; # of spins between measurements in mode
	nenergy = [16,32,24,16,16,32,24,32,16]				; # of energies in mode
	nangle = [1,1,50,1,72,1,50,6,6]					; # of angles in mode
	dat_len = nenergy*nangle					; size data arrays
	spin_decode  = [160,96,3,192,3,96,3,20,38]			; # measurements in packet
	case_decode  = [0,1,5,0,4,1,5,3,2]				; datl[16,32,96,192,1152,1200]==>size[0,1,2,3,4,5]
	angle_decode = [0,1,2,3,5,1,4,6,7]				; angle mode index
	energy_decode = [0,1,2,3,3,4,5,1,3]				; 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) 
	dat_4=bytarr(nmax,1152) 
	dat_5=bytarr(nmax,1200) 
	n_case=lonarr(6) & n_case[*]=0
	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, '455 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.
	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]
	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
	config = [adat[pntr+hdrl],adat[pntr+hdrl+1]]
;	print,'Config bytes = ',config
	mode=mode_decode[config[0],config[1]]
	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]]
	if pktl eq pktl2 then begin
		if mode eq mode2 and mode ne -1 then begin
			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 else last_spin_period=spin_period
		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 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
	pktl=0l
	npkt=0l
	ntot=0l
	hdr  = adat[pntr:pntr+hdrl-1]
	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.




while pntr+1 lt data_size do begin
	hdr  = adat[pntr:pntr+hdrl-1]
	apid_cntr= hdr[3]
	spin_num = hdr[15]
	
	time = (1.d*((hdr[6]*256ul+hdr[7])*256ul+hdr[8])*256ul+hdr[9])+(hdr[10]+hdr[11]/256.d)/256.
	time2 = time
	
;	if spindata then begin
;		if n_elements(spinmodel.s_time) eq 1 then begin
;			spin_offset=long((time-spinmodel.s_time)/spinmodel.spin_period+.5d)
;			time2=spinmodel.s_time+spin_offset*spinmodel.spin_period
;			if time2 gt time then time2=time2-spinmodel.spin_period
;			spin_period=spinmodel.spin_period
;;			offset=2.37+(6./32.)*(spin_period-3.)
;;			if apid_cntr ne last_apid_cntr then   print,'ERROR: peir 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: peir 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
;			mintm=min(abs(spinmodel.s_time-time),ind)
;			if spinmodel.s_time(ind) gt time and ind ne 0 then ind=ind-1
;			spin_offset=long((time-spinmodel.s_time(ind))/spinmodel.spin_period(ind)+.5d)
;			time2=spinmodel.s_time(ind)+spin_offset*spinmodel.spin_period(ind)
;			if time2 gt time then time2=time2-spinmodel.spin_period(ind)
;			spin_period=spinmodel.spin_period(ind)
;;			offset=2.37+(6./32.)*(spin_period-3.)
;;			if apid_cntr ne last_apid_cntr then   print,'ERROR: peir 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: peir 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



	while pktp+4 lt pktl do begin

		config_1 = pkt[pktp] & config_2 = pkt[pktp+1]
		mode = mode_decode[config_1,config_2]
		if mode ne -1 then datl = dat_len[mode]
		if mode ne -1 and pktp+datl+4 le pktl then begin
;			datl = dat_len[mode]
			navg = pkt[pktp+3]
			nsp = nspins[mode]
			config1[ntot] = config_1
			config2[ntot] = config_2
			an_ind[ntot] = angle_decode[mode]
			en_ind[ntot] = energy_decode[mode]
			md_ind[ntot] = mode
			nc = case_decode[mode]
			cs_ind[ntot] = nc
			cs_ptr[ntot] = n_case[nc]

			if mode ne last_mode and pntr+pktl+hdrl+hdrl-1 lt data_size and mode_change ne 1 then begin
				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
				next_time = (1.d*((hdr[6]*256ul+hdr[7])*256ul+hdr[8])*256ul+hdr[9])+(hdr[10]+hdr[11]/256.d)/256.
				n_rem = (pktl - pktp)/(datl+4)
				if n_rem*(datl+4) ne (pktl-pktp) then n_rem=n_rem+1 
				rem_spin = fix((next_time-time)/spin_period+.5) - n_rem*nsp
;				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)
				sub_spin=rem_spin>sub_spin
;				print,'sub_spin = ',sub_spin
			endif

			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
				4: dat_4[n_case[nc],*]=tmp
				5: dat_5[n_case[nc],*]=tmp
			endcase
			ntot = ntot + 1
			n_case[nc] = n_case[nc] + 1

		endif else begin
;			print,'Error: Invalid ESA config header bytes - skipping packet'
;			print,mode,config_1,config_2,pntr,pktl,pktp,datl,pktp+datl+4
			datl=0 & nsp=1
			pktp = pktl
		endelse

		pktp = pktp + 4 + datl
		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
	pntr = pntr + pktl + hdrl
;	print,pntr,'  ',ntot
	last_time=time

endwhile

	tm2=systime(1)
	dprint, dlevel=4, '455 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,*]			
	if n_case[4] gt 0 then dat_4=dat_4[0:n_case[4]-1,*] else dat_4=dat_4[0,*]			
	if n_case[5] gt 0 then dat_5=dat_5[0:n_case[5]-1,*] else dat_5=dat_5[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.i_energy
denergy=emode.i_denergy
nenergy=emode.i_nenergy

;***********************************************************************************
; Calibration: geometric factor and deadtime and efficiency

calib=get_thm_esa_cal(sc=sc,ion=1,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

deadtime=1.7e-7

; Make ion energy efficiency array 
; 	REF: Funsten curves for CEMs, assumes 2 kV pre-acceleration

en_eff=(1.+0.3*((energy+2000.)/30000.)^.5)/1.075

; empirical correction from derived from sheath ion/electron density comparison
; 	was used until 08-07-11: corr=0.85+.35*exp(-energy/400.) & en_eff=en_eff*corr
  corr=1. 						;default
  corr1=1. & corr2=1.					;default
;  corr=0.85+.35*exp(-energy/400.)			;due to leakage fields into ESA
;  corr=0.95+.95*exp(-energy/40.)			;due to leakage fields into ESA
;  corr=.905*(1.+.55*exp(-energy/180.))			;due to leakage fields into ESA
;  corr=.91*(1.+.41*((1.-alog(energy)/7.)>0.))		;due to leakage fields into ESA
  lf=2./3.
;  lf=1.
  corr=.91*(1.+.32*(1.-alog(lf*energy)/8.9)*exp(-lf*energy/350.) + .05*(1.-alog(lf*energy)/12.))	;due to leakage fields into ESA
  corr1=   (1.+.10*((1.-alog(energy/.25)/7.8)>0.))	;due to ESA exit grid 
  corr2=   (1.+.05*((1.-alog(energy/1.5)/7.8)>0.))	;due to MCP grid 


en_eff=en_eff*corr*corr1*corr2

;***********************************************************************************

amode=thm_read_iesa_angle_reduced_mode() 

nbins=amode.i_nbins
theta=amode.i_theta
dtheta=amode.i_dtheta
phi=amode.i_phi
dphi=amode.i_dphi
domega=amode.i_domega
nrg_wt=amode.i_nrg_wt
anodes=amode.i_anodes
an_map=amode.i_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_pei_bkg.pro
;   bkg, bkg_pse, bkg_pei 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_pei is an iESA determined background
;   	bkg is generally the smaller of bkg_pse, bkg_pei -- see keywords in thm_load_pei_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_pei=fltarr(ntot) 	& bkg_pei[*]=!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


ird_dat = 	{project_name:		'THEMIS',				$
		spacecraft:		sc, 					$
		data_name:		'IESA 3D Reduced', 			$
		apid:			'455'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,					$
		geom_factor: 		geom_factor,				$
		dead: 			deadtime,				$

		mass: 			1.0438871e-02, 				$
		charge: 		1., 					$
		sc_pot: 		replicate(0.,ntot),			$

		magf:	 		magf, 					$

		bkg_pse:	 	bkg_pse, 				$
		bkg_pei:	 	bkg_pei, 				$
		bkg:	 		bkg, 					$
		bkg_arr:		bkg_arr,				$

		dat0:			dat_0,					$
		dat1:			dat_1,					$
		dat2:			dat_2,					$
		dat3:			dat_3,					$
		dat4:			dat_4,					$
		dat5:			dat_5}


if sc eq 'a' then begin
	common tha_455,tha_455_ind,tha_455_dat & tha_455_dat=ird_dat & tha_455_ind=0l
endif else if sc eq 'b' then begin
	common thb_455,thb_455_ind,thb_455_dat & thb_455_dat=ird_dat & thb_455_ind=0l
endif else if sc eq 'c' then begin
	common thc_455,thc_455_ind,thc_455_dat & thc_455_dat=ird_dat & thc_455_ind=0l
endif else if sc eq 'd' then begin
	common thd_455,thd_455_ind,thd_455_dat & thd_455_dat=ird_dat & thd_455_ind=0l
endif else if sc eq 'e' then begin
	common the_455,the_455_ind,the_455_dat & the_455_dat=ird_dat & the_455_ind=0l
endif else if sc eq 'f' then begin
	common thf_455,thf_455_ind,thf_455_dat & thf_455_dat=ird_dat & thf_455_ind=0l
endif


;store time range & eclipse correction status
eclipse = undefined(use_eclipse_corrections) ? 0:use_eclipse_corrections
thm_part_trange, sc, 'peir', set={trange:timerange(trange),eclipse:eclipse}


; help,ird_dat,/st
dprint, dlevel=4, 'Data loading complete: th'+sc+' peir - apid 455'

;***********************************************************************************
; 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:15]=reform(decompress[ird_dat.dat0[cs_ptr[i],*]])
		1: en_dat[i,0:31]=reform(decompress[ird_dat.dat1[cs_ptr[i],*]])
		2: en_dat[i,0:15]=total(reform(decompress[ird_dat.dat2[cs_ptr[i],*]],6,16),1)
		3: en_dat[i,0:31]=total(reform(decompress[ird_dat.dat3[cs_ptr[i],*]],6,32),1)
		4: en_dat[i,0:15]=total(reform(decompress[ird_dat.dat4[cs_ptr[i],*]],72,16),1)
		5: en_dat[i,0:23]=total(reform(decompress[ird_dat.dat5[cs_ptr[i],*]],50,24),1)
	endcase
	energy(i,*)=reform(ird_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+'_peir_en_counts'+suffix
store_data,name1,data={x:(ird_dat.time+ird_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','i+ 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+'_peir_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','i+ 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