function hermite,xx,ff,x, FDERIV = fderiv
;+
; NAME:
;       HERMITE
; PURPOSE:
;       To compute Hermite spline interpolation of a tabulated function.
; EXPLANATION:
;       Hermite interpolation computes the cubic polynomial that agrees with 
;       the tabulated function and its derivative at the two nearest 
;       tabulated points.   It may be preferable to Lagrangian interpolation 
;       (QUADTERP) when either (1) the first derivatives are known, or (2)
;       one desires continuity of the first derivative of the interpolated
;       values.    HERMITE() will numerically compute the necessary
;       derivatives, if they are not supplied.
;       
; CALLING SEQUENCE:
;       F = HERMITE( XX, FF, X, [ FDERIV = ])
;
; INPUT PARAMETERS:
;       XX - Vector giving tabulated X values of function to be interpolated
;               Must be either monotonic increasing or decreasing   
;       FF - Tabuluated values of function, same number of elements as X
;       X -  Scalar or vector giving the X values at which to interpolate
;
; OPTIONAL INPUT KEYWORD:
;       FDERIV - function derivative values computed at XX.    If not supplied,
;               then HERMITE() will compute the derivatives numerically.
;               The FDERIV keyword is useful either when (1) the derivative
;               values are (somehow) known to better accuracy than can be 
;               computed numerically, or (2) when HERMITE() is called repeatedly
;               with the same tabulated function, so that the derivatives
;               need be computed only once.
;
; OUTPUT PARAMETER:
;       F - Interpolated values of function, same number of points as X
;
; EXAMPLE:
;       Interpolate the function 1/x at x = 0.45 using tabulated values
;       with a spacing of 0.1
;
;       IDL> x = findgen(20)*0.1 + 0.1
;       IDL> y = 1/x
;       IDL> print,hermite(x,y,0.45)         
;               This gives 2.2188 compared to the true value 1/0.45 = 2.2222
;
;       IDL> yprime = -1/x^2      ;But in this case we know the first derivatives
;       IDL> print,hermite(x,y,0.45,fderiv = yprime)
;             == 2.2219            ;and so can get a more accurate interpolation
; NOTES:
;       The algorithm here is based on the FORTRAN code discussed by 
;       Hill, G. 1982, Publ Dom. Astrophys. Obs., 16, 67.   The original 
;       FORTRAN source is U.S. Airforce. Surveys in Geophysics No 272. 
;
;       HERMITE() will return an error if one tries to interpolate any values 
;       outside of the range of the input table XX
; PROCEDURES CALLED:
;       None
; REVISION HISTORY:
;       Written,    B. Dorman (GSFC) Oct 1993, revised April 1996
;       Added FDERIV keyword,  W. Landsman (HSTX)  April 1996
;       Test for out of range values  W. Landsman (HSTX) May 1996
;       Converted to IDL V5.0   W. Landsman   September 1997
;       Use VALUE_LOCATE instead of TABINV   W. Landsman   February 2001
;-
   On_error,2 

   if N_Params() LT 3 then begin
        print,'Syntax:  f = HERMITE( xx, ff, x, [FDERIV = ] )'
        return,0
   endif

   n = N_elements(xx)           ;Number of knot points
   m = N_elements(x)            ;Number of points at which to interpolate

   l = value_locate(xx,x)       ;Integer index of interpolation points 

   bad = where( (l LT 0) or (l EQ n-1), Nbad)
        if Nbad GT 0 then message, 'ERROR - Valid interpolation range is ' + $
        strtrim(xx[0],2) + ' to ' + strtrim(xx[n-1],2)

   n1 = n - 1
   n2 = n - 2

   l1  = l + 1  
   l2 = l1 + 1
   lm1 = l - 1
   h1 = double(1./(xx[l] - xx[l1]))
   h2 = - h1

; If derivatives were not supplied, then compute numeric derivatives at the 
; two closest knot points

 if N_elements(fderiv) NE 0 then begin
        f2 = fderiv[l1]
        f1 = fderiv[l]

 endif else begin

   f1 = dblarr(m)
   f2 = dblarr(m)
   for i = 0,m-1 do begin
        if l[i] ne 0 then begin
           if l[i] lt n2 then begin
              f2[i] = (ff[l2[i]] - ff[l[i]])/(xx[l2[i]]-xx[l[i]])
           endif else begin
              f2[i] = (ff[n1] - ff[n2])/(xx[n1] - xx[n2])
           endelse
           f1[i] = ( ff[l1[i]] - ff[lm1[i]] )/( xx[l1[i]] - xx[lm1[i]] )
        endif else begin
           f1[i] = (ff[1] - ff[0])/(xx[1] - xx[0])
           f2[i] = (ff[2] - ff[0])/(xx[2] - xx[0])
        endelse
   endfor
 endelse
       
    xl1 = x - xx[l1]
    xl  = x - xx[l]
    s1  = xl1*h1
    s2  = xl*h2

; Now finally the Hermite interpolation formula

    f   = (ff[l]*(1.-2.*h1*xl) + f1*xl)*s1*s1 + $                       
                                        (ff[l1]*(1.-2.*h2*xl1) + f2*xl1)*s2*s2

    if m eq 1 then return,f[0] else return,f

    end