pro polint, xa, ya, x, y, dy ;+ ; NAME: ; POLINT ; PURPOSE: ; Interpolate a set of N points by fitting a polynomial of degree N-1 ; EXPLANATION: ; Adapted from algorithm in Numerical Recipes, Press et al. (1992), ; Section 3.1. ; ; CALLING SEQUENCE ; POLINT, xa, ya, x, y, [ dy ] ; INPUTS: ; XA - X Numeric vector, all values must be distinct. The number of ; values in XA should rarely exceed 10 (i.e. a 9th order polynomial) ; YA - Y Numeric vector, same number of elements ; X - Numeric scalar specifying value to be interpolated ; ; OUTPUT: ; Y - Scalar, interpolated value in (XA,YA) corresponding to X ; ; OPTIONAL OUTPUT ; DY - Error estimate on Y, scalar ; ; EXAMPLE: ; Find sin(2.5) by polynomial interpolation on sin(indgen(10)) ; ; IDL> xa = indgen(10) ; IDL> ya = sin( xa ) ; IDL> polint, xa, ya, 2.5, y ,dy ; ; The above method gives y = .5988 & dy = 3.1e-4 a close ; approximation to the actual sin(2.5) = .5985 ; ; METHOD: ; Uses Neville's algorithm to iteratively build up the correct ; polynomial, with each iteration containing one higher order. ; ; REVISION HISTORY: ; Written W. Landsman January, 1992 ; Converted to IDL V5.0 W. Landsman September 1997 ;- On_error,2 if N_params() LT 4 then begin print,'Syntax - polint, xa, ya, x, y, [ dy ]' print,' xa,ya - Input vectors to be interpolated' print,' x - Scalar specifying point at which to interpolate' print,' y - Output interpolated scalar value' print,' dy - Optional error estimate on y' return endif N = N_elements( xa ) if N_elements( ya ) NE N then message, \$ 'ERROR - Input X and Y vectors must have same number of elements' ; Find the index of XA which is closest to X dif = min( abs(x-xa), ns ) c = ya & d = ya y = ya[ns] ns = ns - 1 for m = 1,n-1 do begin ho = xa[0:n-m-1] - x hp = xa[m:n-1] - x w = c[1:n-m] - d[0:n-m-1] den = ho - hp if min( abs(den) ) EQ 0 then message, \$ 'ERROR - All input X vector values must be distinct' den = w / den d = hp * den c = ho * den if ( 2*ns LT n-m-1 ) then dy = c[ns+1] else begin dy = d[ns] ns = ns - 1 endelse y = y + dy endfor return end