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
|
subroutine tuiep3 (xa, ya, x, y, dy)
C
C Evaluate a cubic polynomial interpolation function at X.
C xa(klo) <= x <= xa(klo+1)
C This routine was copied with slight modifications from the POLINT
C subroutine in Numerical Recipes by Press, Flannery, Teukolsky and
C Vetterling.
C
C XA i: array of four independent-variable values
C YA i: array of four dependent-variable values
C X i: the independent variable
C Y o: the result of interpolations
C DY o: an estimate of the error of interpolation
C
CH Phil Hodge, 18-Apr-1988 Subroutine copied from Numerical Recipes POLINT.
C
double precision xa(4), ya(4), x, y, dy
C--
integer n
C four-point interpolation
parameter (n = 4)
double precision c(n), d(n), dif, dift, den
double precision ho, hp, w
integer i, m, ns
do 10 i = 1, n
if (x .eq. xa(i)) then
y = ya(i)
return
endif
10 continue
ns = 1
dif = abs (x - xa(1))
do 20 i = 1, n
dift = abs (x - xa(i))
if (dift .lt. dif) then
ns = i
dif = dift
endif
c(i) = ya(i)
d(i) = ya(i)
20 continue
y = ya(ns)
ns = ns - 1
do 40 m = 1, n-1
do 30 i = 1, n-m
ho = xa(i) - x
hp = xa(i+m) - x
w = c(i+1) - d(i)
den = w / (ho - hp)
d(i) = hp * den
c(i) = ho * den
30 continue
if (2*ns .lt. n-m) then
dy = c(ns+1)
else
dy = d(ns)
ns = ns - 1
endif
y = y + dy
40 continue
return
end
|
