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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
.help
arival -- random interpolator returns value
No error traps are included -- unreasonable values for the number of
data points or the x position will either produce hard errors or garbage.
This version is written in-line for speed.
.endhelp
real procedure arival(x, datain, n, interptype)
include "interpdef.h"
real x # need 1 <= x <= n
real datain[ARB] # data values
int n # number of data values
int interptype
int i, k, nx, px
real a[6], cd20, cd21, cd40, cd41, s, t, h
real bcoeff[SPLPTS+2], temp[SPLPTS+2], pc[4]
begin
switch (interptype) {
case IT_NEAREST :
return(datain[int(x + 0.5)])
case IT_LINEAR :
nx = x
# protect against x = n case
# else it will reference datain[n + 1]
if(nx >= n)
nx = nx - 1
return((x - nx) * datain[nx+1] + (nx + 1 - x) * datain[nx])
case IT_POLY3 :
nx = x
# The major complication is that near the edge interior polynomial
# must somehow be defined.
k = 0
for(i = nx - 1; i <= nx + 2; i = i + 1){
k = k + 1
# project data points into temporary array
if ( i < 1 )
a[k] = 2. * datain[1] - datain[2 - i]
else if ( i > n )
a[k] = 2. * datain[n] - datain[2 * n - i]
else
a[k] = datain[i]
}
s = x - nx
t = 1. - s
# second central differences
cd20 = 1./6. * (a[3] - 2. * a[2] + a[1])
cd21 = 1./6. * (a[4] - 2. * a[3] + a[2])
return( s * (a[3] + (s * s - 1.) * cd21) +
t * (a[2] + (t * t - 1.) * cd20) )
case IT_POLY5 :
nx = x
# The major complication is that near the edge interior polynomial
# must somehow be defined.
k = 0
for(i = nx - 2; i <= nx + 3; i = i + 1){
k = k + 1
# project data points into temporary array
if ( i < 1 )
a[k] = 2. * datain[1] - datain[2 - i]
else if ( i > n )
a[k] = 2. * datain[n] - datain[2 * n - i]
else
a[k] = datain[i]
}
s = x - nx
t = 1. - s
# second central differences
cd20 = 1./6. * (a[4] - 2. * a[3] + a[2])
cd21 = 1./6. * (a[5] - 2. * a[4] + a[3])
# fourth central differences
cd40 = 1./120. * (a[1] - 4. * a[2] + 6. * a[3] - 4. * a[4] + a[5])
cd41 = 1./120. * (a[2] - 4. * a[3] + 6. * a[4] - 4. * a[5] + a[6])
return( s * (a[4] + (s * s - 1.) * (cd21 + (s * s - 4.) * cd41)) +
t * (a[3] + (t * t - 1.) * (cd20 + (t * t - 4.) * cd40)) )
case IT_SPLINE3 :
nx = x
h = x - nx
k = 0
# maximum number of points used is SPLPTS
for(i = nx - SPLPTS/2 + 1; i <= nx + SPLPTS/2; i = i + 1){
if(i < 1 || i > n)
;
else {
k = k + 1
if(k == 1)
px = nx - i + 1
bcoeff[k+1] = datain[i]
}
}
bcoeff[1] = 0.
bcoeff[k+2] = 0.
# Use special routine for cardinal splines.
call iif_spline(bcoeff, temp, k)
px = px + 1
pc[1] = bcoeff[px-1] + 4. * bcoeff[px] + bcoeff[px+1]
pc[2] = 3. * (bcoeff[px+1] - bcoeff[px-1])
pc[3] = 3. * (bcoeff[px-1] - 2. * bcoeff[px] + bcoeff[px+1])
pc[4] = -bcoeff[px-1] + 3. * bcoeff[px] - 3. * bcoeff[px+1] +
bcoeff[px+2]
return(pc[1] + h * (pc[2] + h * (pc[3] + h * pc[4])))
}
end
|
