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
129
130
131
132
133
134
135
136
137
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
.help
procedures to evaluate interpolants
These procedures are seperated from the other programs in order
to make it easier to optimise these in case it becomes necessary to
obtain the fastest possible interpolants. The interpolation package
spends most of its time in these routines.
.endhelp
procedure iievne(x,y,n,a) # nearest neighbor
real x[ARB] # x values, must be within [1,n]
real y[ARB] # interpolated values returned to user
int n # number of x values
real a[ARB] # data to be interpolated
int i
begin
do i = 1, n
y[i] = a[int(x[i] + 0.5)]
return
end
procedure iievli(x,y,n,a) # linear
real x[ARB] # x values, must be within [1,n]
real y[ARB] # interpolated values returned to user
int n # number of x values
real a[ARB] # data to be interpolated
int i,nx
begin
do i = 1,n {
nx = x[i]
y[i] = (x[i] - nx) * a[nx + 1] + (nx + 1 - x[i]) * a[nx]
}
return
end
procedure iievp3(x,y,n,a) # interior third order polynomial
real x[ARB] # x values, must be within [1,n]
real y[ARB] # interpolated values returned to user
int n # number of x values
real a[ARB] # data to be interpolated from a[0] to a[n+2]
int i,nx,nxold
real s,t,cd20,cd21
begin
nxold = -1
do i = 1,n {
nx = x[i]
s = x[i] - nx
t = 1. - s
if (nx != nxold) {
# second central differences:
cd20 = 1./6. * (a[nx+1] - 2. * a[nx] + a[nx-1])
cd21 = 1./6. * (a[nx+2] - 2. * a[nx+1] + a[nx])
nxold = nx
}
y[i] = s * (a[nx+1] + (s * s - 1.) * cd21) +
t * (a[nx] + (t * t - 1.) * cd20)
}
return
end
procedure iievp5(x,y,n,a) # interior fifth order polynomial
real x[ARB] # x values, must be within [1,n]
real y[ARB] # interpolated values returned to user
int n # number of x values
real a[ARB] # data to be interpolated - from a[-1] to a[n+3]
int i,nx,nxold
real s,t,cd20,cd21,cd40,cd41
begin
nxold = -1
do i = 1,n {
nx = x[i]
s = x[i] - nx
t = 1. - s
if (nx != nxold) {
cd20 = 1./6. * (a[nx+1] - 2. * a[nx] + a[nx-1])
cd21 = 1./6. * (a[nx+2] - 2. * a[nx+1] + a[nx])
# fourth central differences
cd40 = 1./120. * (a[nx-2] - 4. * a[nx-1] +
6. * a[nx] - 4. * a[nx+1] + a[nx+2])
cd41 = 1./120. * (a[nx-1] - 4. * a[nx] +
6. * a[nx+1] - 4. * a[nx+2] + a[nx+3])
nxold = nx
}
y[i] = s * (a[nx+1] + (s * s - 1.) *
(cd21 + (s * s - 4.) * cd41)) +
t * (a[nx] + (t * t - 1.) *
(cd20 + (t * t - 4.) * cd40))
}
return
end
procedure iievs3(x,y,n,a) # cubic spline evaluator
real x[ARB] # x values, must be within [1,n]
real y[ARB] # interpolated values returned to user
int n # number of x values
real a[ARB] # basis spline coefficients - from a[0] to a[n+1]
int i,nx,nxold
real s,c0,c1,c2,c3
begin
nxold = -1
do i = 1,n {
nx = x[i]
s = x[i] - nx
if (nx != nxold) {
# convert b-spline coeff's to poly. coeff's
c0 = a[nx-1] + 4. * a[nx] + a[nx+1]
c1 = 3. * (a[nx+1] - a[nx-1])
c2 = 3. * (a[nx-1] - 2. * a[nx] + a[nx+1])
c3 = -a[nx-1] + 3. * a[nx] - 3. * a[nx+1] + a[nx+2]
nxold = nx
}
y[i] = c0 + s * (c1 + s * (c2 + s * c3) )
}
return
end
|
