;+
;PROCEDURE:	thm_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
;		09/04/16	bkg keyword removed for Themis electron sst, assume bkg removal is part of get_th?_pse?
; 2014/06/30  Adding thm prefix to differentiate from identically named FAST routines
; 2016/06/30  re-enabling time range constraints, minor changes for spedas compatibility
;
;NOTES:	  
;	Current version only works for FAST and THEMIS
;-

pro thm_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, $
	_extra=_extra

  compile_opt strictarr

;	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

	if ndimen(dat.data) eq ndimen(dat.bins) then dat.data=dat.data*dat.bins
	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' and dat.charge gt 0. then dat=thm_pei_bkg_sub(dat) $ 
		else if dat.UNITS_PROCEDURE eq 'thm_convert_esa_units' and dat.charge lt 0. then dat=thm_pee_bkg_sub(dat) 
;		else if strmid(dat.DATA_NAME,0,12) eq 'SST Electron' then dat=thm_pse_bkg_sub(dat) ; added to get_routine
	endif
	if keyword_set(units) then dat = conv_units(dat,units,_extra=_extra)

	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
	dprint, 'Invalid packet, dat.valid ne 1, at: ',time_string(dat.time), dlevel=2
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:


;make sure valid data was found -af
if undefined(count) then begin
  no_data = 1
  dprint, 'No valid data found for "'+get_dat+'"', dlevel=1
  return
endif



;	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 dprint, 'all_same=',strtrim(all_same,2),dlevel=4
;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]

dprint, 'time range = ',time_string(minmax(time)), dlevel=4

if keyword_set(t1) then begin
	ind=where(time ge t1,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)
	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
dprint, strtrim(ex_time,2)+' seconds execution time.', dlevel=2
dprint, 'Number of data points = ',strtrim(n,2), dlevel=4

return

end