; $Id: c_correlate.pro,v 1.20 2004/01/21 15:54:48 scottm Exp $
;
; Copyright (c) 1995-2004, Research Systems, Inc.  All rights reserved.
;       Unauthorized reproduction prohibited.
;+
; NAME:
;       C_CORRELATE
;
; PURPOSE:
;       This function computes the cross correlation Pxy(L) or cross
;       covariance Rxy(L) of two sample populations X and Y as a function
;       of the lag (L).
;
; CATEGORY:
;       Statistics.
;
; CALLING SEQUENCE:
;       Result = C_correlate(X, Y, Lag)
;
; INPUTS:
;       X:    An n-element vector of type integer, float or double.
;
;       Y:    An n-element vector of type integer, float or double.
;
;     LAG:    A scalar or n-element vector, in the interval [-(n-2), (n-2)],
;             of type integer that specifies the absolute distance(s) between
;             indexed elements of X.
;
; KEYWORD PARAMETERS:
;       COVARIANCE:    If set to a non-zero value, the sample cross
;                      covariance is computed.
;
;       DOUBLE:        If set to a non-zero value, computations are done in
;                      double precision arithmetic.
;
; EXAMPLE
;       Define two n-element sample populations.
;         x = [3.73, 3.67, 3.77, 3.83, 4.67, 5.87, 6.70, 6.97, 6.40, 5.57]
;         y = [2.31, 2.76, 3.02, 3.13, 3.72, 3.88, 3.97, 4.39, 4.34, 3.95]
;
;       Compute the cross correlation of X and Y for LAG = -5, 0, 1, 5, 6, 7
;         lag = [-5, 0, 1, 5, 6, 7]
;         result = c_correlate(x, y, lag)
;
;       The result should be:
;         [-0.428246, 0.914755, 0.674547, -0.405140, -0.403100, -0.339685]
;
; PROCEDURE:
;       See computational formula published in IDL manual.
;
; REFERENCE:
;       INTRODUCTION TO STATISTICAL TIME SERIES
;       Wayne A. Fuller
;       ISBN 0-471-28715-6
;
; MODIFICATION HISTORY:
;       Written by:  GGS, RSI, October 1994
;       Modified:    GGS, RSI, August 1995
;                    Corrected a condition which excluded the last term of the
;                    time-series.
;       	     - GGS, RSI, April 1996
;                    Simplified CROSS_COV function. Added DOUBLE keyword.
;                    Modified keyword checking and use of double precision.
;       	     - W. Biagiotti,  Advanced Testing Technologies
;       	     Inc., Hauppauge, NY, July 1997, Moved all
;       	     constant calculations out of main loop for
;       	     greatly reduced processing time.
;   CT, RSI, September 2002. Further speed improvements, per W. Biagiotti.
;                Now handles large vectors and complex inputs.
;-
function C_Correlate, X, Y, Lag, Covariance = Covariance, Double = doubleIn

    compile_opt idl2

    ; Compute the sample cross correlation or cross covariance of
    ; (Xt, Xt+l) and (Yt, Yt+l) as a function of the lag (l).

    ON_ERROR, 2

    typeX = SIZE(X, /TYPE)
    typeY = SIZE(Y, /TYPE)
    nX = N_ELEMENTS(x)

    if (nX ne N_ELEMENTS(y)) then $
        MESSAGE, "X and Y arrays must have the same number of elements."

    ;Check length.
    if (nX lt 2) then $
        MESSAGE, "X and Y arrays must contain 2 or more elements."

    isComplex = (typeX eq 6) or (typeX eq 9) or $
        (typeY eq 6) or (typeY eq 9)

    ;If the DOUBLE keyword is not set then the internal precision and
    ;result are identical to the type of input.
    useDouble = (N_ELEMENTS(doubleIn) eq 1) ? $
        KEYWORD_SET(doubleIn) : $
        (typeX eq 5 or typeY eq 5) or (typeX eq 9 or typeY eq 9)

    ; This will now be in double precision if Double is set.
    Xd = x - TOTAL(X, Double = useDouble) / nX ;Deviations
    Yd = y - TOTAL(Y, Double = useDouble) / nX

    nLag = N_ELEMENTS(Lag)

    Cross = useDouble ? $
        (isComplex ? DCOMPLEXARR(nLag) : DBLARR(nLag)) : $
        (isComplex ? COMPLEXARR(nLag) : FLTARR(nLag))

    M = ABS(Lag)
    for k = 0L, nLag-1 do begin
        ; Note the reversal of the variables for negative lags.
        Cross[k] = (Lag[k] ge 0) ? $
            TOTAL(Xd[0:nX - Lag[k] - 1L] * Yd[Lag[k]:*]) : $
            TOTAL(Yd[0:nX + Lag[k] - 1L] * Xd[-Lag[k]:*])
    endfor

    ; Divide by N for covariance, or divide by variance for correlation.
    Cross = TEMPORARY(Cross) / $
        (KEYWORD_SET(Covariance) ? nX : SQRT(TOTAL(Xd^2)*TOTAL(Yd^2)))

    return, useDouble ? Cross : $
        (isComplex ? COMPLEX(Cross) : FLOAT(Cross))

end