Function temp_t_integration, array, n
;simulate a time integration using the smooth function
;at each point, 
;result = n*smooth(array, n)/(n-1)-shift(array,
;-n/2)/(2.0*(n-1))-shift(array, n/2)/(2.0*(n-1))
;put n values on either side of the array to avoid edge issues
  narr = n_elements(array)
  array_x = [replicate(array[0], n), array, replicate(array[narr-1], n)]
  array_x = n*smooth(array_x, n)/(n-1)-$
    shift(array_x,-n/2)/(2.0*(n-1))-$
    shift(array_x, n/2)/(2.0*(n-1))
  return, array_x[n:n+narr-1]

End
;+
;NAME:
; smooth_in_time
;PURPOSE:
; Runs ts_smooth for irregular grids, after regularising grid
;CALLING SEQUENCE:
; ts = smooth_in_time(array, time_array, dt, /backward, /forward,
;                  /double, /no_time_interp)
;INPUT:
; array = a data array, can be 2-d (ntimes, n_something_else), the
;         first index is smoothed or averaged.
; time_array = a time array (in units of seconds)
; dt = the averaging time (in seconds)
; backward = if set, perform an average over the previous dt, the
;            default is to average from t-dt/2 to t_dt/2
; forward = if set, perform an average over the next dt
; double = if set, do calculation in double precision
; no_time_interp = if set, do *not* interpolate the data to the
;                  minimum time resolution. The default procedure is
;                  to interpolate the data to a regularly spaced grid,
;                  and then use ts_smooth to get the running
;                  average. This alternative can be slow.
; smooth_nans = if set, replace Nan values in the input array with the
;               average values calculated using the ts_smooth
;               process. This has not been implemented for the
;               no_time_interp option.
; true_t_integration = if set, subtract 1/2 of the end points of the
;                      integration from each value, to obtain the
;                      value for an integration over time of the
;                      appropriate interval. This has not been
;                      implemented for the no_time_interp option.
;                      Ths is created for the high_pass_filter.
;OUTPUT:
; ts = the data array smoothed or averaged
;HISTORY:
; 13-mar-2008, jmm, jimm@ssl.berkeley.edu, hacked from
; high_pass_filter.pro and added ts_smooth as the default
; 13-mar-2008, ts_smooth is way too slow, just uses smooth.pro now
; 6-may-2008, jmm, added sort for input data for cases with
;                  non-monotonic time_arrays
;$LastChangedBy: jimm $
;$LastChangedDate: 2008-05-13 18:09:41 -0700 (Tue, 13 May 2008) $
;$LastChangedRevision: 3081 $
;$URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/ssl_general/tags/tdas_5_00/misc/smooth_in_time.pro $
;-

Function smooth_in_time, array, time_array, dt, $
                         backward = backward, forward = forward, $
                         double = double, no_time_interp = no_time_interp, $
                         smooth_nans = smooth_nans, $
                         true_t_integration = true_t_integration, $
                         _extra = _extra
  out_array = -1                ;initialize
;; determine number of rows in input array
;; Note: this is a tplot array, reversed from
;; idl array
  n =  n_elements(array[*, 0])
;; Make sure time values exist for each entry
  If(n_elements(time_array) Ne n) Then Begin
    message, /info, 'Array mismatch'
    return,  out_array
  Endif
;; Produces array of values, the first being the dimension of the array
;; which will later be used as a check
  sz = size(array)
  If(sz[0] Eq 2) Then nv = sz[2] Else nv = 1 ;the 2nd index will be looped over
;; Now declare output array, fill with NaN's
  If(keyword_set(double)) Then Begin
    out_array = double(array) & out_array[*] = !values.d_nan
  Endif Else Begin
    out_array = float(array) & out_array[*] = !values.f_nan
  Endelse
;; Do the calculation
  If(keyword_set(no_time_interp)) Then Begin
;; This for loop will take us through the full array of values; this
;; can be very slow
    For j = 0l, n-1 Do Begin 
;; Get subscripts of group to take running average over
;; nss is the number values returned
      If(keyword_set(backward)) Then Begin
        t0 = time_array[j]-dt
        t1 = time_array[j]
      Endif Else If(keyword_set(forward)) Then Begin
        t0 = time_array[j]
        t1 = time_array[j]+dt
      Endif Else Begin
        t0 = time_array[j]-dt/2.0
        t1 = time_array[j]+dt/2.0
      Endelse
      ss = where(time_array Lt t1 And time_array Ge t0, nss)
;; Check if subscripts available
      If(nss Gt 0) Then Begin
        For k = 0l, nv-1 Do Begin
          ok = where(finite(array[ss, k]), nok) ;Do not include NaN's
          If(nok Gt 0) Then out_array[j, k] = total(array[ss[ok]])/nok
        Endfor
      Endif
    Endfor
  Endif Else Begin              ;default behavior is to interpolate
    For k = 0, nv-1 Do Begin
      ok = where(finite(array[*, k]), nok) ;Do not include NaN's
      If(nok Gt 0) Then Begin
        tx = time_array[ok]     ;ok times
        ax = array[ok, k]       ;ok data points
        dtx = tx[1:*]-tx
        bad_dtx = where(dtx Le 0.0, nbad_dtx)
        If(nbad_dtx Gt 0) Then Begin ;sort the data
          message, /info, 'Data is non-monotonic, Sorting...'
          ss_tx = sort(tx)
          tx = tx[ss_tx]
          ax = ax[ss_tx]
          dtx = tx[1:*]-tx
        Endif
        dtx0 = min(dtx[where(dtx Gt 0.0)])         ;min value of t resolution
        not_min = where(abs(dtx-dtx0) Gt dtx0/100.0, cnot_min) ;small allowance
        nrv = ceil(dt/dtx0)
;Note that for non-forward or backwards, this value must be and odd
;number gt 3
        If(nrv Lt 3) Then message, /info, $
          'Number of smoothing points is LT 3, Smoothing over 3*minimum resolution'
        nrv = nrv > 3
        If(nrv Mod 2 Eq 0) Then Begin 
          message, /info, 'Even number of smoothing points:'+$
            strcompress(string(nrv))+', Adding 1'
          nrv = nrv+1
        Endif
;Now do the smoothing      
        If(cnot_min Ne 0) Then Begin
;Create the regular grid
          nr = ceil((tx[nok-1]-tx[0])/dtx0)
          t1 = tx[0]+dtx0*dindgen(nr)
        Endif Else Begin
          t1 = tx
          nr = nok
        Endelse
        out1 = interpol(temporary(ax), temporary(tx), t1) ;interp to hi-res
;get the average, pad backward and forwards if needed
        If(keyword_set(backward)) Then Begin
          out1 = [fltarr(nrv/2)+out1[0], out1]
          If(keyword_set(true_t_integration)) Then Begin
            rout1 = temp_t_integration(out1, nrv)
          Endif Else rout1 = smooth(out1, nrv, /edge_truncate)
          rout1 = rout1[0:nr-1]
        Endif Else If(keyword_set(forward)) Then Begin
          out1 = [out1, fltarr(nrv/2)+out1[nr-1]]
          If(keyword_set(true_t_integration)) Then Begin
            rout1 = temp_t_integration(out1, nrv)
          Endif Else rout1 = smooth(out1, nrv, /edge_truncate)
          rout1 = rout1[nrv/2:*]
        Endif Else Begin
          If(keyword_set(true_t_integration)) Then Begin
            rout1 = temp_t_integration(out1, nrv)
          Endif Else rout1 = smooth(out1, nrv, /edge_truncate)
        Endelse
;And interpolate back to the full time_array
        If(keyword_set(smooth_nans)) Then Begin
          out_array[*, k] = interpol(rout1, t1, time_array)
        Endif Else out_array[ok, k] = interpol(rout1, t1, time_array[ok])
      Endif
    Endfor
  Endelse
  Return, out_array
End