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)
  first = array[0]
  last = array[narr-1]
  array_x = [replicate(first, n), temporary(array), replicate(last, 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
;                  regardless of input type. (If input data is double
;                  calculation is always done 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.
;                      
; interp_resolution = If time interpolation is being used, set this
;                     option to control the number of seconds between
;                     interpolated samples.  The default is to use
;                     the value of the smallest separation between 
;                     samples.  Any number higher than this will sacrifice
;                     output resolution to save memory. (NOTE: This option
;                     will not be applied if no interpolation is being
;                     performed because either (1) no_time_interp is set or
;                     (2) the sample rate of the data is constant)
;                     
; interactive_warning = if keyword is set pops up a message box if there are memory problems and asks
;                     the user if they would like to continue
;                     
; interactive_varname = set this to a string indicating the name of the quantity to be used in the warning message.                 
;                     
; warning_result = assign a named variable to this keyword to determine the result of the computation                 
;  
;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
; 23-apr-2008, pcruce, Added padding for no_time_interp option, added _extra keyword
; 28-apr-2008, pcruce, Added interp_resolution option, added memory warning, 
;                        mod to guarantee that precision of output is at least as 
;                        large as precision of input
;$LastChangedBy: jimmpc $
;$LastChangedDate: 2009-05-29 15:09:16 -0700 (Fri, 29 May 2009) $
;$LastChangedRevision: 6003 $
;$URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/ssl_general/tags/tdas_5_1/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,$
                         interp_resolution = interp_resolution,$
                         interactive_warning=interactive_warning,$
                         interactive_varname=interactive_varname,$
                         warning_result=warning_result
                         
  out_array = -1                ;initialize
  warning_result = 1
;; 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)) || is_num(array,/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


    ;Note:  The loop below could probably be speed-optimized by use of the value_locate routine
    ;which would prevent the where function from being called on every iteration
    ;This might entail the need to allocate copies of the inputs for the duration
    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
      
      ;padding is done in-place.  This probably entails a speed hit because the operation is repeated,
      ;But it is assumed that the /no_time_interp is being used because the user values space over time
      ss = where([time_array[0]-dt/2.0,time_array,time_array[n-1]+dt/2.0] Lt t1 And [time_array[0]-dt/2.0,time_array,time_array[n-1]+dt/2.0] Ge t0, nss)
      
      ;; Check if subscripts available
      
      If(nss Gt 0) Then Begin
       
        For k = 0l, nv-1 Do Begin
        
          ok = where(finite(([array[0,k],array,array[n-1,k]])[ss, k]), nok) ;Do not include NaN's
          
          If(nok Gt 0) Then out_array[j, k] = total(([array[0,k],array,array[n-1,k]])[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

          if keyword_set(interp_resolution) then begin
            dtx0 = interp_resolution
          endif
          
          nr = ceil((tx[nok-1]-tx[0])/dtx0)
          t1 = tx[0]+dtx0*dindgen(nr)
        Endif Else Begin
          t1 = tx
          nr = nok
        Endelse
        
        ;first loop warning on memory allocation
        if k eq 0 then begin
          
          if is_num(out_array,/double) then begin
            vector_mem_factor = 2
          endif else begin
            vector_mem_factor = 1
          endelse
          
          ;elements in target * 4 bytes/word * ((2 words per time element)+(1 or 2 words per data element))
          ;divided by 1024 bytes per kB and 1024 kB per mB
          mem_usage_mb_interp = 4D*n_elements(t1)*(2D + vector_mem_factor)/(1024D^2)
          
          if keyword_set(true_t_integration) then begin
            ;(padding elements(front and back) + n of data elements)*4 bytes/word*(1 or 2 bytes per data element)
            ;divided by 1024 bytes per kB and 1024 kB per mB
            mem_usage_mb_smooth = (nrv * 2 + n_elements(t1))*4*vector_mem_factor*2/(1024D^2)
          endif else begin
            mem_usage_mb_smooth = 0
          endelse
          
          mem_usage_total = max([mem_usage_mb_smooth,mem_usage_mb_interp])
          
          ;only warn if memory usage is significant
          if mem_usage_total gt 100 then begin
            ;because of temporary allocation during operations, memory allocation may be as much as double the declared usage
            if keyword_set(interactive_warning) then begin
              if keyword_set(interactive_varname) then begin
                 s = 'WARNING: Operation on ' + interactive_varname + ' will take between ' + strtrim(mem_usage_total,2) + ' MB and ' + strtrim(2*mem_usage_total,2) + ' MB of memory. Do you want to continue?'
              endif else begin
                 s = 'WARNING: Operation will take between ' + strtrim(mem_usage_total,2) + ' MB and ' + strtrim(2*mem_usage_total,2) + ' MB of memory. Do you want to continue?'
              endelse
              ok = dialog_message(s,/question,/center)
              if strlowcase(ok) eq 'no' then begin
                warning_result = 0
                return,out_array     
              endif
            endif else begin
              message,/info,'WARNING: Operation will take between ' + strtrim(mem_usage_total,2) + ' MB and ' + strtrim(2*mem_usage_total,2) + ' MB of memory'
            endelse
          endif
          
        endif
        
        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