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