;+

;PROCEDURE:	get_sp_spec

;PURPOSE:	

;	Generates spin phase angle vs. time spectrogram data structures for tplot

;INPUT:		

;	data_str, 	a string (either 'c1_hia_p6','c3_hia_p6', ...)

;			where get_'string' returns a 2D  

;			pitch angle data structure

;KEYWORDS:		

;	T1:		start time, seconds since 1970

;	T2:		end time, seconds since 1970		

;	ENERGY:		fltarr(2)		energy range to sum over

;	ERANGE:		intarr(2)		energy bin range to sum over

;	EBINS:		bytarr(dat.nenergy)	energy 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

;	BKG:  		A 3d data structure containing the background counts.

;	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_"get_dat"_ts

;	EWEIGHT		Weights the energy bin average by bin energy width.

;	SHIFT90		If set true, pitch angle range is -90 to 270

;

;CREATED BY:	J.McFadden

;VERSION:	1

;LAST MODIFICATION:  97/07/11

;MOD HISTORY:

;		96/08/14  	Added EWEIGHT keyword  - mcfadden

;		97/03/04	T1,T2 keywords added

;		97/07/11	CALIB keyword added

;		97/07/23	fixed bug - ind1=fix((max(ind)+min(ind))/2) 

;

;NOTES:	  

;	Current version only works for FAST

;-



pro get_sp_spec,data_str,  $

	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, $

	eweight = eweight, $

	shift90 = shift90



;	Time how long the routine takes

ex_start = systime(1)



;	Set defaults for keywords, etc.



n = 0

max = 30000        ; this could be improved

all_same = 1



routine = 'get_'+data_str



;t = 0             ; get first sample

;dat = call_function(routine,t,CALIB=calib,index=0)

; FAST uses "start" keyword



if keyword_set(t1) then begin

	t=t1

	if routine eq 'get_fa_sebs' then dat = call_function(routine,t,/first) else dat = call_function(routine,t,CALIB=calib)

endif else begin

	t = 1000             ; get first sample

	dat = call_function(routine,t,CALIB=calib, /start)

endelse



if dat.valid eq 0 then begin no_data = 1 & return & end $

else no_data = 0



ytitle = data_str + '_pa_spec'

last_time = (dat.time+dat.end_time)/2.

nenergy = dat.nenergy

nmaxvar = dat.nbins



default_gap_time = 200

if dat.project_name eq 'FAST' then begin

	nmaxvar=64

	default_gap_time = 8.

endif

if not keyword_set(gap_time) then gap_time = default_gap_time



time   = dblarr(max)

data   = fltarr(max,nmaxvar)

var   = fltarr(max,nmaxvar)

nvar = dat.nbins

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 max) then  goto, continue  ; goto YUCK!



if keyword_set(t2) then tmax=t2 else tmax=1.e30



while (dat.valid ne 0) and (n lt max) do begin

if (dat.valid eq 1) then begin



	count = dat.nenergy



	if keyword_set(en) then begin

		er=energy_to_ebin(dat,en)

	endif

	if keyword_set(er) then begin

		ebins=bytarr(dat.nenergy)

		if er(0) gt er(1) then er=reverse(er)

		ebins(er(0):er(1))=1

	endif

	if not keyword_set(ebins) then begin

		ebins=bytarr(dat.nenergy)

		ebins(*)=1

	endif



	if keyword_set(retrace) then ebins(0:retrace-1)=0

	ind=where(ebins,count)



	if units eq 'Counts' then norm = 1 else norm = count

	if ind(0) ne -1 and keyword_set(eweight) then norm = total(dat.denergy(ind,0),1)



	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 dat = sub3d(dat,bkg)

	if keyword_set(units) then dat = conv_units(dat,units)



	nvar = dat.nbins

	if nvar gt nmax then nmax = nvar

	time(n)   = (dat.time+dat.end_time)/2.

	if ind(0) eq -1 then begin

		tmpdata = total( dat.data(0,*), 1)

		tmpdata(*) = 0

		tmpvar = total( dat.theta(0,*), 1)

		if keyword_set(shift90) then tmpvar = $

		((360.*(tmpvar/360.-floor(tmpvar/360.)) + 90.) mod 360.) -90.

	endif else begin

		if keyword_set(eweight) then begin

			tmpdata = total( dat.data(ind,*)*dat.denergy(ind,*), 1)/norm

		endif else begin

			tmpdata = total( dat.data(ind,*), 1)/norm

		endelse

; The following line will need to change if theta depends upon the energy index

;  Some of this section could be cleaned up 

		ind1= fix((max(ind)+min(ind))/2)

		tmpvar = total( dat.theta(ind1,*), 1)

		if keyword_set(shift90) then tmpvar = $

		((360.*(tmpvar/360.-floor(tmpvar/360.)) + 90.) mod 360.) -90.

	endelse

; Shift the angle array so that angles increase monotonicly with index 

; -- needed for tplot

	minvar = min(tmpvar,indminvar)

	if (indminvar gt 1) then begin

		if tmpvar(0) gt tmpvar(1) then begin

			tmpvar=reverse(shift(tmpvar,-indminvar-1))

			tmpdata=reverse(shift(tmpdata,-indminvar-1))

		endif else begin

			tmpvar=shift(tmpvar,-indminvar)

			tmpdata=shift(tmpdata,-indminvar)

		endelse

	endif else begin

		if tmpvar(2) gt tmpvar(3) then begin

			tmpvar=reverse(shift(tmpvar,-indminvar-1))

			tmpdata=reverse(shift(tmpdata,-indminvar-1))

		endif else begin

			tmpvar=shift(tmpvar,-indminvar)

			tmpdata=shift(tmpdata,-indminvar)

		endelse

	endelse

	var(n,0:nvar-1) = tmpvar

	data(n,0:nvar-1) = tmpdata



; 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



	dat = call_function(routine,t,CALIB=calib,/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 nenergy then ytitle = ytitle+'_'+strtrim(count,2)

	if keyword_set(name) eq 0 then name=ytitle else ytitle = name

	ytitle = ytitle+' ('+units+')'



data = data(0:n-1,0:nmax-1)

var = var(0:n-1,0:nmax-1)















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)+ ' deg'

;endif



time = time(0:n-1)



if keyword_set(t1) then begin

	ind=where(time ge t1)

	time=time(ind)

	data=data(ind,*)

	var=var(ind,*)

endif

if keyword_set(t2) then begin

	ind=where(time le t2)

	time=time(ind)

	data=data(ind,*)

	var=var(ind,*)

endif



;datastr = {ztitle:units,x:time,y:data,v:var,  $

;	labels:labels,	$

;    ylog:0,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

