;+
; NAME:
;   CUBETERP
;
; 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:
;   Cubic spline interpolation with known derivatives
;
; MAJOR TOPICS:
;   Interpolation
;
; CALLING SEQUENCE:
;   CUBETERP, XTAB, YTAB, YPTAB, XINT, YINT, YPINT=, YPPINT=, EXTRAP_ORDER=
;
; DESCRIPTION:
;
;   CUBETERP performs cubic spline interpolation of a function.  This
;   routine is different from the many other spline interpolation
;   functions for IDL in that it allows you to choose the slope of the
;   spline at each control point.  I.e. it is not forced to be a
;   "natural" spline.
;
;   The user provides a tabulated set of data, whose (X,Y) positions
;   are (XTAB, YTAB), and whose derivatives are YPTAB.  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 slopes at 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 derivatives ( = dY/dX, evaluated at XTAB).
;
;   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.
;
;   EXTRAP_ORDER - technique used to extrapolate beyond the tabulated
;                  values.  Allowed values:
;                    -1 - extrapolated points are set to NaN (not a number)
;                     0 - constant extrapolation, equal to the value
;                         at the nearest tabulated point
;                     1 - linear extrapolation, based on slope at
;                         nearest tabulated value
;                     2 - quadratic extrapolation, based on slope and
;                         second derivative at nearest tabulated value
;                     3 - cubic extrapolation.
;                  DEFAULT: 2 (quadratic extrapolation)
;
;
; EXAMPLE:
;
;  ;; Set up some fake data
;  xtab = [0D,2,5,10]
;  ytab = [2D,4,-3,-5]
;  yptab = [-1D,0.5,2.3,-4]
;
;  ;; Interpolate to a finer grid
;  xint = dindgen(1001)/100
;  cubeterp, xtab, ytab, yptab, xint, yint
;
;  ;; 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, July 2003
;   Added EXTRAP_ORDER = -1 option, CM, 15 May 2005
;   Syntax error fix, CM, 07 Mar 2007
;   Clarified documentation a bit, CM, 12 Nov 2007
;   Small documentation changes, CM, 16 Apr 2009
;
;  $Id: cubeterp.pro,v 1.8 2009/05/05 04:59:57 craigm Exp $
;
;-
; Copyright (C) 2003, 2005, 2007, 2009, 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.
;-

; Cubic interpolation with known slopes


pro cubeterp, xtab, ytab, yptab, xint, yint, $
              extrap_order=extr0, $
              ypint=ypint, yppint=yppint

  ntab = n_elements(xtab)
  if n_elements(extr0) EQ 0 then extr = 2 else extr = round(extr0(0))

  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)) then begin
      message, 'ERROR: Number of elements of XTAB, YTAB and YPTAB must agree'
  endif
  if n_elements(xint) EQ 0 then begin
      message, 'ERROR: XINT must be passed'
  endif

  ;; 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
  whll = where(ii LT 0, ctll)
  if ctll GT 0 then iis(whll) = iis(whll) + 1
  whgg = where(ii GE (ntab-1), ctgg)
  if ctgg GT 0 then iis(whgg) = iis(whgg) - 1

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

  ;; Compute polynomial coefficients
  a = y0
  b = yp0
  c = 3*ys - 2*yp0 - yp1
  d = yp0 + yp1 - 2*ys

  ;; Extrapolate only quadratically
  if extr EQ 2 then begin
      if ctll GT 0 then begin
          ;; Lower end of extrapolation
          d(whll) = 0
      endif
      if ctgg GT 0 then begin
          ;; Upper end of extrapolation
          dgg = d(whgg)
          a(whgg) = a(whgg) + dgg
          b(whgg) = b(whgg) - 3*dgg
          c(whgg) = c(whgg) + 3*dgg
          d(whgg) = 0
      endif
  endif

  ;; Extrapolate only linearly
  if extr EQ 1 then begin

      if ctll GT 0 then begin
          ;; Lower end of extrapolation
          c(whll) = 0
          d(whll) = 0
      endif
      if ctgg GT 0 then begin
          ;; Upper end of extrapolation
          dgg = d(whgg)
          cgg = c(whgg)
          a(whgg) = a(whgg) - cgg - 2*dgg
          b(whgg) = b(whgg) + 2*cgg + 3*dgg
          c(whgg) = 0
          d(whgg) = 0
      endif
  endif

  if extr EQ 0 then begin
      if ctll GT 0 then begin
          ;; Lower end of extrapolation
          b(whll) = 0
          c(whll) = 0
          d(whll) = 0
      endif
      if ctgg GT 0 then begin
          ;; Upper end of extrapolation
          a(whgg) = a(whgg) + b(whgg) + c(whgg) + d(whgg)
          b(whgg) = 0
          c(whgg) = 0
          d(whgg) = 0
      endif
  endif      

  if extr EQ -1 then begin
      sz = size(y0) & tp = sz(sz(0)+1)
      if sz EQ 4 OR sz EQ 6 then nanv = !values.f_nan $
      else nanv = !values.d_nan

      if ctll GT 0 then begin
          ;; Lower end of extrapolation
          a(whll) = nanv
      endif
      if ctgg GT 0 then begin
          ;; Upper end of extrapolation
          a(whgg) = nanv
      endif
  endif

  yint = a + dx*(b + dx*(c + dx*d))

  ;; Compute derivatives if requested
  if arg_present(ypint) then ypint = (b + dx*(2*c + dx*3*d))/xs
  if arg_present(yppint) then yppint = (2*c + 6*d*dx)/xs^2


  return
end