pro gauss_avg, data_in, time_in, data_out, time_out, $
                       dt=dt, sigma=sigma, verbose=verbose

if (n_params() eq 0) then begin
    Message,'NAME: GAUSS_AVG.PRO',/info
    Message,'',/info
    Message,'PURPOSE: Uses weighted time-domain averaging to convert',/info
    Message,'         irregular-cadence data to regular-cadence data. ',/info
    Message,'',/info
    Message,'USAGE: gauss_avg, data_in, time_in, data_out, time_out, $',/info
    Message,'                  dt=dt, sigma=sigma, verbose=verbose',/info
    Message,'',/info
    Message,'ARGUMENTS: ',/info
    Message,'  DATA_IN: fltarr(nx,ny,nt0) of data irregularly sampled in',/info
    Message,'        time, to be averaged into nt1 regular data slices',/info
    Message,'  TIME_IN: fltarr(nt0) of irregular data times [seconds]',/info
    Message,'  DATA_OUT: fltarr(nx,ny,nt1) regularly sampled data slices',/info
    Message,'  TIME_OUT: fltarr(nt1) of regular data sample times',/info
    Message,'',/info
    Message,'KEYWORDS: ',/info
    Message,'  DT: Regular cadence of DATA_OUT; either user set, or the',/info
    Message,'      computed median of time intervals in TIME_IN',/info
    Message,'  SIGMA: Width of Gaussian averaging window [seconds]',/info
    Message,'  VERBOSE: Set to print progress, other info.',/info
    Message,'',/info
    Message,'BLAME: Written 30 Apr 2007, B.T. Welsch, UCB/SSL',/info
    return
endif

nx = n_elements(data_in(*,0,0))
ny = n_elements(data_in(0,*,0))
nt0 = n_elements(data_in(0,0,*))

dt_in = time_in(1:nt0-1) - time_in(0:nt0-2) ; arr of input time differences 
if not(keyword_set(dt)) then dt = median(dt_in) ; regular (output) time diff.

if not(keyword_set(sigma)) then sigma = 10*dt ; width of averaging window

mint = time_in(0)
maxt = time_in(nt0-1)

nt1 = floor((maxt - mint)/dt) + 1. ; # regular time steps (extrapolates endpt)
time_out = dt*findgen(nt1)      ; array of regular times

if (keyword_set(verbose)) then begin
    print,'Sigma, dt, and # of steps:'
    print,sigma,dt,nt1
endif

; Scan data set to find max. # of time steps within 2*sigma
;============================================================
n_twght = 0.                    ; will be max # time steps w/in window
for i=0,nt1-1 do begin
    tnear = where((time_in ge time_out(i)-sigma) and $
                  (time_in le time_out(i)+sigma), n_tnear) 
    if (n_tnear gt n_twght) then n_twght = n_tnear
endfor

; Define (x,y,n_twght) array of weights
;=======================================================
twght = fltarr(nx,ny,n_twght)   
tgauss = fltarr(n_twght)        ; 1-D array of temporal weights

; Do 0th time step first, then append succeeding steps.
;=======================================================

; Find of data slices "near" time_in(0)
tnear = where(time_in le time_out(0) + sigma, n_tnear) 
near_arr = time_in(tnear)       ; times of nearby data

; Fill in 1-D array of temporal weights
tgauss(0:n_tnear-1) = exp(-(time_out(0) - near_arr)^2/(2*sigma^2))
if (n_tnear lt n_twght) then tgauss(n_tnear:n_twght-1) = 0.

; Normalize the 1-D array of temporal weights
tgauss(0:n_tnear-1) = tgauss(0:n_tnear-1)/total(tgauss(0:n_tnear-1))

; Replicate the 1-D array in time over two spatial dimensions (nx,ny)
for j=0,n_twght-1 do replicate_inplace, $
  twght, tgauss(j), 1, [0,0,j], 2, indgen(ny)

data_out =  total(data_in(*,*,0:n_tnear-1)*twght(*,*,0:n_tnear-1),3)

; If run in verbose mode, display info
;===================================================
if (keyword_set(verbose)) then begin
    print,'Step, # avg. pts., start pt., end pt.'
    print,0,n_tnear,tnear(0),tnear(n_tnear-1)
endif

; nt1 = number of regularly-spaced temporal points
for i=1,nt1-1 do begin

    ; find m'grams "near" in t
    tnear = where((time_in ge time_out(i)-sigma) and $
                  (time_in le time_out(i)+sigma), n_tnear) 

    if (n_tnear gt n_twght) then print,fail_variable ; n_twght is wrong!

    time_arr = time_in(tnear)   ; find times of nearby m'grams

    ; fill in 1-D array of temporal weights
    tgauss(0:n_tnear-1) = exp(-(time_out(i) - time_arr)^2/(2*sigma^2))
    if (n_tnear lt n_twght) then tgauss(n_tnear:n_twght-1) = 0.

    ;normalize the 1-D array of temporal weights
    tgauss(0:n_tnear-1) = tgauss(0:n_tnear-1)/total(tgauss(0:n_tnear-1))

    ; replicate the 1-D array in time over two spatial dimensions (nx,ny)
    for j=0,n_twght-1 do replicate_inplace, $
      twght, tgauss(j), 1, [0,0,j], 2, indgen(ny)

    ; for i near nt0, choose correct end point of temporal window 
    imax = min([tnear(0)+n_tnear-1,nt0-1]) 
    di = imax - tnear(0)

    data_out_i =  total(data_in(*,*,tnear(0):imax)*twght(*,*,0:di),3)

    data_out = [[[data_out]],[[data_out_i]]]

    ; print status if in verbose mode
    if (keyword_set(verbose)) then print,i,n_tnear,tnear(0),tnear(n_tnear-1)


endfor



end