;+
;PROCEDURE:	thm_load_peeb
;PURPOSE:	
;	Decommutates the raw THEMIS electron esa burst data packets, (peeb - apid 459), and puts data structure in a common for access by get_th*_peeb.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
;
;
;CREATED BY:	J. McFadden
;VERSION:	1
;LAST MODIFICATION:  07/03/18
;MOD HISTORY:
;
;NOTES:	  
;	
;-

pro thm_load_peeb,file=file,sc=sc,themishome=themishome,gap_time=gap_time

; sc default
	if not keyword_set(sc) then begin
		dprint,dlevel=3,'S/C number not set, default = a'
		sc='a'
	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()
		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.))
		dprint,dlevel=4,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_459_'+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,dlevel=3,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,dlevel=3,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
	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





; define mode variables for different modes
	nspins = [1,1,1,1,1]				; # of spins between measurements in mode
	nenergy = [32,32,32,32,32]				; # of energies in mode
	nangle = [88,88,88,88,88]				; # of angles in mode
	dat_len = nenergy*nangle				; size data arrays
	spin_decode  = [1,1,1,1,1]				; # measurements in packet
	case_decode  = [0,0,0,0,0]				; datl[16,32,96,192,1152,1200]==>size[0,1,2,3,4,5]
	angle_decode = [0,0,0,0,0]				; angle mode index
	energy_decode = [0,1,1,2,2]				; energy mode index

; initialize arrays
	ndays=n_elements(file)
	nmax=40000l*ndays
	dat_0=bytarr(nmax/4,2816) 





	n_case=lonarr(1) & 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=3,'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=3,'Data download complete'
	dprint,dlevel=3,'459 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
;	dprint,dlevel=3,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]]
;	dprint,dlevel=3,'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
;	dprint,dlevel=3,pntr,'  ',pktl,'  ',time_string(time),'  ',adat[pntr+hdrl],'  ',adat[pntr+hdrl+1],'  spin period= ',spin_period
endwhile

;	if tf_spin then dprint,dlevel=0,'ERROR : spin period calculation error, setting spin period to 3. seconds'
	if tf_spin then spin_period=3.
	dprint,dlevel=3,'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)

; Start decommutation

;	help,adat
	dprint,dlevel=3,'Decommutating data ...  Data volume (bytes) =',data_size
;	dprint,dlevel=3,'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=1.83+(6./32.)*(spin_period-3.)
;			if apid_cntr ne last_apid_cntr then   print,'ERROR: peeb 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: peeb 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=1.83+(6./32.)*(spin_period-3.)
;			if apid_cntr ne last_apid_cntr then   print,'ERROR: peeb 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: peeb 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
	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)



	while pktp+datl+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 and mode eq last_mode 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]

















			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





			endcase
			ntot = ntot + 1
			n_case[nc] = n_case[nc] + 1

		endif else begin
			pktp = pktl
			dprint,dlevel=3,'Error: Invalid ESA config header bytes or mode change in packet - skipping packet'
			dprint,dlevel=3,mode,config_1,config_2,pntr,pktl,pktp
			datl=0 & nsp=1
		endelse

		pktp = pktp + 4 + datl

		sub_spin = sub_spin + nsp
		last_mode = mode
;		if ntot-1000l*ceil(ntot/1000l) eq 0 then dprint,dlevel=3,pntr,'  ',pktl,'  ',time_string(time),'  ',adat[pntr+hdrl],'  ',adat[pntr+hdrl+1]
	endwhile
	pntr = pntr + pktl + hdrl
;	dprint,dlevel=3,pntr,'  ',ntot
	last_time=time

endwhile

	tm2=systime(1)
	dprint,dlevel=3,'459 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,*]






; 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_burst_mode(fm)

energy=emode.e_energy
denergy=emode.e_denergy
nenergy=emode.e_nenergy

; 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 = 600. ;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
eff = (1-exp(-k*delta/delta_max))/(1-exp(-k))

; the following needs some work
; corr=(1.0-.23*exp(-(energy/60.)^2)+.05*exp(-(energy/40.)^2))
; eff=eff*corr*1.135*.94

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

amode=thm_read_esa_angle_burst_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

;***********************************************************************************
; to get magnetometer data, use thm_load_esa_mag.pro

	magf=fltarr(ntot,3)
	b_gse=fltarr(ntot,3)

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

;calib=thm_read_esa_burst_calib()

;gf=calib.gf
;geom_factor=calib.geom_fact

;***********************************************************************************
; Make 3D data structures
; Note: sc_pot, b_gse, magf, rot, phi_offset - are to be added later

	dt_arr=(dphi/11.25)*(dtheta/22.5)*anodes*nrg_wt
	gf=dphi/22.5*dtheta/22.5*nrg_wt
	geom_factor=.0011						; geom_factor of 22.5x22.5 angle bin
	dt=e_time-s_time

ebd_dat = 	{project_name:		'THEMIS',				$
		spacecraft:		sc, 					$
		data_name:		'EESA 3D BURST', 			$
		apid:			459,					$
		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: 			eff,	 				$

		nbins: 			nbins, 					$
		theta: 			theta,  				$
		dtheta: 		dtheta,  				$
		phi: 			phi,  					$
		dphi: 			dphi,  					$
		phi_offset: 		replicate(0.,ntot),  			$
		domega: 		domega,  				$
		gf: 			gf, 					$

		geom_factor: 		geom_factor,				$
		dead: 			1.e-7,					$
		mass: 			5.68566e-06, 				$
		charge: 		-1., 					$
		sc_pot: 		replicate(0.,ntot),			$

		b_gse:	 		b_gse, 					$
		magf:	 		magf, 					$
;		rot:	 		replicate(1.,ntot)#[[1.,0.,0.],[0.,1.,0.],[0.,0.,1.]],		$

		dat0:			dat_0}







if sc eq 'a' then begin
	common tha_459,tha_459_ind,tha_459_dat & tha_459_dat=ebd_dat & tha_459_ind=0l
endif else if sc eq 'b' then begin
	common thb_459,thb_459_ind,thb_459_dat & thb_459_dat=ebd_dat & thb_459_ind=0l
endif else if sc eq 'c' then begin
	common thc_459,thc_459_ind,thc_459_dat & thc_459_dat=ebd_dat & thc_459_ind=0l
endif else if sc eq 'd' then begin
	common thd_459,thd_459_ind,thd_459_dat & thd_459_dat=ebd_dat & thd_459_ind=0l
endif else if sc eq 'e' then begin
	common the_459,the_459_ind,the_459_dat & the_459_dat=ebd_dat & the_459_ind=0l
endif else if sc eq 'f' then begin
	common thf_459,thf_459_ind,thf_459_dat & thf_459_dat=ebd_dat & thf_459_ind=0l
endif

; help,ebd_dat,/st
dprint,dlevel=3,'Data loading complete: th'+sc+' peeb - apid 459'

;***********************************************************************************
; Make tplot structures

cols=get_colors()

en_dat=fltarr(ntot,32)
energy=fltarr(ntot,32)

for i=0l,ntot-1l do begin
;	dprint,dlevel=3,ntot,i,cs_ptr[i],cs_ind[i]
	case cs_ind[i] of
		0: en_dat[i,0:31]=total(reform(decompress(ebd_dat.dat0[cs_ptr[i],*]),88,32),1)





	endcase
	energy[i,*]=reform(ebd_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+'_peeb_en_counts'
store_data,name1,data={x:(ebd_dat.time+ebd_dat.end_time)/2.,y:en_dat,v:energy}
	zlim,name1,0,0,1
	ylim,name1,3.,40000.,1
	options,name1,'ztitle','Counts'
	options,name1,'ytitle','e- th'+sc+'!C!C eV'
	options,name1,'spec',1
	options,name1,'x_no_interp',1
	options,name1,'y_no_interp',1

name1='th'+sc+'_peeb_mode'
store_data,name1,data={x:(s_time+e_time)/2.,y:[[config1-.1],[config2+.1],[md_ind]]}
	ylim,name1,-1,10,0
	options,name1,'ytitle','i+ th'+sc+'!C!C mode'
	options,name1,'colors',[cols.red,cols.green,cols.blue]
	options,name1,labels=['config1', 'config2', 'mode'],constant=0.
	options,name1,'labflag',3
	options,name1,'labpos',[0,3,6]

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

end