;+
; 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.
;-


pro quinterp, xtab, ytab, yptab, ypptab, xint, yint, $
              ypint=ypint, yppint=yppint, $
              missing=missing0

  ntab = n_elements(xtab)

  if n_elements(xtab) EQ 0 OR n_elements(ytab) EQ 0 then begin
      message, 'ERROR: XTAB and YTAB must be passed'
  endif
  if (n_elements(xtab) NE n_elements(ytab) OR $
      n_elements(xtab) NE n_elements(yptab) OR $
      n_elements(xtab) NE n_elements(ypptab)) then begin
      message, $
        'ERROR: Number of elements of XTAB, YTAB, YPTAB and YPPTAB must agree'
  endif
  if n_elements(xint) EQ 0 then begin
      message, 'ERROR: XINT must be passed'
  endif
  if n_elements(missing0) EQ 0 then begin
      miss = 0d
  endif else begin
      miss = missing0(0)
  endelse

  ;; Locate previous tabulated value
  ii = value_locate(xtab, xint)

  ;; Here we make a safety check, in case the desired point(s) is
  ;; above or below the interior of the interpolation range.  In that
  ;; case, we will need to extrapolate, based on the next nearest
  ;; interval.
  iis = ii > 0 < (ntab-2)
  whbad = where(xint LT xtab(0) OR xint GT xtab(ntab-1), ctbad)

  ;; Distance from interpolated abcissae to previous tabulated abcissa
  dx = (xint - xtab(iis))
  
  ;; Distance between adjoining tabulated abcissae and ordinates
  xs = xtab(iis+1) - xtab(iis)
  ys = ytab(iis+1) - ytab(iis) 

  ;; Rescale or pull out quantities of interest
  dx   = dx/xs              ;; Rescale DX
  y0   = ytab(iis)          ;; No rescaling of Y - start of interval
  y1   = ytab(iis+1)        ;; No rescaling of Y - end of interval
  yp0  = yptab(iis)*xs      ;; Rescale tabulated derivatives - start of interval
  yp1  = yptab(iis+1)*xs    ;; Rescale tabulated derivatives - end of interval
  ypp0 = ypptab(iis)*xs*xs  ;; Rescale tabulated 2nd der. - start of interval
  ypp1 = ypptab(iis+1)*xs*xs;; Rescale tabulated 2nd der. - end of interval


  ;; Compute values of t^n for quintic (n = 0 .. 5)
  t0 = 1d
  t1 = dx
  t2 = dx*dx
  t3 = dx*t2
  t4 = dx*t3
  t5 = dx*t4

  ;; Quintic Hermite polynomial
  yint = ((-6*t5 + 15*t4 - 10*t3 + 1 )*y0 + $
          ( 6*t5 - 15*t4 + 10*t3     )*y1 + $
          (-3*t5 +  8*t4 -  6*t3 + t1)*yp0 + $
          (-3*t5 +  7*t4 -  4*t3     )*yp1 + $
          (  -t5 +  3*t4 -  3*t3 + t2)*ypp0/2d + $
          (   t5 -  2*t4 +    t3     )*ypp1/2d)
  if ctbad GT 0 then yint(whbad) = miss
  
  if arg_present(ypint) then begin
      ypint = ((-30*t4 + 60*t3 - 30*t2       )*y0 + $
               ( 30*t4 - 60*t3 + 30*t2       )*y1 + $
               (-15*t4 + 32*t3 - 18*t2 + 1   )*yp0 + $
               (-15*t4 + 28*t3 - 12*t2       )*yp1 + $
               ( -5*t4 + 12*t3 -  9*t2 + 2*t1)*ypp0/2d + $
               (  5*t4 -  8*t3 +  3*t2       )*ypp1/2d) / xs
      if ctbad GT 0 then ypint(whbad) = miss
  endif
  
  if arg_present(yppint) then begin
      yppint = ((-120*t3 + 180*t2 - 60*t1     )*y0 + $
                ( 120*t3 - 180*t2 + 60*t1     )*y1 + $
                ( -60*t3 +  96*t2 - 36*t1     )*yp0 + $
                ( -60*t3 +  84*t2 - 24*t1     )*yp1 + $
                ( -20*t3 +  36*t2 - 18*t1 + 2d)*ypp0/2d + $
                (  20*t3 -  24*t2 +  6*t1     )*ypp1/2d) / (xs*xs)
      if ctbad GT 0 then yppint(whbad) = miss
  endif

  return
end