;+
; NAME: rbsp_sample_rate.pro
; SYNTAX: 
; PURPOSE: Determines the sample rate of an input time series of data with the
;                 ability to set gap thresholds to avoid including them in the
;                 calculation.
; INPUT: time :  [N]-element array of time series times
; OUTPUT: 
; KEYWORDS: 
;               GAP_THRESH   :  Scalar defining the maximum data gap [s] allowed in
;                                 the calculation
;                                 [Default = MAX(TIME) - MIN(TIME)]
;               AVERAGE      :  If set, routine returns the scalar average sample rate
;                                 [Default = 0, which returns an array of sample rates]
;               OUT_MED_AVG  :  Set to a named variable to return the median and average
;                                 values of the sample rate if the user wants all values
;                                 as well
; HISTORY:  ;   CREATED:  03/28/2012  Lynn B. Wilson III
;    
; NOTES:          
;               1)  The output is the sample rate in [# samples per unit time]
;               2)  If GAP_THRESH is set too small, then the returned result is a 
;
; VERSION: 
;   $LastChangedBy: aaronbreneman $
;   $LastChangedDate: 2014-10-15 12:38:05 -0700 (Wed, 15 Oct 2014) $
;   $LastChangedRevision: 15997 $
;   $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/spdsoft/tags/spedas_6_0/general/missions/rbsp/efw/utils/rbsp_sample_rate.pro $
;-

FUNCTION rbsp_sample_rate,time,GAP_THRESH=gap_thresh,AVERAGE=average,OUT_MED_AVG=medavg

;;----------------------------------------------------------------------------------------
;;  Define dummy variables
;;----------------------------------------------------------------------------------------
f        = !VALUES.F_NAN
d        = !VALUES.D_NAN
medavg   = [d,d]
;;----------------------------------------------------------------------------------------
;;  Check input
;;----------------------------------------------------------------------------------------
IF (N_PARAMS() NE 1) THEN RETURN,0

tt       = REFORM(time)
nt       = N_ELEMENTS(tt)
;;  Sort input just in case
sp       = SORT(tt)
tt       = tt[sp]
;;  Define the total time between the first and last data point
trange   = MAX(tt,/NAN) - MIN(tt,/NAN)
;;  Define shifted difference
lower    = LINDGEN(nt - 1L)
upper    = lower + 1L
sh_diff  = [d,(tt[upper] - tt[lower])]
;;----------------------------------------------------------------------------------------
;;  Determine the maximum allowable gap and remove "bad" values
;;----------------------------------------------------------------------------------------
;IF NOT KEYWORD_SET(gap_thresh) THEN mx_gap = trange ELSE mx_gap = gap_thresh[0]
IF (N_ELEMENTS(gap_thresh) EQ 0) THEN mx_gap = trange[0] ELSE mx_gap = gap_thresh[0]
bad      = WHERE(ABS(sh_diff) GT mx_gap[0] OR ABS(sh_diff) EQ 0,bd,COMPLEMENT=good,NCOMPLEMENT=gd)
IF (bd GT 0 AND bd LT nt) THEN BEGIN
  sh_diff[bad] = d
ENDIF ELSE BEGIN
  IF (bd EQ nt) THEN RETURN,d
ENDELSE
;;----------------------------------------------------------------------------------------
;;  Check for outliers and remove if necessary
;;----------------------------------------------------------------------------------------
avg_dt   = MEAN(sh_diff,/NAN,/DOUBLE)
med_dt   = MEDIAN(sh_diff,/DOUBLE)
top      = ABS(avg_dt[0]) > ABS(med_dt[0])
bot      = ABS(avg_dt[0]) < ABS(med_dt[0])
;;  Make sure there is not a significant difference between Avg. and Med.
ratio    = ABS(top[0]/bot[0])*1d1
test     = (ratio[0] GE 1) AND (ABS(avg_dt[0]) NE ABS(med_dt[0]))  ;;  If true, then start loop to eliminate outliers
;test     = 1
jj       = 0L
IF (test) THEN BEGIN
  std_dt   = STDDEV(sh_diff,/NAN,/DOUBLE)
  tb       = avg_dt[0] + std_dt[0]*3d0*[-1d0,1d0]
  ;;  Check for outliers
  temp     = ((sh_diff GE tb[1]) OR (sh_diff LE tb[0]))
  bad      = WHERE(temp,bd)
  IF (bd GT 0 AND bd LT nt) THEN BEGIN
    sh_diff[bad] = d
    avg_dt       = MEAN(sh_diff,/NAN,/DOUBLE)
    med_dt       = MEDIAN(sh_diff,/DOUBLE)
  ENDIF
ENDIF
;WHILE (test) DO BEGIN
;  ;;  Calculate average and median ∆t
;  min_dt   = MIN(ABS(sh_diff),/NAN)
;  max_dt   = MAX(ABS(sh_diff),/NAN)
;  avg_dt   = MEAN(sh_diff,/NAN,/DOUBLE)
;  std_dt   = STDDEV(sh_diff,/NAN,/DOUBLE)
;  med_dt   = MEDIAN(sh_diff,/DOUBLE)
;  tb       = avg_dt[0] + std_dt[0]*3d0*[-1d0,1d0]
;  ;;  Make sure there is not a significant difference between Avg. and Med.
;  ratio    = ABS(avg_dt[0]/med_dt[0]*1d2)
;;  temp     = ((ABS(sh_diff) GE tb[1]) OR (ABS(sh_diff) LE tb[0]))
;  temp     = ((sh_diff GE tb[1]) OR (sh_diff LE tb[0]))
;  bad      = WHERE(temp,bd)
;  jj      += test
;  test     = (bd GT 0) AND (jj LT 10)  ;;  keep from iterating too long...
;  IF (test) THEN sh_diff[bad] = d
;  PRINT,';;  bd = ', bd
;  ration   = ABS(med_dt[0]/min_dt[0]*1d2)
;  ratiox   = ABS(max_dt[0]/med_dt[0]*1d2)
;  test     = ((ratiox[0] GE 1) OR (ration[0] GE 1)) AND (jj LT 100)  ;;  keep from iterating too long...
;  jj      += test
;  IF (test) THEN BEGIN
;    tb  = avg_dt[0] + std_dt[0]*3d0*[-1d0,1d0]
;    bad = WHERE(ABS(sh_diff) GE tb[1] OR ABS(sh_diff) LE tb[0],bd)
;    IF (bd GT 0) THEN sh_diff[bad] = d
;  ENDIF
;ENDWHILE
;;----------------------------------------------------------------------------------------
;;  Calculate the sample rate
;;----------------------------------------------------------------------------------------
samrates = 1d0/sh_diff
avgsmrt  = MEAN(samrates,/NAN,/DOUBLE)
medsmrt  = MEDIAN(samrates,/DOUBLE)
IF KEYWORD_SET(average) THEN sam_rate = avgsmrt[0] ELSE sam_rate = samrates
medavg   = [medsmrt[0],avgsmrt[0]]
;;----------------------------------------------------------------------------------------
;;  Return to user
;;----------------------------------------------------------------------------------------

RETURN,sam_rate
END