;+
; NAME:
;   QUINTERP
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE: 
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Quintic spline interpolation from tabulated first and second derivatives
;
; MAJOR TOPICS:
;   Interpolation, Hermite Interpolation
;
; CALLING SEQUENCE:
;   QINTERP, XTAB, YTAB, YPTAB, YPPTAB, $
;      XINT, YINT, YPINT=, YPPINT=, MISSING=
;
; DESCRIPTION:
;
;   QUINTERP performs quintic spline interpolation of a function.
;   This routine is a natural extension of CUBETERP, in that it meant
;   for interpolation where the tabulated function has known values,
;   first derivatives *and* second derivatives at each point.  Given
;   that there are six known values for each interpolation interval,
;   the resulting interpolation function is a quintic polynomial (one
;   of a class of Hermite interpolating splines).
;
;   The user provides a tabulated set of data, whose (X,Y) positions
;   are (XTAB, YTAB), and whose first and second derivatives are YPTAB
;   and YPPTAB.  The user also provides a set of desired "X" abcissae
;   for which interpolants are requested.  The interpolated spline
;   values are returned in YINT.  The interpolated curve will smoothly
;   pass through the control points, and have the requested
;   derivatives at those points.
;
;   Note that the user must provide both derivatives (they are not
;   optional).  If you don't have one or more derivatives, then you
;   should use the IDL spline functions SPL_INIT/SPL_INTERP, or the
;   functions CUBETERP, QUADTERP or LINTERP instead.  Unlike CUBETERP,
;   if the requested point is outside of the tabulated range, the
;   function is not extrapolated.  Instead the value provided by the
;   MISSING keyword is returned for those points.
;
;   The user may also optionally request the first and second
;   derivatives of the function with the YPINT and YPPINT keywords.
;
; INPUTS:
;
;   XTAB - tabulated X values.  Must be sorted in increasing order.
;
;   YTAB - tabulated Y values.
;
;   YPTAB - tabulated first derivatives ( = dY/dX ). Not optional
;   YPPTAB - tabulated second derivatives ( = d(YPTAB)/dX ).  Not optional.
;
;   XINT - X values of desired interpolants.
;
; OUTPUTS:
;
;   YINT - Y values of desired interpolants.
;
; OPTIONAL KEYWORDS:
;
;   YPINT - upon return, the slope (first derivative) at the
;           interpolated positions.
;
;   YPPINT - upon return, the second derivative at the interpolated
;            positions.
;
;   MISSING - a value to report for "missing" data.  This function
;             does not perform extrapolation; any requested point
;             outside the range [MIN(XTAB),MAX(XTAB)] is considered
;             missing.
;             Default: 0
;
; EXAMPLE:
;
;  ;; Set up some fake data, a sinusoid
;  xtab = dindgen(101)/100d * 2d*!dpi    ; 100 points from 0 -> 2*!dpi
;  ytab = sin(xtab)    ;; values
;  yptab = cos(xtab)   ;; 1st deriv
;  ypptab = -sin(xtab) ;; 2nd deriv
;
;  ;; Interpolate to a finer grid
;  xint = dindgen(1001)/1000 * 2d*!dpi   ;; 1000 points from 0->2*!dpi
;  quinterp, xtab, ytab, yptab, ypptab, xint, yint, ypint=ypint, yppint=yppint
;
;  ;; Plot it
;  plot, xint, yint
;  oplot, xtab, ytab, psym=1, symsize=2
;  for i = 0, n_elements(xtab)-1 do $   ;; Also plot slopes
;    oplot, xtab(i)+[-0.5,0.5], ytab(i)+[-0.5,0.5]*yptab(i)
; 
;
; MODIFICATION HISTORY:
;   Written and documented, CM, 08 Oct 2008
;
;  $Id: quinterp.pro,v 1.2 2009/04/15 04:17:30 craigm Exp $
;
;-
; Copyright (C) 2008, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-