;+
;
; @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
; <pro>spl_incr</pro> 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 example to build daily data
; from monthly mean and keep the monthly mean of the computed daily
; data equal to the original values)
;
; @param x {in}{required}
; An n-elements (at least 2) input vector that specifies the tabulate points in
; a strict ascending order.
;
; @param yin {in}{required}{type=array}
; 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 mean 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:
;
; IDL> yr1 = 1990
; IDL> yr2 = 1992
; IDL> nyr = yr2-yr1+1
; IDL> n1 = 12*nyr+1
; IDL> x = julday(1+findgen(n1), replicate(1, n1) $
; IDL>        , replicate(yr1, n1), fltarr(n1))
; IDL> n2 = 365*nyr + total(leapyr(yr1+indgen(nyr))) + 1
; IDL> x2 = julday(replicate(1, n2), 1+findgen(n2) $
; IDL>            , replicate(yr1, n2), fltarr(n2))
; IDL> y = abs(randomn(0, n1-1))
; IDL> y2 = spl_keep_mean(x, y, x2, /ge0)

; IDL> print, min(x, max = ma), ma
; IDL> print, min(x2, max = ma), ma
; IDL> print, vairdate([min(x, max = ma), ma])
; IDL> print, total(y*(x[1:n1-1]-x[0:n1-2]))
; IDL> print, total(y2*(x2[1:n2-1]-x2[0:n2-2]))
;
; @history
;  Sebastien Masson (smasson\@lodyc.jussieu.fr): May 2005
;
; @version
; $Id$
;
;-
FUNCTION spl_keep_mean, x, yin, x2, YP0 = yp0, YPN_1 = ypn_1, GE0 = ge0
;
  compile_opt idl2, strictarrsubs
;
;---------------------------------
; check and initialization ...
;---------------------------------
;
  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