pro lagcovar_window,data1,data2,hz,lcov,tlag,tpeak,maxwin,minwin,badval
;========================================================
; Rewritten by Ken Davis in Feb 95. Modified June 95.
; Modified by K. Davis Nov. 1998 to add a weighting window for selecting
; the peak of the lagged covariance.
; S Mayor, Sept 94 
;========================================================
; VERSION SCREENS FOR BAD DATA VALUES.
; Program uses ffts to compute the lagged covariance
; between two input data arrays.
; Program then selects the peak of the lagged covariance and
; returns the peak value and the time lag which corresponds to
; this peak value.  
; If data2 lags data1, the lag is defined as negative and
;   TPEAK should have a negative value.
; Applies a narrow window for the peak search based on analysis
; of a search under unrestricted conditions.  
; INPUTS
;  DATA1 - first data array.  One dimensional.
;  DATA2 - second data array.  Compute lagged covariance 
;    between DATA1 and DATA2.  Must be the same size.
;  HZ - data frequency (1/time).
;  MINWIN - the minimum allowed lag time.
;  MAXWIN - the maximum allowed lag time.
;  BADVAL - the value given invalid data points in DATA1 and 2.
; OUTPUTS
;  LCOV - lagged covariance between DATA1 and DATA2.  Zero
;    lag is found at nrecs/2, where nrecs is computed below.
;    0 location is the most negative lag, nrecs-1 the most
;    positive lag.  Zero lag value = total(data1*data2)/nrecs.
;  TLAG - time lag (units of 1/HZ) corresponding to the
;    data in LCOV.  plot,tlag,lcov.
;  TPEAK - the time (units of 1/HZ) where the absolute value
;    of LCOV is a maximum.  Use to find lag between two variables.

; Count the number of data points in DATA1 (=DATA2)
nrecs=min([n_elements(data1),n_elements(data2)])
; Assure an even number of data points
nrecs = 2*(nrecs/2)
dat1=data1(0:nrecs-1)
dat2=data2(0:nrecs-1)

; Screen out bad values in DAT1 and DAT2
; At this point, replace bad data with the array mean.  In the future, 
;  interpolated to get the missing data points.
b=where((dat1 eq badval) or (dat2 eq badval))
g=where((dat1 ne badval) and (dat2 ne badval),count)
if (count gt (nrecs/4)) then begin
  m1=total(dat1(g))/count
  m2=total(dat2(g))/count
endif else begin
  print,'no good data - skipping lagged covariance'
  printf,8,'no good data - skipping lagged covariance'

lcov = fltarr(nrecs)
tlag = fltarr(nrecs)
for i=0,nrecs-1 do begin 
lcov(i)= badval
tlag(i)= badval
endfor
tpeak=badval
  return
endelse
if (b(0) ne -1 ) then dat1(b)=m1
if (b(0) ne -1 ) then dat2(b)=m2

fdata1 = (fft(dat1,-1))
fdata2 = (fft(dat2,-1))
; compute the lagged covariance via inverse FFT of the forward FFTs
lcov = float(fft(fdata1*conj(fdata2),1))

; Rearrange so that zero lag is the central data point (nrecs/2) with
; the most negative lag at 0 and the most postive lag at nrecs-1.
; Zero time lag point is equal to total(data1*data2)/nrecs.
lcov=[lcov(nrecs/2:nrecs-1),lcov(0:nrecs/2-1)]
; Compute a time lag array (units of 1/HZ) to accompany LCOV.
tlag=(findgen(nrecs)-nrecs/2)*1./hz

; NO - DON'T use this - the triangle weighting distorts the time for peaks 
; which aren't exactly equal to the lagguess.  Use a narrow
; rectangular window.
;
; Apply a triangular weighting to the lagged covariance array centered
; about a best-guess lag value.  This is intended to improve our
; ability to identify our true lag value each hour.
;pts=hwidth*hz
;slope=findgen(pts)/float(pts)
;triangle=[slope,1,reverse(slope)]
;offset=lagguess*hz
;zeros=fltarr(nrecs/2-pts)
;nzeros=n_elements(zeros)
;weight=[zeros,triangle,zeros(0:nzeros-2)]
;weight=shift(weight,offset)
;wlcov=lcov*weight

; Find the maximum correlation to identify the lag between the two time series
; Limit the search to minwin to maxwin values of the time lag.
rpeak=where((tlag ge minwin) and (tlag le maxwin),num)
if ( num gt 0) then begin
    peak=max(abs(lcov(rpeak)))
    ipeak=where(abs(lcov) eq peak,n)
; Return the time lag between data1 and data2 in units of 1/HZ.
    tpeak=tlag(ipeak)

; Limit the search to +/- nrecs/6 about zero lag.
;peak=max(abs(lcov(nrecs/2-nrecs/6:nrecs/2+nrecs/6)))
;ipeak=where(abs(lcov) eq peak,n)
; Return the time lag between data1 and data2 in units of 1/HZ.
;tpeak=tlag(ipeak)

;  The time lag
;  is the location of the peak in the lagged covariance array.  Note that if 
;  more than one element in lcov = max value, then we assume that the peak in 
;  the lagged covariance is not a valid point.

; We also check to see if the peak is significant by comparing the the std.
;  deviation of the lcov array.

; Make sure the peak is unique, not one of many values of lcov, and
; that the peak is not located at the ends of the window.
iminwin=nrecs/2+minwin*5
imaxwin=nrecs/2+maxwin*5
    if ((ipeak(0) ge 0) and (n eq 1) and (ipeak(0) gt iminwin+2) and (ipeak(0) lt imaxwin-2)) then begin

; Compare magnitude of peak to std dev of the lcov array to check significance.
; EXPERIMENTAL - hey, it seems to work just fine.  Keep it.
        out=moment(lcov)	; compute statistics of lcov array
        std=sqrt(out(1))	; compute standard deviation of lcov array
        if (peak le abs(out(0))+2.5*std) then begin
            tpeak=badval
            ipeak=badval
        endif
        
    endif else begin
        tpeak=badval
    endelse

endif else begin   
tpeak=badval 
endelse

return
end