source: trunk/SRC/Interpolation/spl_fstdrv.pro @ 232

;+
;
; @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
; IDL> 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
;
; @version
; $Id$
;
;-
;
FUNCTION spl_fstdrv, x, y, yscd, x2
;
; compute the first derivative of the spline function
;
  compile_opt idl2, strictarrsubs
;
  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