Index: trunk/ToBeReviewed/STATISTICS/c_timecorrelate.pro
===================================================================
--- trunk/STATISTICS/c_timecorrelate.pro	(revision 2)
+++ trunk/ToBeReviewed/STATISTICS/c_timecorrelate.pro	(revision 21)
@@ -97,27 +97,25 @@
 ;
 ; MODIFICATION HISTORY:
-;       01/03/2000 Sebastien Masson (smasson@lodyc.jussieu.fr)
+;       - 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, X, Y, M, nT, Double = Double, ZERO2NAN = zero2nan
-;
-   if double then one = 1.0d ELSE one = 1.0
-; Sample cross covariance function.
-   TimeDim = size(X, /n_dimensions)
-   Xmean = TOTAL(X, TimeDim, Double = Double) / nT
-   Xmean = Xmean[*]#replicate(one,  nT - M)
-   Ymean = TOTAL(Y, TimeDim, Double = Double) / nT
-   Ymean = Ymean[*]#replicate(one,  nT - M)
-;
-   case TimeDim of
-      1:res = TOTAL((X[0:nT - M - 1L] - Xmean) * (Y[M:nT - 1L] - Ymean) $
+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((X[*, 0:nT - M - 1L] - Xmean) * (Y[*, M:nT - 1L] - Ymean) $
-                    , TimeDim, Double = Double)
-      3:res = TOTAL((X[*, *, 0:nT - M - 1L] - Xmean) * (Y[*, *, M:nT - 1L] - Ymean) $
-                    , TimeDim, Double = Double)
-      4:res = TOTAL((X[*, *, *, 0:nT - M - 1L] - Xmean) * (Y[*, *, *, M:nT - 1L] - Ymean) $
-                    , TimeDim, 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
@@ -157,64 +155,59 @@
    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 Correl = FLTARR(nLag) else Correl = DBLARR(nLag)
-      2:if Double eq 0 then Correl = FLTARR(Xsize[1], nLag) else Correl = DBLARR(Xsize[1], nLag)
-      3:if Double eq 0 then Correl = FLTARR(Xsize[1], Xsize[2], nLag) $
-      else Correl = DBLARR(Xsize[1], Xsize[2], nLag)
-      4:if Double eq 0 then Correl = FLTARR(Xsize[1], Xsize[2], Xsize[3], nLag) $
-      else Correl = DBLARR(Xsize[1], Xsize[2], Xsize[3], nLag)
+      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 Correlation.
+   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:Correl[k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / $
-                sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $
-                     TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan))
-               2:Correl[*, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / $
-                sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $
-                     TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan))
-               3:Correl[*, *, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / $
-                sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $
-                     TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan))
-               4:Correl[*, *, *, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / $
-                sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $
-                     TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan))
-            endcase
+               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:Correl[k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / $
-                sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $
-                     TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan))
-               2:Correl[*, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / $
-                sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $
-                     TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan))
-               3:Correl[*, *, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / $
-                sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $
-                     TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan))
-               4:Correl[*, *, *, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / $
-                sqrt(TimeCross_Cov(X,X, 0L, nT, Double = Double, /zero2nan) * $
-                     TimeCross_Cov(Y,Y, 0L, nT, Double = Double, /zero2nan))
-            endcase
+               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
+       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:Correl[k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / nT
-               2:Correl[*, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / nT
-               3:Correl[*, *, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / nT
-               4:Correl[*, *, *, k] = TimeCross_Cov(X, Y, Lag[k], nT, Double = Double) / nT
+               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:Correl[k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / nT
-               2:Correl[*, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / nT
-               3:Correl[*, *, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / nT
-               4:Correl[*, *, *, k] = TimeCross_Cov(y, x, ABS(Lag[k]), nT, Double = Double) / nT
+               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 
@@ -222,6 +215,5 @@
    endelse
 
-   if Double eq 0 then RETURN, FLOAT(Correl) else $
-    RETURN, Correl
+   if Double eq 0 then RETURN, FLOAT(Cross) else RETURN,  Cross
 
 END
