PRO quadterp, xtab, ytab, xint, yint, MISSING = MISSING
;+
; NAME:
; QUADTERP
; PURPOSE:
; Quadratic interpolation of X,Y vectors onto a new X grid
; EXPLANATION:
; Interpolate a function Y = f(X) at specified grid points using an
; average of two neighboring 3 point quadratic (Lagrangian) interpolants.
; Use LINTERP for linear interpolation
;
; CALLING SEQUENCE:
; QUADTERP, Xtab, Ytab, Xint, Yint, [ MISSING = ]
;
; INPUT:
; Xtab - Vector (X TABle) containing the current independent variable
; Must be either monotonic increasing or decreasing
; Ytab - Vector (Y TABle) containing the dependent variable defined
; at each of the points of XTAB.
; Xint - Scalar or vector giving the values of X for which interpolated
; Y values are sought
;
; OUTPUT:
; Yint - Interpolated value(s) of Y, same number of points as Xint
;
; OPTIONAL INPUT KEYWORD:
; MISSING - Scalar specifying Yint value(s) to be assigned, when Xint
; value(s) are outside of the range of Xtab. Default is to
; truncate the out of range Yint value(s) to the nearest value
; of Ytab. See the help for the INTERPOLATE function.
; METHOD:
; 3-point Lagrangian interpolation. The average of the two quadratics
; derived from the four nearest points is returned in YTAB. A single
; quadratic is used near the end points. VALUE_LOCATE is used
; to locate center point of the interpolation.
;
; NOTES:
; QUADTERP provides one method of high-order interpolation. The
; RSI interpol.pro function includes the following alternatives:
;
; interpol(/LSQUADRATIC) - least squares quadratic fit to a 4 pt
; neighborhood
; interpol(/QUADRATIC) - quadratic fit to a 3 pt neighborhood
; interpol(/SPLINE) - cubic spline fit to a 4 pt neighborhood
;
; Also, the IDL Astro function HERMITE fits a cubic polynomial and its
; derivative to the two nearest points.
; RESTRICTIONS:
; Unless MISSING keyword is set, points outside the range of Xtab in
; which valid quadratics can be computed are returned at the value
; of the nearest end point of Ytab (i.e. Ytab[0] and Ytab[NPTS-1] ).
;
; EXAMPLE:
; A spectrum has been defined using a wavelength vector WAVE and a
; flux vector FLUX. Interpolate onto a new wavelength grid, e.g.
;
; IDL> wgrid = [1540.,1541.,1542.,1543.,1544.,1545.]
; IDL> quadterp, wave, flux, wgrid, fgrid
;
; FGRID will be a 5 element vector containing the quadratically
; interpolated values of FLUX at the wavelengths given in WGRID.
;
; EXTERNAL ROUTINES:
; ZPARCHECK
; REVISION HISTORY:
; 31 October 1986 by B. Boothman, adapted from the IUE RDAF
; 12 December 1988 J. Murthy, corrected error in Xint
; September 1992, W. Landsman, fixed problem with double precision
; August 1993, W. Landsman, added MISSING keyword
; June, 1995, W. Landsman, use single quadratic near end points
; Converted to IDL V5.0 W. Landsman September 1997
; Fix occasional problem with integer X table,
; YINT is a scalar if XINT is a scalar W. Landsman Dec 1999
; Use VALUE_LOCATE instead of TABINV W. Landsman Feb. 2000
;-
On_error,2
compile_opt idl2
if N_params() LT 4 then begin
print,'Syntax - QUADTERP, xtab, ytab, xint, yint, [ MISSING = ]'
return
endif
zparcheck,'QUADTERP',xtab,1,[1,2,3,4,5],1,'Independent (X) vector'
zparcheck,'QUADTERP',ytab,2,[1,2,3,4,5],1,'Dependent (Y) vector'
npts = min( [N_elements(xtab), N_elements(ytab) ] )
m = n_elements(xint)
if size(xtab,/TNAME) NE 'DOUBLE' then xt = float(xtab) else xt = xtab
Xmin = min( [ Xtab[0],Xtab[npts-1] ], max = Xmax)
u = xint > Xmin < Xmax
if npts LT 3 then $
message,' ERROR - At least 3 points required for quadratic interpolation'
; Determine index of data-points from which interpolation is made
index = value_locate(xtab,xint) > 0L < (npts-2)
; First quadratic
i0 = (index-1) > 0 & i1 = i0+1 & i2 = (i1 +1)
x0 = xt[i0] & x1 = xt[i1] & x2 = xt[i2]
p1 = ytab[i0] * (u-x1) * (u-x2) / ((x0-x1) * (x0-x2)) + $
ytab[i1] * (u-x0) * (u-x2) / ((x1-x0) * (x1-x2)) + $
ytab[i2] * (u-x0) * (u-x1) / ((x2-x0) * (x2-x1))
; Second Quadratic
i2 = (index+2) < (npts-1) & i1 = i2-1 & i0 = (i1-1)
x0 = xt[i0] & x1 = xt[i1] & x2 = xt[i2]
p2 = ytab[i0] * (u-x1) * (u-x2) / ((x0-x1) * (x0-x2)) + $
ytab[i1] * (u-x0) * (u-x2) / ((x1-x0) * (x1-x2)) + $
ytab[i2] * (u-x0) * (u-x1) / ((x2-x0) * (x2-x1))
yint = (p1 + p2) / 2. ;Average of two quadratics
if N_elements(missing) EQ 1 then begin
bad = where( (Xint LT Xmin) or (Xint GT Xmax ), Nbad)
if Nbad GT 0 then Yint[bad] = missing
endif
return
end