;---------------------------------------------------------------------------
;+
; NAME:
;     avsig_1
; PURPOSE:
;     Average and dispersion of an array, zeros can be not included,
; CALLING SEQUENCE:
;     xbar=avsig_1(x, sigma, no_zeros=no_zeros, sig_mean=sig_mean, dimension=dimension, /fractional)
; INPUT:
;     x = an array
; OUTPUT:
;     xbar = mean, total(x)/n_elements(x)
;     sigma = standard deviation, sqrt(total((x-xbar)^2/(nx-1)))
; KEYWORDS:
;     no_zeros= if set, strip out zeros, note that this option does
;               not work if the dimension keyword is used. At least not as of
;               2-13-95...
;     sig_mean = sigma/sqrt(nx), the standard deviation of the mean of the array
;     dimension = the dimension of the array to find the mean in,
;                 passed into the total command, it must be a scalar.
;     fractional = the fractional error is passed out as sigma,
;                  don't use this if zero is a valid value of xbar...
; HISTORY:
;     12-9-94, jmm
;     2-13-95, jmm added dimension keyword, switched from ok_zeros to no_zeros
;     5-sep-1996, jmm, switched to double precision
;-
FUNCTION Avsig_1, x0, sigma, no_zeros=no_zeros, sig_mean=sig_mean, $
                  dimension=dimension, fractional=fractional

   x = double(x0)
   IF(KEYWORD_SET(dimension)) THEN BEGIN
;get the size of the array, if it's only 1d, goto the else part
      size_x = size(x)
      IF(size_x(0) LE 1) THEN GOTO, its_a_vector
;Ok, now be sure that the given dimension exists
      d0 = long(dimension(0))
      IF(d0 GT size_x(0)) THEN BEGIN
         message, 'NOT ENOUGH DIMENSIONS IN ARRAY, USING WHOLE ARRAY', /info
         GOTO, its_a_vector
      ENDIF
;N-elements
      nx = float(size_x(d0))
      xbar = total(x, d0)/nx ;the mean is easy
      sigma = total(x^2, d0) ;gotta use two sums

      IF(nx GT 1.0) THEN sigma = sqrt((sigma-nx*xbar^2)/(nx-1.0)) $
        ELSE sigma = make_array(size = size(xbar)) + 0.0
      sig_mean = sigma/sqrt(nx)
   ENDIF ELSE BEGIN
      its_a_vector: IF(KEYWORD_SET(no_zeros)) THEN BEGIN
         ok = where(x GT 0.0, nx)
         nx = float(nx)
         IF(nx GT 0) THEN BEGIN
            xbar = total(x(ok))/nx
            IF(nx GT 1) THEN sigma = sqrt(total((x(ok)-xbar)^2)/(nx-1.0)) $
              ELSE sigma = 0.0
            sig_mean = sigma/sqrt(nx)
         ENDIF ELSE BEGIN
            xbar = 0.0
            sigma = 0.0
            sig_mean = 0.0
         ENDELSE
      ENDIF ELSE BEGIN
         nx = float(N_ELEMENTS(x))
         xbar = total(x)/nx
         IF(nx GT 1) THEN sigma = sqrt(total((x-xbar)^2)/(nx-1.0)) $
           ELSE sigma = 0.0
         sig_mean = sigma/sqrt(nx)
      ENDELSE
   ENDELSE

   IF(KEYWORD_SET(fractional)) THEN BEGIN
      xxx = where(xbar NE 0.0)
      fraction = xbar & fraction(*) = 0.0
      fraction(xxx) = sigma(xxx)/xbar(xxx)
      sigma = fraction
   ENDIF

   RETURN, xbar
END

;+
; NAME:
;     hsi_obs_summ_bck
; PURPOSE:
;     Given a stream of data, find the appropriate background
;     levels, and the times, for the background interval and
;     for where the data is a certain specified level above
;     the background. Data is assumed to have the same accumulation
;     time for each interval.
; CALLING SEQUENCE:
;     hsi_obs_summ_bck, data, dt, bck_ss, ratio_f, $
;        nbmax = nbmax, nbmin = nbmin, nbinc=nbinc, $
;        dtb_max = dtb_max, dtb_min=dtb_min,  dtb_inc=dtb_inc, $
;        zero_ok=zero_ok, two=two, use_minimum=use_minimum, $
;        quiet=quiet, sig_ratio=sig_ratio
;        data_is_a_rate=data_is_a_rate, $
;        b_ratio=b_ratio, stops=stops, plot = plot
; INPUT:
;     data = an array of data, a vector
;     dt = the accumulation interval for each data point
; OUTPUT:
;     ratio_f = is the final value of the ratio of sigma(counts)/sqrt(counts),
;               -1 is passed out, if we're not satisfied with the result.
;     bck_ss = start and end subscripts for background interval(s)
; KEYWORDS:
;     nbmax = the interval in number of data points, that you want the
;             background accumulated for, the default is 30.
;     nbmin = minimum interval, in number of data points, for the background
;               accumulation, the default is 5. 
;     nbinc = interval increment, in number of datapoints, for the background
;               accumulation, the default is 5.
;     dtb_max = the interval in seconds, that you want the background
;           accumulated for, the default is to use nbmax=30
;           But the actual interval time may be somewhat less, depending
;           on how good the data is, also the background is accumulated
;           for an integral number of intervals, This number of
;           intervals, when summed, may not equal the input dtb.
;           Only used if set, with no nbmax, the default is to use 30 pts
;     dtb_min = minimum interval, in seconds, for the background
;               accumulation, default is to use nbmin. This is automatically
;               reset to max(interval time) if it is smaller than
;               the largest interval. HXT-LO data almost invariably
;               gets to this point.
;     dtb_inc = interval increment, in seconds, for the background
;               accumulation, the default is to use nbinc
;     zero_ok = if set, this will allow zero values in the data.
;     two = if set, find two intervals, one before the max value
;           and one after. dtb_max and dtb_min may be given two values
;           the background level will be found by interpolation.
;     quiet = if set, don't plot the data and show intervals.
;      plot = if set, then show plots, if not set, and quiet
;             is not set there will be some info on the screen
;     data_is_a_rate = set this if you pass in a count rate, and not counts
;     sig_ratio = the routine will reject intervals if
;                 sigma(counts)/sqrt(counts) is greater than this value,
;                 and will reduce dtb to try again. Default is 1.25
;     b_ratio = the routine will reject intervals if the chosen background
;               level over the minimum is more than this value. Default is
;               1.25
;     use_minimum= short circuits the whole program, just use the interval
;                  that has the minimum counts for the given time interval,
;                  This is good for GOES and HXT.
; RESTRICTIONS:
;     Use some sense, this is designed to find the background levels near
;     a flare, don't use it on a whole days worth of data. Don't use it
;     if there are big gaps in the data, if there are long stretches of
;     zeros, don't set the zero_ok keyword. If it looks like the flare
;     starts before the data, as if you only have flare mode data,
;     don't use the two keyword.
; HISTORY:
;     27-jul-94, JMM
;     12-jul-95, added the /two option, finally, nbcks and nsigmas keywords
;                and trange output are gone, they don't make sense without
;                doing the interpolation...
;     26-jul-95, changed the name from bck_finder to auto_bck_find
;     12-feb-96, jmm, added use_minimum keyword
;     21-jun-1999, jmm, changed time formats, and default ratios
;                changed ()'s to []'s for IDL 5 where applicable
;     23-jun-1999, jmm, slightly different, now the approved background
;                       is the min value with sigma<sig_ratio, and
;                       not the min value with sigma<sig_ratio*min_sigma
;     13-sep-1999, jmm, changed time formats to use ANYTIM format
;     14-Apr-2000, jmm, changed name again, Assumes constant interval
;                       times, added nbmin, nbmax, nbinc keywords,
;                       Now a procedure
;     6-sep-2000, jmm, added plot keyword, the default is now to
;                     to show some diagnostic info, but not to plot
;     30-jan-2002, jmm, resets 0 data to 1.0e20 to avoid data gaps
;                       being used as background intervals.
;-
PRO Hsi_obs_bck0, info, data, dt, bck_rate, bck_ss, ratio_f, ratio_b, $
                  use_minimum=use_minimum, quiet=quiet, stops=stops, $
                  Nbmax=nbmax, Nbmin=nbmin, Nbinc=nbinc, $
                  Max_sig_ratio=max_sig_ratio, Max_b_ratio=max_b_ratio


   IF(KEYWORD_SET(quiet)) THEN qt = 0 ELSE qt = 1
   ndset = N_ELEMENTS(data)
   nds1 = ndset-1
   dt0 = dt[0]
   
   ratio_f = -1
   ratio_b = -1
   bck_rate = 0.0               ;don't return this undefined
   b_lvl_f = 0.0                ;don't return this undefined
   bck_ss = [0, nds1]           ;don't return this undefined
   
;Set trial interval sizes
   IF(datatype(info) EQ 'STC') THEN BEGIN
      nbmx = info.nbmax
      nbmn = info.nbmin
      nbnc = info.nbinc
      ratio0 = info.max_sig_ratio
      ratiob = info.max_b_ratio
      info.time_intv = dt0
   ENDIF ELSE BEGIN
      IF(KEYWORD_SET(nbmax)) THEN nbmx = nbmax ELSE nbmx = 15
      IF(KEYWORD_SET(nbmin)) THEN nbmn = nbmin ELSE nbmn = 5
      IF(KEYWORD_SET(nbinc)) THEN nbnc = nbinc ELSE nbnc = 5
      IF(KEYWORD_SET(max_sig_ratio)) THEN ratio0 = max_sig_ratio $
        ELSE ratio0 = 1.25
      IF(KEYWORD_SET(max_b_ratio)) THEN ratiob = max_b_ratio $
        ELSE ratiob = 1.25
   ENDELSE

;Checking
   nbmn = nbmn > 1
   nbmx = nbmx > nbmn
   nbnc = nbnc > 1

;prepare the data
   y0 = data*dt0                ;data_is_a_rate
;Zero data will bomb this routine, set zero points to 1.0e20
   zss0 = where(y0 EQ 0, nzss0)
   IF(nzss0 GT 0) THEN BEGIN
      IF(nzss0 GT ndset-nbmx) THEN BEGIN
         IF(qt) THEN message, /info, 'Not enough non-zero data points'
         RETURN
      ENDIF
      y0[zss0] = 1.0e20      ;oh my, how long has this been nzss0 here
    ENDIF
;Ok, first find the minimum ratio of sigma^2/counts
;and get the minimum count rate over the dtb sized interval
   nbi = nbmx                   ;the starting interval size
   IF(nbi GT ndset) THEN BEGIN  ;reduce the nbi if necessary
      WHILE(nbi GT ndset) DO nbi = nbi - nbnc
      IF(nbi LT nbmn) THEN BEGIN
         IF(qt) THEN message, /info, 'NOT enough data points'
         RETURN
      ENDIF
   ENDIF
   bloop: IF(qt) THEN message, '***CALCULATING RATIOS AND BACKGROUNDS', /info

   temp_array = fltarr(nbi, ndset-nbi)
   FOR j = 0, nds1-nbi DO temp_array[*, j] = y0[j:j+nbi-1]
   y00b = avsig_1(temporary(temp_array), sy002, dimen=1)

   ratio_all = fltarr(ndset-nbi)+1.0e20
   b_lvl_all = ratio_all
   b0_all = lindgen(ndset-nbi) & b1_all = b0_all+nbi-1 ;oops, 12-dec-2003
   xaa = where((y00b GT 0.0) AND (sy002 GT 0.0))
   IF(xaa[0] NE -1) THEN BEGIN
      ratio_all[xaa] = sy002[xaa]/sqrt(y00b[xaa])
      b_lvl_all[xaa] = y00b[xaa]
   ENDIF
;Ok, that's fine, but the background interval may contain some flare-like
;data, how do i tell? 1st guess, minimum background, ratio must be less than
;ratio0*(minimum ratio), if not then try the 2nd min background, etc.
;until you get an acceptable answer, if no answers come, then you change
;the value of dtb      
   ratio_min = min(ratio_all, min_r)
   ratio_tst = ratio0
   b_lvl_min = min(b_lvl_all, min_b)
   ratio_bmin = ratio_all[min_b]
   b_lvl_tst = ratiob*b_lvl_min
;   plt = 0
;   IF(plt) THEN BEGIN
;      print, ' NDSET =', nds1+1
;      print, 'b0_minb =', b0_all[min_b], '   b1_minb =', b1_all[min_b]
;      print, '*** RATIO_MIN =', ratio_min
;      print, '*** RATIO_TST =', ratio_tst
;      print, '*** B_MIN = ', b_lvl_min
;      print, '*** RATIO AT B_MIN = ', ratio_bmin
;      plot, y0/dt0, psym = 10, title = 'COUNT RATE', xstyle = 1
;      yarr = [!Cymin, !Cymax]
;      oplot, [b0_all[min_b], b0_all[min_b]], yarr
;      oplot, [b1_all[min_b], b1_all[min_b]], yarr
;   ENDIF

   IF(KEYWORD_SET(use_minimum)) THEN BEGIN
      IF(qt) THEN BEGIN
         print, '*** MINIMUM B_LEVEL WILL BE USED'
         print, '*** FINAL NBI =', nbi, ' PTS'
      ENDIF
      b0_f = b0_all[min_b] & b1_f = b1_all[min_b]
      b_lvl_f = b_lvl_min
      ratio_f = ratio_bmin
      ratio_b = 1.0
   ENDIF ELSE BEGIN
      IF(ratio_tst EQ 0) THEN BEGIN ;data gap!!!
         IF(qt) THEN BEGIN
            print, '*** MINIMUM RATIO IS ZERO'
            print, '*** MINIMUM B_LEVEL WILL BE USED'
            print, '*** FINAL NBI =', nbi, ' PTS'
         ENDIF
         b0_f = b0_all[min_b] & b1_f = b1_all[min_b]
         b_lvl_f = b_lvl_min
         ratio_f = ratio_bmin
         ratio_b = 1.0
      ENDIF ELSE BEGIN
         aa = where(ratio_all LT ratio_tst)
         IF(aa[0] EQ -1) THEN BEGIN
;No good intervals, change nbi and start over              
            IF(nbi GT nbmn) THEN BEGIN
               IF(nbi GE nbmn+nbnc) THEN nbi = nbi-nbnc ELSE nbi = nbmn
               IF(qt) THEN BEGIN
                  print, '*** ALL RATIOS ARE GREATER THAN', ratio_tst
                  print, '*** TRYING NBI =', nbi, ' PTS'
               ENDIF
               ratio_f = -1
               GOTO, bloop
            ENDIF ELSE BEGIN
               IF(qt) THEN BEGIN
                  print, '*** ALL RATIOS GREATER THAN', ratio_tst
                  print, '*** MINIMUM B_LEVEL WILL BE USED'
                  print, '*** FINAL NBI =', nbi, ' PTS'
               ENDIF
               b0_f = b0_all[min_b] & b1_f = b1_all[min_b]
               b_lvl_f = b_lvl_min
               ratio_f = ratio_bmin
               ratio_b = 1.0
            ENDELSE
         ENDIF ELSE BEGIN
            b_lvl_f = min(b_lvl_all[aa], min_br)
            ratio_f = ratio_all[aa[min_br]]
            IF(b_lvl_min GT 0.) THEN ratio_b = b_lvl_f/b_lvl_min ELSE BEGIN
               ratio_b = 1.0e20
               IF(qt) THEN BEGIN
                  print, '*** MINIMUM B_LEVEL IS ZERO'
                  print, '*** FINAL NBI =', nbi, ' PTS'
               ENDIF
            ENDELSE
            b0_f = b0_all[aa[min_br]]
            b1_f = b1_all[aa[min_br]]
         ENDELSE
      ENDELSE
   ENDELSE
   IF(qt) THEN BEGIN
      print, '*** FINAL RATIO = ', ratio_f
      print, '*** FINAL B_LEVEL = ', b_lvl_f
      print, '*** FINAL NBI = ', nbi, ' PTS'
   ENDIF
;Ok, now you're out, obtain the intervals
   bck_ss = [b0_f, b1_f]
   bck_rate = b_lvl_f/dt
         
   IF(KEYWORD_SET(stops)) THEN stop
   RETURN
END