Index: trunk/ToBeReviewed/STATISTICS/a_correlate2d.pro
===================================================================
--- trunk/ToBeReviewed/STATISTICS/a_correlate2d.pro	(revision 21)
+++ trunk/ToBeReviewed/STATISTICS/a_correlate2d.pro	(revision 21)
@@ -0,0 +1,113 @@
+;+
+; 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
Index: trunk/ToBeReviewed/STATISTICS/a_timecorrelate.pro
===================================================================
--- trunk/ToBeReviewed/STATISTICS/a_timecorrelate.pro	(revision 21)
+++ trunk/ToBeReviewed/STATISTICS/a_timecorrelate.pro	(revision 21)
@@ -0,0 +1,167 @@
+;+
+; NAME:
+;       A_TIMECORRELATE
+;
+; PURPOSE:
+;       Same function as A_CORRELATE but accept array (until 4
+;       dimension) for input and do the autocorrelation or the
+;       autocovariance along the time dimension which must be the last
+;       one of the input array.
+;
+;       This function computes the autocorrelation Px(L) or autocovariance
+;       Rx(L) of a sample population X as a function of the lag (L).
+;
+; CATEGORY:
+;       Statistics.
+;
+; CALLING SEQUENCE:
+;       Result = a_timecorrelate(X, Lag)
+;
+; INPUTS:
+;       X:    an Array which last dimension is the time dimension os
+;       size n.
+;
+;     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 autocovariance
+;                      is computed.
+;
+;       DOUBLE:        If set to a non-zero value, computations are done in
+;                      double precision arithmetic.
+;
+; EXAMPLE
+;       Define an n-element sample population.
+;         x = [3.73, 3.67, 3.77, 3.83, 4.67, 5.87, 6.70, 6.97, 6.40, 5.57]
+;
+;       Compute the autocorrelation of X for LAG = -3, 0, 1, 3, 4, 8
+;         lag = [-3, 0, 1, 3, 4, 8]
+;         result = a_correlate(x, lag)
+;
+;       The result should be:
+;         [0.0146185, 1.00000, 0.810879, 0.0146185, -0.325279, -0.151684]
+;
+; PROCEDURE:
+;
+;
+;                         n-L-1 
+;                         sigma  (X[k]-Xmean)(X[k+L]-Xmean)
+;                          k=0
+;   correlation(X,L)=----------------------------------------
+;                            n-1                   
+;                           sigma  (X[k]-Xmean)^2
+;                            k=0                   
+;
+;
+;
+;                        n-L-1 
+;                        sigma  (X[k]-Xmean)(X[k+L]-Xmean)
+;                         k=0
+;   covariance(X,L)=-------------------------------------------
+;                                     n
+;
+;   Where Xmean is the Time mean of the sample population
+;   x=(x[t=0],x[t=1],...,x[t=n-1])
+;
+;
+; REFERENCE:
+;       INTRODUCTION TO STATISTICAL TIME SERIES
+;       Wayne A. Fuller
+;       ISBN 0-471-28715-6
+;
+; MODIFICATION HISTORY:
+;       24/2/2000 Sebastien Masson (smasson@lodyc.jussieu.fr)
+;       Based on the A_CORRELATE procedure of IDL
+;-
+
+FUNCTION TimeAuto_Cov, X, M, nT, Double = Double, zero2nan = zero2nan
+;Sample autocovariance function
+   TimeDim = size(X, /n_dimensions)
+   Xmean = TOTAL(X, TimeDim, Double = Double) / nT
+   if double then one = 1.0d ELSE one = 1.0
+   Xmean = Xmean[*]#replicate(one, nT - M)
+;
+;
+   case TimeDim of
+      1:res = TOTAL((X[0:nT - M - 1L] - Xmean) * (X[M:nT - 1L] - Xmean), $
+                    TimeDim, Double = Double)
+      2:res = TOTAL((X[*, 0:nT - M - 1L] - Xmean) * $
+                    (X[*, M:nT - 1L] - Xmean) $
+                    , TimeDim, Double = Double)
+      3:res = TOTAL((X[*, *, 0:nT - M - 1L] - Xmean) * $
+                    (X[*, *, M:nT - 1L] - Xmean) $
+                    , TimeDim, Double = Double)
+      4:res = TOTAL((X[*, *, *, 0:nT - M - 1L] - Xmean) * $
+                    (X[*, *, *, M:nT - 1L] - Xmean) $
+                    , TimeDim, Double = Double)
+   ENDCASE
+   if keyword_set(zero2nan) then begin
+      zero = where(res EQ 0)
+      if zero[0] NE -1 then res[zero] = !values.f_nan
+   endif
+   RETURN, res
+
+END
+
+FUNCTION A_TimeCorrelate, 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)
+   nT = XDim[XNDim-1]
+                                ;Check length.
+   if nT lt 2 then $
+    MESSAGE, "Time axis of X array must contain 2 or more 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 n_elements(lag) EQ 0 then lag = 0
+   nLag = N_ELEMENTS(Lag)
+
+   if nLag eq 1 then Lag = [Lag] ;Create a 1-element vector.
+
+   case XNDim of
+      1:if Double eq 0 then Auto = FLTARR(nLag) else Auto = DBLARR(nLag)
+      2:if Double eq 0 then Auto = FLTARR(XDim[0], nLag) else Auto = DBLARR(XDim[0], nLag)
+      3:if Double eq 0 then Auto = FLTARR(XDim[0], XDim[1], nLag) $
+      else Auto = DBLARR(XDim[0], XDim[1], nLag)
+      4:if Double eq 0 then Auto = FLTARR(XDim[0], XDim[1], XDim[2], nLag) $
+      else Auto = DBLARR(XDim[0], XDim[1], XDim[2], nLag)
+   endcase
+
+   if KEYWORD_SET(Covariance) eq 0 then begin ;Compute Autocorrelation.
+      for k = 0, nLag-1 do $
+       case XNDim of
+         1:Auto[k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / $
+          TimeAuto_Cov(X, 0L, nT, Double = Double, /zero2nan)
+         2:Auto[*, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / $
+          TimeAuto_Cov(X, 0L, nT, Double = Double, /zero2nan)
+         3:Auto[*, *, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / $
+          TimeAuto_Cov(X, 0L, nT, Double = Double, /zero2nan)
+         4:Auto[*, *, *, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / $
+          TimeAuto_Cov(X, 0L, nT, Double = Double, /zero2nan)
+      endcase
+   endif else begin             ;Compute Autocovariance.
+      for k = 0, nLag-1 do $ 
+       case XNDim of
+         1:Auto[k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / nT
+         2:Auto[*, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / nT
+         3:Auto[*, *, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / nT
+         4:Auto[*, *, *, k] = TimeAuto_Cov(X, ABS(Lag[k]), nT, Double = Double) / nT
+      endcase
+   endelse
+
+   if Double eq 0 then RETURN, FLOAT(Auto) else $
+    RETURN, Auto
+
+END
Index: trunk/ToBeReviewed/STATISTICS/c_timecorrelate.pro
===================================================================
--- trunk/ToBeReviewed/STATISTICS/c_timecorrelate.pro	(revision 21)
+++ trunk/ToBeReviewed/STATISTICS/c_timecorrelate.pro	(revision 21)
@@ -0,0 +1,220 @@
+;+
+; NAME:
+;       C_TIMECORRELATE
+;
+; PURPOSE:
+;       This function computes the "time cross correlation" Pxy(L) or
+;       the "time cross covariance" between 2 arrays (this is some
+;       kind of c_correlate but for multidimenstionals arrays) as a
+;       function of the lag (L).
+;
+; CATEGORY:
+;       Statistics.
+;
+; CALLING SEQUENCE:
+;       Result = c_timecorrelate(X, Y, Lag)
+;
+; INPUTS:
+;       X:   an Array which last dimension is the time dimension of
+;       size n, float or double.
+;
+;       Y:    an Array which last dimension is the time dimension of
+;       size n, 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_timecorrelate(x, y, lag)
+;
+;       The result should be:
+;         [-0.428246, 0.914755, 0.674547, -0.405140, -0.403100, -0.339685]
+;
+; PROCEDURE:
+;
+;
+;   FOR L>=0 
+;
+;                              n-L-1 
+;                              sigma  (X[k]-Xmean)(Y[k+L]-Ymean)
+;                               k=0
+;   correlation(X,Y,L)=------------------------------------------------------
+;                              n-1                     n-1
+;                      sqrt( (sigma  (X[k]-Xmean)^2)*(sigma  (Y[k]-Ymean)^2))
+;                              k=0                     k=0
+;
+;
+;
+;                              n-L-1 
+;                              sigma  (X[k]-Xmean)(Y[k+L]-Ymean)
+;                               k=0
+;   covariance(X,Y,L)=------------------------------------------------------
+;                                            n
+;
+;   FOR L<0
+;
+;
+;                              n-L-1 
+;                              sigma  (X[k+L]-Xmean)(Y[k]-Ymean)
+;                               k=0
+;   correlation(X,Y,L)=------------------------------------------------------
+;                              n-1                     n-1
+;                      sqrt( (sigma  (X[k]-Xmean)^2)*(sigma  (Y[k]-Ymean)^2))
+;                              k=0                     k=0
+;
+;
+;
+;                              n-L-1 
+;                              sigma  (X[k+L]-Xmean)(Y[k]-Ymean)
+;                               k=0
+;   covariance(X,Y,L)=------------------------------------------------------
+;                                            n
+;
+;   Where Xmean and Ymean are the time means of the sample populations
+;   x=(x[t=0],x[t=1],...,x[t=n-1]) and y=(y[t=0],y[t=1],...,y[t=n-1]),
+;   respectively.
+;
+;
+;
+; REFERENCE:
+;       INTRODUCTION TO STATISTICAL TIME SERIES
+;       Wayne A. Fuller
+;       ISBN 0-471-28715-6
+;
+; MODIFICATION HISTORY:
+;       - 01/03/2000 Sebastien Masson (smasson@lodyc.jussieu.fr)
+;       Based on the C_CORRELATE procedure of IDL
+;       - August 2003 Sebastien Masson 
+;       update according to the update made in C_CORRELATE by
+;       W. Biagiotti and available in IDL 5.5
+;-
+
+FUNCTION TimeCross_Cov, Xd, Yd, M, nT, Ndim, Double = Double, ZERO2NAN = zero2nan
+  ;Sample cross covariance function.
+
+  compile_opt hidden
+;
+   case Ndim OF
+      1:res = TOTAL(Xd[0:nT - M - 1L] * Yd[M:nT - 1L] $
+                    , Double = Double)
+      2:res = TOTAL(Xd[*, 0:nT - M - 1L] * Yd[*, M:nT - 1L] $
+                    , Ndim, Double = Double)
+      3:res = TOTAL(Xd[*, *, 0:nT - M - 1L] * Yd[*, *, M:nT - 1L] $
+                    , Ndim, Double = Double)
+      4:res = TOTAL(Xd[*, *, *, 0:nT - M - 1L] * Yd[*, *, *, M:nT - 1L] $
+                    , Ndim, Double = Double)
+   ENDCASE
+   if keyword_set(zero2nan) then begin
+      zero = where(res EQ 0)
+      if zero[0] NE -1 then res[zero] = !values.f_nan
+   ENDIF
+;
+   RETURN, res
+
+END
+
+FUNCTION C_Timecorrelate, X, Y, Lag, Covariance = Covariance, Double = Double
+
+;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
+
+   xsize = SIZE(X)
+   ysize = SIZE(Y)
+   nt = float(xsize[xsize[0]])
+   NDim = xsize[0]
+
+   if total(xsize[0:xsize[0]] NE ysize[0:ysize[0]]) NE 0 then $
+    MESSAGE, "X and Y arrays must have the same size and the same dimensions"
+
+;Check length.
+   if nt lt 2 then $
+    MESSAGE, "Time dimension of X and Y arrays must contain 2 or more 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 = (Xsize[Xsize[0]+1] eq 5 or ysize[ysize[0]+1] eq 5)
+
+   if n_elements(lag) EQ 0 then lag = 0
+   nLag = N_ELEMENTS(Lag)
+
+;Deviations
+   if double then one = 1.0d ELSE one = 1.0
+   Ndim = size(X, /n_dimensions)
+   Xd = TOTAL(X, Ndim, Double = Double) / nT
+   Xd = X - Xd[*]#replicate(one,  nT)
+   Yd = TOTAL(Y, Ndim, Double = Double) / nT
+   Yd = Y - Yd[*]#replicate(one,  nT)
+
+   if nLag eq 1 then Lag = [Lag] ;Create a 1-element vector.
+
+   case NDim of
+      1:if Double eq 0 then  Cross = FLTARR(nLag) else  Cross = DBLARR(nLag)
+      2:if Double eq 0 then  Cross = FLTARR(Xsize[1], nLag) else  Cross = DBLARR(Xsize[1], nLag)
+      3:if Double eq 0 then  Cross = FLTARR(Xsize[1], Xsize[2], nLag) $
+      else  Cross = DBLARR(Xsize[1], Xsize[2], nLag)
+      4:if Double eq 0 then  Cross = FLTARR(Xsize[1], Xsize[2], Xsize[3], nLag) $
+      else  Cross = DBLARR(Xsize[1], Xsize[2], Xsize[3], nLag)
+   endcase
+
+   if KEYWORD_SET(Covariance) eq 0 then begin ;Compute Cross  Crossation.
+      for k = 0, nLag-1 do begin
+         if Lag[k] ge 0 then BEGIN 
+            case NDim of
+               1: Cross[k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double)
+               2: Cross[*, k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double)
+               3: Cross[*, *, k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double)
+               4: Cross[*, *, *, k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double) 
+             endcase
+         ENDIF else BEGIN 
+            case NDim of
+               1: Cross[k] = TimeCross_Cov(Yd, Xd, ABS(Lag[k]), nT, Ndim, Double = Double)
+               2: Cross[*, k] = TimeCross_Cov(Yd, Xd, ABS(Lag[k]), nT, Ndim, Double = Double)
+               3: Cross[*, *, k] = TimeCross_Cov(Yd, Xd, ABS(Lag[k]), nT, Ndim, Double = Double)
+               4: Cross[*, *, *, k] = TimeCross_Cov(Yd, Xd, ABS(Lag[k]), nT, Ndim, Double = Double)
+             endcase
+         ENDELSE 
+       ENDFOR
+       div = sqrt(TimeCross_Cov(Xd, Xd, 0L, nT, Ndim, Double = Double, /zero2nan) * $
+                  TimeCross_Cov(Yd, Yd, 0L, nT, Ndim, Double = Double, /zero2nan))
+       Cross = temporary(Cross)/((temporary(div))[*]#replicate(one, nLag))
+   endif else begin             ;Compute Cross Covariance.
+      for k = 0, nLag-1 do begin
+         if Lag[k] ge 0 then BEGIN
+            case NDim of
+               1: Cross[k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double) / nT
+               2: Cross[*, k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double) / nT
+               3: Cross[*, *, k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double) / nT
+               4: Cross[*, *, *, k] = TimeCross_Cov(Xd, Yd, Lag[k], nT, Ndim, Double = Double) / nT
+            ENDCASE
+         ENDIF else BEGIN 
+            case NDim of
+               1: Cross[k] = TimeCross_Cov(yd, xd, ABS(Lag[k]), nT, Ndim, Double = Double) / nT
+               2: Cross[*, k] = TimeCross_Cov(yd, xd, ABS(Lag[k]), nT, Ndim, Double = Double) / nT
+               3: Cross[*, *, k] = TimeCross_Cov(yd, xd, ABS(Lag[k]), nT, Ndim, Double = Double) / nT
+               4: Cross[*, *, *, k] = TimeCross_Cov(yd, xd, ABS(Lag[k]), nT, Ndim, Double = Double) / nT
+            ENDCASE
+         ENDELSE 
+      endfor
+   endelse
+
+   if Double eq 0 then RETURN, FLOAT(Cross) else RETURN,  Cross
+
+END
+  
