;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; @file_comments
;
; Given the arrays X and Y, which tabulate a function (with the X[i]
; AND Y[i] in ascending order), and given an input value X2, the
; SPL_INCR function returns an interpolated value for the given values
; of X2. The interpolation method is based on cubic spline, corrected
; in a way that integral of the interpolated values is the same as the
; integral of the input values. (-> for exemple to build daily data
; from monthly mean and keep the monthly mean of the computed daily
; data equa to the original values)
;
; @examples  y2 =  spl_keep_mean(x, y, x2)
;
;    @param x {in}{required}  An n-element (at least 2) input vector that specifies the
;    tabulate points in a strict ascending order.
;
;    @param yin {in}{required}  an array with one element less than x. y[i] represents the
;    mean value between x[i] and x[i+1]. if /GE0 is activated, y must
;    have positive values.
;
;    @param x2 {in}{required}  The input values for which the interpolated values are
;    desired. Its values must be strictly monotonically increasing. 
;
;
; @keyword    /GE0 to force that y2 is always GE than 0. In that case, y must
;    also be GE than 0.
;
; @keyword    YP0 The first derivative of the interpolating function at the
;    point X0. If YP0 is omitted, the second derivative at the
;    boundary is set to zero, resulting in a "natural spline."
;
; @keyword    YPN_1 The first derivative of the interpolating function at the
;    point Xn-1. If YPN_1 is omitted, the second derivative at the
;    boundary is set to zero, resulting in a "natural spline." 
;
; @returns 
;
;    y2: the meean value between two consecutive values of x2. This
;    array has one element less than y2. y2 has double precision.
;
; @restrictions
;   It might be possible that y2 has very small negative values
;   (amplitude smaller than 1.e-6)... 
;
;
; @examples 
;
;    12 monthly values of precipitations into daily values:
;
;    yr1 = 1990
;    yr2 = 1992
;    nyr = yr2-yr1+1
;    n1 = 12*nyr+1
;    x = julday(1+findgen(n1), replicate(1, n1) $
;           , replicate(yr1, n1), fltarr(n1))
;    n2 = 365*nyr + total(leapyr(yr1+indgen(nyr))) + 1
;    x2 = julday(replicate(1, n2), 1+findgen(n2) $
;               , replicate(yr1, n2), fltarr(n2))
;    y = abs(randomn(0, n1-1))
;    y2 = spl_keep_mean(x, y, x2, /ge0)

;    print, min(x, max = ma), ma
;    print, min(x2, max = ma), ma
;    print, vairdate([min(x, max = ma), ma])
;    print, total(y*(x[1:n1-1]-x[0:n1-2]))
;    print, total(y2*(x2[1:n2-1]-x2[0:n2-2]))
;
; @history
;  Sebastien Masson (smasson\@lodyc.jussieu.fr): May 2005
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION spl_keep_mean, x, yin, x2, YP0 = yp0, YPN_1 = ypn_1, GE0 = ge0
;
;---------------------------------
; check and initialisation ...
;---------------------------------
;
  nx = n_elements(x)
  ny = n_elements(yin)
  nx2 = n_elements(x2)
; x must have at least 2 elements
  IF nx LT 2 THEN stop 
; x2 must have at least 2 elements
  IF nx2 LT 2 THEN stop 
; x be monotonically increasing
  IF min(x[1:nx-1]-x[0:nx-2]) LE 0 THEN stop 
; x2 be monotonically increasing
  IF min(x2[1:nx2-1]-x2[0:nx2-2])  LE 0 THEN stop 
;
;---------------------------------
; compute the integral of y
;---------------------------------
; if spl_keep_mean is called by the user (and not by itself) we must compute
; the integral of y. yin must have one element less than x
  IF nx NE ny+1 THEN stop
  y = double(yin)*double(x[1:nx-1]-x[0:nx-2])
  y = [0.0d, temporary(y)]
  y = total(temporary(y), /cumulative, /double)
;
;---------------------------------
; compute the "spline" interpolation
;---------------------------------
;
  IF keyword_set(ge0) THEN BEGIN
; if the want that the interpolated values are always >= 0, we must
; have yin >= 0.0d
    IF min(yin) LT 0 THEN stop
; call spl_incr
    y2 = spl_incr(x, temporary(y), x2, yp0 = yp0, ypn_1 = ypn_1)
  ENDIF ELSE BEGIN
    yscd = spl_init(x, y, yp0 = yp0, ypn_1 = ypn_1, /double)
    y2 = spl_interp(x, y, temporary(yscd), x2, /double)
  ENDELSE
;                ;
;---------------------------------
; Compute the derivative of y
;---------------------------------
;
  yfrst = (y2[1:nx2-1]-y2[0:nx2-2])/(x2[1:nx2-1]-x2[0:nx2-2])
; it can happen that we have very small negative values (-1.e-6 for ex)
;  yfrst = 0.0d > temporary(yfrst)
  RETURN, yfrst
  
;------------------------------------------------------------------
;------------------------------------------------------------------
;
 END