/usr/share/gnudatalanguage/astrolib/hermite.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 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 | function hermite,xx,ff,x, FDERIV = fderiv
;+
; NAME:
; HERMITE
; PURPOSE:
; To compute Hermite spline interpolation of a tabulated function.
; EXPLANATION:
; Hermite interpolation computes the cubic polynomial that agrees with
; the tabulated function and its derivative at the two nearest
; tabulated points. It may be preferable to Lagrangian interpolation
; (QUADTERP) when either (1) the first derivatives are known, or (2)
; one desires continuity of the first derivative of the interpolated
; values. HERMITE() will numerically compute the necessary
; derivatives, if they are not supplied.
;
; CALLING SEQUENCE:
; F = HERMITE( XX, FF, X, [ FDERIV = ])
;
; INPUT PARAMETERS:
; XX - Vector giving tabulated X values of function to be interpolated
; Must be either monotonic increasing or decreasing
; FF - Tabulated values of function, same number of elements as X
; X - Scalar or vector giving the X values at which to interpolate
;
; OPTIONAL INPUT KEYWORD:
; FDERIV - function derivative values computed at XX. If not supplied,
; then HERMITE() will compute the derivatives numerically.
; The FDERIV keyword is useful either when (1) the derivative
; values are (somehow) known to better accuracy than can be
; computed numerically, or (2) when HERMITE() is called repeatedly
; with the same tabulated function, so that the derivatives
; need be computed only once.
;
; OUTPUT PARAMETER:
; F - Interpolated values of function, same number of points as X
;
; EXAMPLE:
; Interpolate the function 1/x at x = 0.45 using tabulated values
; with a spacing of 0.1
;
; IDL> x = findgen(20)*0.1 + 0.1
; IDL> y = 1/x
; IDL> print,hermite(x,y,0.45)
; This gives 2.2188 compared to the true value 1/0.45 = 2.2222
;
; IDL> yprime = -1/x^2 ;But in this case we know the first derivatives
; IDL> print,hermite(x,y,0.45,fderiv = yprime)
; == 2.2219 ;and so can get a more accurate interpolation
; NOTES:
; The algorithm here is based on the FORTRAN code discussed by
; Hill, G. 1982, Publ Dom. Astrophys. Obs., 16, 67. The original
; FORTRAN source is U.S. Airforce. Surveys in Geophysics No 272.
;
; HERMITE() will return an error if one tries to interpolate any values
; outside of the range of the input table XX
; PROCEDURES CALLED:
; None
; REVISION HISTORY:
; Written, B. Dorman (GSFC) Oct 1993, revised April 1996
; Added FDERIV keyword, W. Landsman (HSTX) April 1996
; Test for out of range values W. Landsman (HSTX) May 1996
; Converted to IDL V5.0 W. Landsman September 1997
; Use VALUE_LOCATE instead of TABINV W. Landsman February 2001
;-
On_error,2
if N_Params() LT 3 then begin
print,'Syntax: f = HERMITE( xx, ff, x, [FDERIV = ] )'
return,0
endif
n = N_elements(xx) ;Number of knot points
m = N_elements(x) ;Number of points at which to interpolate
l = value_locate(xx,x) ;Integer index of interpolation points
bad = where( (l LT 0) or (l EQ n-1), Nbad)
if Nbad GT 0 then message, 'ERROR - Valid interpolation range is ' + $
strtrim(xx[0],2) + ' to ' + strtrim(xx[n-1],2)
n1 = n - 1
n2 = n - 2
l1 = l + 1
l2 = l1 + 1
lm1 = l - 1
h1 = double(1./(xx[l] - xx[l1]))
h2 = - h1
; If derivatives were not supplied, then compute numeric derivatives at the
; two closest knot points
if N_elements(fderiv) NE 0 then begin
f2 = fderiv[l1]
f1 = fderiv[l]
endif else begin
f1 = dblarr(m)
f2 = dblarr(m)
for i = 0,m-1 do begin
if l[i] ne 0 then begin
if l[i] lt n2 then begin
f2[i] = (ff[l2[i]] - ff[l[i]])/(xx[l2[i]]-xx[l[i]])
endif else begin
f2[i] = (ff[n1] - ff[n2])/(xx[n1] - xx[n2])
endelse
f1[i] = ( ff[l1[i]] - ff[lm1[i]] )/( xx[l1[i]] - xx[lm1[i]] )
endif else begin
f1[i] = (ff[1] - ff[0])/(xx[1] - xx[0])
f2[i] = (ff[2] - ff[0])/(xx[2] - xx[0])
endelse
endfor
endelse
xl1 = x - xx[l1]
xl = x - xx[l]
s1 = xl1*h1
s2 = xl*h2
; Now finally the Hermite interpolation formula
f = (ff[l]*(1.-2.*h1*xl) + f1*xl)*s1*s1 + $
(ff[l1]*(1.-2.*h2*xl1) + f2*xl1)*s2*s2
if m eq 1 then return,f[0] else return,f
end
|