/usr/share/gnudatalanguage/astrolib/polint.pro is in gdl-astrolib 2018.02.16+dfsg-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 | 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
|