;+
;
; @file_comments
; 
; @categories
; Statistics
;
; @param X {in}{required}{type=array}
; An array which last dimension is the time dimension so size n.
;
; @param NT
;
; @keyword DOUBLE
; If set to a non-zero value, computations are done in
; double precision arithmetic.
;
; @keyword NAN
;
; @hidden
;
; @version
; $Id$
;
;-
FUNCTION skewness_num, x, nt, DOUBLE = double,  NAN =  nan
; Compute the numerator of the skewness expression
;

  compile_opt idl2, strictarrsubs
;  
  TimeDim = size(X, /n_dimensions)
  Xmean = NAN ? TOTAL(X, TimeDim, Double = Double,  NAN =  nan) / TOTAL(FINITE(X), TimeDim) : $
   TOTAL(X, TimeDim, Double = Double) / nT    
  one = double ? 1.0d : 1.0
  Xmean = Xmean[*]#replicate(one, nT)
  res = TOTAL( (X-Xmean)^3, TimeDim, Double = Double, NAN =  nan)

  RETURN, res

END
;+
; @file_comments
; Same function as SKEWNESS but accept array (until 4
; dimension) for input and perform  the skewness 
; along the time dimension which must be the last
; one of the input array.
;
; @categories
; Statistics
;
; @param X {in}{required}{type=array}
; An array which last dimension is the time dimension so size n.
;
; @keyword DOUBLE
; If set to a non-zero value, computations are done in
; double precision arithmetic.
;
; @keyword NVAL
; A named variable that, on exit, contains the number of valid
; observations (not NAN)
;
; @examples
;
; @history
; May 2007 Michel Kolasinski (michel.kolasinski@locean-ipsl.upmc.fr)
;
; Based on the  a_timecorrelate procedure of IDL
; INTRODUCTION TO STATISTICAL TIME SERIES
; Wayne A. Fuller
; ISBN 0-471-28715-6
;
; @version
; $Id$
;
;-

FUNCTION skewness_4d, x, DOUBLE = double,  NVAL =  nval
;
   compile_opt idl2, strictarrsubs
;
; Compute the skewness from 1d to 4d vectors

   ON_ERROR, 2

   XDim = SIZE(X, /dimensions)
   XNDim = SIZE(X, /n_dimensions)
   nT = XDim[XNDim-1]

; Keyword NAN activated if needed
; Keyword NVAL not compulsory.

   NAN = ( (WHERE(FINITE(X) EQ 0 ))[0] NE -1 ) ? 1 : 0
;We can retrieve the matrix of real lenghts of time-series   
   nTreal = ( (WHERE(FINITE(X) EQ 0 ))[0] NE -1 ) ?  TOTAL(FINITE(X),  XNDim) : nT 

   IF ARG_PRESENT(NVAL) THEN nval = nTreal

; Check length.
   IF (WHERE(nTreal LE 1))[0] NE -1 THEN $
    MESSAGE,  "Matrix of length of time-series 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.
   type = SIZE(X, /TYPE)
   useDouble = (N_ELEMENTS(Double) eq 1) ? KEYWORD_SET(Double) : $
    (type eq 5)

; Type of outputs according to the type of data input
   case XNDim of
      1: Skew = useDouble ? DBLARR(1) : FLTARR(1)
      2: Skew = useDouble ? DBLARR(XDim[0]) : FLTARR(XDim[0])
      3: Skew = useDouble ? DBLARR(XDim[0], XDim[1]) : FLTARR(XDim[0], XDim[1])
      4: Skew = useDouble ? DBLARR(XDim[0], XDim[1], XDim[2]) : FLTARR(XDim[0], XDim[1], XDim[2])
   endcase
; Compute standard deviation
; nTreal might be a matrix
   std = a_timecorrelate(X, 0, /covariance)
   std = sqrt(std)
   zero = where(std EQ 0)

   if zero[0] NE -1 then STOP,  $
    'Cannot compute skewness since there are zeros in the matrix of standard deviations !'

; Problem with high masked values (x^3 makes NaN when x is high)
; Threshold put on the values of the tab
   idx_std = WHERE (std GT 1.0e+10)
   X = X < 1.0e+10
   std = std < 1.0e+10

; Compute skewness
   Skew = Skewness_Num(X, nT, Double = useDouble, NAN = nan) / (nTreal*std^3)
   IF idx_std[0] NE -1 THEN Skew[idx_std] = valmask

   return, useDouble ? Skew : FLOAT(Skew)

END