;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
;
; @file_comments SPL_FSTDRV returns the values of the first derivative of
; the interpolating function at the points X2i. it is a double
; precision array.
;
; 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 interpolated value are also in ascending order
;
; @examples  y2 =  spl_fstdrv(x, y, yscd, x2)
;
;    @param x {in}{required}  An n-element (at least 2) input vector that specifies the
;    tabulate points in ascending order.
;
;    @param y {in}{required}  f(x) = y. An n-element input vector that specifies the values
;    of the tabulated function F(Xi) corresponding to Xi.
;
;    @param yscd {in}{required}  The output from SPL_INIT for the specified X and Y.
;
;    @param x2 {in}{required}  The input values for which the first derivative values are
;    desired. X can be scalar or an array of values.

; @returns 
;
;    y2: f'(x2) = y2. 
;
;
; @history
;  Sebastien Masson (smasson\@lodyc.jussieu.fr): May 2005
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION spl_fstdrv, x, y, yscd, x2
;
; compute the first derivative of the spline function;
;
  nx = n_elements(x)
  ny = n_elements(y)
; x must have at least 2 elements
  IF nx LT 2 THEN stop   
; y must have the same number of elements than x
  IF nx NE ny THEN stop
; define loc in a way that 
;  if loc[i] eq -1   :                 x2[i] <  x[0]
;  if loc[i] eq nx2-1:                 x2[i] >= x[nx-1]
;  else              :    x[loc[i]] <= x2[i] <  x[loc[i]+1]
  loc = value_locate(x, x2)
; change loc in order to 
; use x[0]    and x[1]    even if x2[i] <  x[0]    -> extrapolation
; use x[nx-2] and x[nx-1] even if x2[i] >= x[nx-1] -> extrapolation
  loc = 0 > temporary(loc) < (nx-2) 
; distance between to consecutive x
  deltax = x[loc+1]-x[loc]
; distance between to consecutive y
  deltay = y[loc+1]-y[loc]
; relative distance between x2[i] and x[loc[i]+1]
  a = (x[loc+1]-x2)/deltax
; relative distance between x2[i] and x[loc[i]]
  b = 1.0d - a
; compute the first derivative on x (see numerical recipes Chap 3.3)
  yfrst = temporary(deltay)/deltax $
    - 1.0d/6.0d * (3.0d*a*a - 1.0d) * deltax * yscd[loc] $
    + 1.0d/6.0d * (3.0d*b*b - 1.0d) * deltax * yscd[loc+1]
; beware of the computation precision... 
; force near zero values to be exactly 0.0
  zero = where(abs(yfrst) LT 1.e-10)
  IF zero[0] NE -1 THEN yfrst[zero] = 0.0d

  RETURN, yfrst
END