;+
; NAME:
;       A_CORRELATE2d
;
; PURPOSE:
;
;       This function computes the autocorrelation Px(K,L) or
;       autocovariance Rx(K,L) of a sample population X[nx,ny] as a
;       function of the lag (K,L).
;
; CATEGORY:
;       Statistics.
;
; CALLING SEQUENCE:
;       Result = a_correlate2d(X, Lag)
;
; INPUTS:
;       X:    an 2 dimension Array [nx,ny]
;
;     LAG:    2-element vector, in the intervals [-(nx-2), (nx-2)],[-(ny-2), (ny-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 autocovariance
;                      is computed.
;
;       DOUBLE:        If set to a non-zero value, computations are done in
;                      double precision arithmetic.
;
; EXAMPLE:
;
; PROCEDURE:
;
;
;                          nx-k-1  ny-l-1 
;                          sigma   sigma  (X[i,j]-Xmean)(X[i+k,j+l]-Ymean)
;                           i=0     j=0
;   correlation(X,[k,l])=------------------------------------------------------
;                                  nx-1   ny-1                  
;                                 sigma  sigma  (X[i,j]-Xmean)^2)
;                                  i=0    j=0                      
;
;
;                          nx-k-1  ny-l-1 
;                          sigma   sigma  (X[i,j]-Xmean)(Y[i+k,j+l]-Ymean)
;                           i=0     j=0
;   covariance(X,[k,l])=------------------------------------------------------
;                                            nx*ny
;
;   Where Xmean is the mens of the sample population
;   x=(x[0,0],x[1,0],...,x[nx-1,ny-1]).
;
;
; REFERENCE:
;
; MODIFICATION HISTORY:
;       28/2/2000 Sebastien Masson (smasson@lodyc.jussieu.fr)
;       Based on the A_CORRELATE procedure of IDL
;-

FUNCTION Auto_Cov2d, X, Lag, Double = Double, zero2nan = zero2nan
   XDim = SIZE(X, /dimensions)
   nx = XDim[0]
   ny = XDim[1]
;Sample autocovariance function
   Xmean = TOTAL(X, Double = Double) / (1.*nx*ny)
;
   res = TOTAL( (X[0:nx-1-lag[0], 0:ny-1-lag[1]] - Xmean) * $
                (X[lag[0]:nx-1, lag[1]:ny-1] - Xmean) $
                , Double = Double )
   if keyword_set(zero2nan) AND res EQ 0 then res = !values.f_nan
   RETURN, res

END

FUNCTION A_Correlate2d, X, Lag, Covariance = Covariance, Double = Double

;Compute the sample-autocorrelation or autocovariance of (Xt, Xt+l)
;as a function of the lag (l).

   ON_ERROR, 2

   XDim = SIZE(X, /dimensions)
   XNDim = SIZE(X, /n_dimensions)
   nx = XDim[0]
   ny = XDim[1]
   if XNDim NE 2 then $
    MESSAGE, "X array must contain 2 dimensions."
;Check length.
   if nx lt 2 then $
    MESSAGE, "first dimension of X array must contain 2 or more elements."
   if ny lt 2 then $
    MESSAGE, "second dimension of X array must contain 2 or more elements."
   if n_elements(Lag) NE 2 THEN $
    MESSAGE, "Lag array must contain 2 elements."
   
;If the DOUBLE keyword is not set then the internal precision and
;result are identical to the type of input.
   if N_ELEMENTS(Double) eq 0 then $
    Double = (SIZE(X, /type) eq 5)

   if KEYWORD_SET(Covariance) eq 0 then begin ;Compute Autocorrelation.
      Auto = Auto_Cov2d(X, ABS(Lag), Double = Double) / $
          Auto_Cov2d(X, [0L, 0L], Double = Double, /zero2nan)
   endif else begin             ;Compute Autocovariance.
      Auto = Auto_Cov2d(X, ABS(Lag), Double = Double) / n_elements(X) 
   endelse

   if Double eq 0 then RETURN, FLOAT(Auto) else $
    RETURN, Auto

END