;+
;  Name:
;       AVERAGE_STEPS
;  Category:
;	BATSE, OCCULTATION
;	
; PURPOSE: average occultation steps using resistant mean and best
;	step model based on standard deviation to fit of step
;	
; CALLING SEQUENCE:
;	average_steps, steps_in, average_step, sigma_cut = sigma_cu
; INPUTS: 
;	steps_in- occultation steps dimensioned (2, nchan, n_durations, n_det, n_obs)
;		or                  dimensioned (2, nchan, n_durations, n_obs)
;		first row is step value and standard deviation from plot_occs.pro
;		nchan - number of energy bins
;		n_durations - number of durations used for models and/or fit order (0,1,2)
;		n_det - number of detectors fit's saved, up to 4 for BATSE
;		n_obs - number of occultations steps fit
;		
;	
; OPTIONAL INPUTS:
;	sigma_cut - sigma_cut used in resistant mean, default 2
;	
; OUTPUTS:
;	Average_step - resultant average
; OPTIONAL OUTPUTS:
;
; PROCEDURE:
;	Preps the occultation results for use with RESISTANT_MEAN.
;
; CALLS: resistant_mean
;
; COMMON BLOCKS:
;
; RESTRICTIONS:
;
; MODIFICATION HISTORY:
;	
;-
;
pro average_steps, steps_in, average_step, sigma_cut = sigma_cut

dims = size( steps_in)
if dims(0) eq 4 then steps = reform( steps_in, dims(1), dims(2), dims(3), 1, dims(4) ) else $
	steps = steps_in

dims = size(steps)
average_step = fltarr( 2, dims(2), dims(4) )


;work with 1 detector at a time
for k= 0, dims(4)-1 do begin ;loop over detectors
checkvar, sigma_cut, 2.

wuse = where( steps(1,0,0,k,*) ne 0, nwuse)

if nwuse ge 1 then begin 
	for m=0,dims(2)-1 do begin ;loop over energy bands
	  steps_mk = reform( steps(*,m,*,k,wuse) )
	  for i=0,nwuse-1 do begin ;select best on std dev for each obs interv.
		steps_mk(0,0,i) = steps_mk(*,sort(steps_mk(1,*,i)),i)
	  endfor
	  resistant_mean, (steps_mk(0,0,*))(*), sigma_cut, mean, sigma ;average over obs ints
	  average_step(0, m, k) = [mean, sigma]
	endfor
endif
endfor
end