;+
;
; @file_comments
; 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
;
; IDL> y2 = spl_fstdrv(x, y, yscd, x2)
;
; @param x {in}{required}
; An n-elements (at least 2) input vector that specifies the
; tabulate points in ascending order.
;
; @param y {in}{required}
; f(x) = y. An n-elements 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} {type= scalar or array}
; The input values for which the first derivative values are desired.
;
; @returns
; y2: f'(x2) = y2.
;
; @history
; Sebastien Masson (smasson\@lodyc.jussieu.fr): May 2005
;
; @version
; $Id$
;
;-
FUNCTION spl_fstdrv, x, y, yscd, x2
;
compile_opt idl2, strictarrsubs
;
; 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