;+
;PROCEDURE:	get_en_spec
;PURPOSE:	
;	Generates energy-time spectrogram data structures for tplot
;INPUT:		
;	get_dat, 	a string that defines the 'routine' that returns a 2D or 3D data structure
;			if 'probe' keyword set, assumes routine = get_dat
;			if 'probe' keyword is not set, assumes routine = 'get_'+ get_dat
;				unless get_dat already includes 'get_', in which case routine = get_dat
;			for example, get_dat='thc_peir' then routine = 'get_thc_peir'
;			for example, get_dat='thm_sst_pser' then routine = 'thm_sst_pser'
;
;KEYWORDS:
;	T1:		start time, seconds since 1970
;	T2:		end time, seconds since 1970		
;	ANGLE:		fltarr(2),fltarr(4)	angle range to sum over
;	ARANGE:		intarr(2)		bin range to sum over
;	BINS:		bytarr(dat.nbins)	bins to sum over
;	gap_time: 	time gap big enough to signify a data gap 
;			(default 200 sec, 8 sec for FAST)
;	NO_DATA: 	returns 1 if no_data else returns 0
;	UNITS:		convert to these units if included
;	NAME:  		New name of the Data Quantity, default = get_dat + '_en_spec'
;	BKG:  		A 3d data structure containing the background counts 
;			If bkg=1, it will run themis background subtraction routine
;	FLOOR:  	Sets the minimum value of any data point to sqrt(bkg).
;	MISSING: 	value for bad data.
;	RETRACE: 	Set to number of retrace energy steps to be eliminated starting at energy step 0
;	CALIB:		Calib keyword passed on to get_dat -- no longer used
;	PROBE:		probe keyword passed on to get_dat -- used by themis sst routines
;
;
;
;CREATED BY:	J.McFadden
;VERSION:	1
;LAST MODIFICATION:  97/03/04
;MOD HISTORY:
;		97/03/04	T1,T2 keywords added
;		97/05/22	CALIB keyword added
;		08/07/10	bkg keyword changed to work for Themis ion ESA
;		08/07/10	averaging routine changed to weight different angle bins by dat.gf 
;		08/12/31	bkg keyword changed to work for Themis electron sst
;		08/12/31	probe keyword added for themis sst routines
;		09/01/05	changed loop to use index to allow for dat.valid=0 data in stored data
;		09/01/05	calib keyword no longer passed to routine
;
;
;NOTES:	  
;	Current version only works for FAST and THEMIS
;-

pro get_en_spec,get_dat,  $
	T1=t1, $
	T2=t2, $
;	ENERGY=en, $
;	ERANGE=er, $
;	EBINS=ebins, $
	ANGLE=an, $
	ARANGE=ar, $
	BINS=bins, $
	gap_time=gap_time, $ 
	no_data=no_data, $
	units = units,  $
        name  = name, $
	bkg = bkg, $
        missing = missing, $
        floor = floor, $
        retrace = retrace, $
        CALIB = calib, $
	probe = probe



;	Time how long the routine takes
ex_start = systime(1)

;	Set defaults for keywords, etc.

n = 0l
all_same = 1


if strlowcase(strmid(get_dat,0,4)) eq 'get_' then begin
	routine = get_dat
	get_dat = strlowcase(strmid(get_dat,4,strlen(get_dat)))
endif else if keyword_set(probe) then begin
	routine = get_dat
endif else routine = 'get_'+get_dat 

if keyword_set(probe) then times = call_function(routine,probe=probe,/times) else times = call_function(routine,/times)
maxind = n_elements(times)-1

if keyword_set(t1) then tmpmax=min(abs(times-time_double(t1)),idx) else idx=1l
if keyword_set(t2) then tmpmax=min(abs(times-time_double(t2)),idxmax) else idxmax=maxind

if keyword_set(probe) then begin
	dat = call_function(routine,probe=probe,index=idx) 
	while (dat.valid eq 0 and idx lt idxmax) do begin
		idx = idx + 1
		dat = call_function(routine,probe=probe,index=idx)
	endwhile
endif else begin
	dat = call_function(routine,index=idx)
	while (dat.valid eq 0 and idx lt idxmax) do begin
		idx = idx + 1
		dat = call_function(routine,index=idx)
	endwhile
endelse

;if keyword_set(t1) then begin
;	t=t1
;	if keyword_set(calib) then dat = call_function(routine,t,CALIB=calib) $
;		else if keyword_set(probe) then dat = call_function(routine,t,probe=probe,index=idx) $
;		else if routine eq 'get_fa_sebs' then dat = call_function(routine,t,/first) $
;		else dat = call_function(routine,t)
;endif else begin
;	t = 1000             ; get first sample
;	idx=0
;	if keyword_set(calib) then dat = call_function(routine,t,CALIB=calib,/start) $
;		else if keyword_set(probe) then dat = call_function(routine,probe=probe,index=idx) $
;		else if routine eq 'get_fa_sebs' then dat = call_function(routine,t,/first) $
;		else dat = call_function(routine,t,/start)
;endelse
;if dat.valid eq 0 then begin no_data = 1 & return & end $
;else no_data = 0


ytitle = get_dat + '_en_spec'
last_time = (dat.time+dat.end_time)/2.
nbins = dat.nbins
nmaxvar = dat.nenergy

default_gap_time = 200.
if dat.project_name eq 'FAST' then begin
	nmaxvar=96
	default_gap_time = 8.
endif
if strupcase(strmid(dat.project_name,0,6)) eq 'THEMIS' then begin
	nmaxvar=32
	default_gap_time = 13.
endif
if not keyword_set(gap_time) then gap_time = default_gap_time

maxpts = idxmax-idx+1000 < 300000l			; this limits a tplot str to < 77 MBytes
time   = dblarr(maxpts)
data   = fltarr(maxpts,nmaxvar)
var   = fltarr(maxpts,nmaxvar)
nvar = dat.nenergy
nmax=nvar

if not keyword_set(units) then units = 'Counts'
if not keyword_set(missing) then missing = !values.f_nan



;	Collect the data - Main Loop
    
; May want to use the following lines when "index" is operational in FAST get* routines    
;times = call_function(routine,t,CALIB=calib, /times)
;for idx=0,n_elements(times)-1 do begin
;    if (dat.valid eq 0) or (n ge maxpts) then  goto, continue  ; goto YUCK!


;while (dat.valid ne 0) and (n lt maxpts) do begin
while (idx lt idxmax) and (n lt maxpts) do begin

if (dat.valid eq 1) then begin

	count = dat.nbins
	if keyword_set(an) then begin
		str_element,dat,'PHI',INDEX=tf_phi
		if tf_phi lt 0 then bins=angle_to_bins(dat,an)
		if tf_phi gt 0 then begin
			th=reform(dat.theta(0,*)/!radeg)
			ph=reform(dat.phi(fix(dat.nenergy/2),*)/!radeg)
			xx=cos(ph)*cos(th)
			yy=sin(ph)*cos(th)
			zz=sin(th)
			Bmag=(dat.magf(0)^2+dat.magf(1)^2+dat.magf(2)^2)^.5
			pitch=acos((dat.magf(0)*xx+dat.magf(1)*yy+dat.magf(2)*zz)/Bmag)*!radeg
			if an(0) gt an(1) then an=reverse(an)
			bins= pitch gt an(0) and pitch lt an(1)
			if total(bins) eq 0 then begin
				tmp=min(abs(pitch-(an(0)+an(1))/2.),ind)
				bins(ind)=1
			endif
		endif
	endif
	if keyword_set(ar) then begin
		nb=dat.nbins
		bins=bytarr(nb)
		if ar(0) gt ar(1) then begin
			bins(ar(0):nb-1)=1
			bins(0:ar(1))=1
		endif else begin
			bins(ar(0):ar(1))=1
		endelse
	endif
; Set the "count" to the number of bins summed over

	if not keyword_set(bins) then ind=indgen(dat.nbins) else ind=where(bins,count)
	if max(ind) gt dat.nbins-1 then ind=-1
	if units eq 'Counts' or units eq 'counts' then norm = 1 else norm = count

	if abs((dat.time+dat.end_time)/2.-last_time) ge gap_time then begin
		if n ge 2 then dbadtime = time(n-1) - time(n-2) else dbadtime = gap_time/2.
		time(n) = (last_time) + dbadtime
		data(n,*) = missing
		var(n,*) = missing
		n=n+1
		if (dat.time+dat.end_time)/2. gt time(n-1) + gap_time then begin
			time(n) = (dat.time+dat.end_time)/2. - dbadtime
			data(n,*) = missing
			var(n,*) = missing
			n=n+1
		endif
	endif

	if keyword_set(bkg) then begin
; sub3d requires bkg to be the subtraction array -- sub3d is a general routine and works similar to the pse routine
		if ndimen(bkg) ne 0 then dat = sub3d(dat,bkg) $
;		else if dat.UNITS_PROCEDURE eq 'thm_convert_esa_units' then dat=thm_esa_bgnd_sub(dat) $ ; old routine
		else if dat.UNITS_PROCEDURE eq 'thm_convert_esa_units' then dat=thm_pei_bkg_sub(dat) $
		else if strmid(dat.DATA_NAME,0,12) eq 'SST Electron' then dat=thm_pse_bkg_sub(dat) 
	endif
	if keyword_set(units) then dat = conv_units(dat,units)

	nvar = dat.nenergy
	if nvar gt nmax then nmax = nvar
	time(n)   = (dat.time+dat.end_time)/2.
	if ind(0) ne -1 then begin
		if ndimen(dat.data) gt 1 and n_elements(ind) gt 1 then begin
			if units eq 'Counts' or units eq 'counts' then begin
				data(n,0:nvar-1) = total( dat.data(*,ind), 2)/norm
				var(n,0:nvar-1) = total( dat.energy(*,ind), 2)/count
			endif else begin
				data(n,0:nvar-1) = total( dat.data(*,ind)*dat.gf(*,ind), 2)/total(dat.gf(*,ind),2)
				var(n,0:nvar-1) = total( dat.energy(*,ind), 2)/count
			endelse
		endif else if ndimen(dat.data) gt 1 then begin
			data(n,0:nvar-1) = dat.data(*,ind)/norm
			var(n,0:nvar-1) = dat.energy(*,ind)/count
		endif else begin
			data(n,0:nvar-1) = dat.data(*)
			var(n,0:nvar-1) = dat.energy(*)
		endelse
	endif else begin
		data(n,0:nvar-1) = 0
		var(n,0:nvar-1) =  dat.energy(*,0)
	endelse

; test the following lines, the 96-6-19 version of tplot did not work with !values.f_nan
;	if nvar lt nmaxvar then data(n,nvar:nmaxvar-1) = !values.f_nan
;	if nvar lt nmaxvar then var(n,nvar:nmaxvar-1) = !values.f_nan
	if nvar lt nmaxvar then data(n,nvar:nmaxvar-1) = data(n,nvar-1)
	if nvar lt nmaxvar then var(n,nvar:nmaxvar-1) = 1.5*var(n,nvar-1)-.5*var(n,nvar-2)

	if (all_same eq 1) then begin
		if dimen1(where(var(n,0:nvar-1) ne var(0,0:nvar-1))) gt 1 then all_same = 0
	endif
	last_time = time(n)
	n=n+1

endif else begin
	print,'Invalid packet, dat.valid ne 1, at: ',time_to_str(dat.time)
endelse
	idx=idx+1
	if keyword_set(probe) then dat = call_function(routine,probe=probe,index=idx) else dat = call_function(routine,index=idx)

;	if keyword_set(probe) then idx=idx+1
;	if keyword_set(calib) then dat = call_function(routine,t,CALIB=calib,/ad) $
;		else if keyword_set(probe) then dat = call_function(routine,probe=probe,index=idx) $
;		else dat = call_function(routine,t,/ad)
;;	dat = call_function(routine,t,CALIB=calib,index=idx)
;	if dat.valid ne 0 then if dat.time gt tmax then dat.valid=0

endwhile
;endfor
;continue:






;	Store the data

	if count ne nbins then ytitle = ytitle+'_'+strtrim(count,2)
	if keyword_set(name) eq 0 then name=ytitle else ytitle = name
	ytitle = ytitle+' ('+units+')'

if not keyword_set(retrace) then begin
;	If you want to plot the retrace, set the retrace flag to 1.
	data = data(0l:n-1,0:nmax-1)
	var = var(0l:n-1,0:nmax-1)
endif else begin
	data = data(0l:n-1,retrace:nmax-1)
	var = var(0l:n-1,retrace:nmax-1)
endelse


if not all_same then print,'all_same=',all_same
;labels=''
; The following has be removed so that FAST summary cdf files contain both arrays
;if all_same then begin
;	var = reform(var(0,*))
;	labels = strtrim( round(var) ,2)+ ' eV'
;endif

time = time(0l:n-1)

print,'get_en_spec time range = ',time_string(minmax(time))
;if keyword_set(t1) then begin
;	ind=where(time ge t1,count)
;	print,count
;	if count ne 0 then begin
;		time=time(ind) 
;		data=data(ind,*)
;		var=var(ind,*)
;	endif else return
;endif
;if keyword_set(t2) then begin
;	ind=where(time le t2,count)
;	print,count
;	if count ne 0 then begin
;		time=time(ind)
;		data=data(ind,*)
;		var=var(ind,*)
;	endif else return
;endif

; remove any negative values caused by background subtractions
	data=data>0.

;datastr = {ztitle:units,x:time,y:data,v:var,  $
;	labels:labels,	$
;    ylog:1,panel_size:2.}
datastr = {x:time,y:data,v:var}
store_data,name,data=datastr

ex_time = systime(1) - ex_start
message,string(ex_time)+' seconds execution time.',/cont,/info
print,'Number of data points = ',n

return

end