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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <math/iminterp.h>
include "im1interpdef.h"
# ASIDER -- Calculate nder derivatives assuming that x lands in the region
# 1 <= x <= npts.
procedure asider (asi, x, der, nder)
pointer asi # interpolant descriptor
real x[ARB] # x value
real der[ARB] # derivatives, der[1] is value der[2] is f prime
int nder # number items returned = 1 + number of derivatives
int nearx, i, j, k, nterms, nd
pointer c0ptr, n0
real deltax, accum, tmpx[2], pcoeff[MAX_NDERIVS], diff[MAX_NDERIVS]
begin
# Return zero for derivatives that are zero.
do i = 1, nder
der[i] = 0.
# Nterms is number of terms in case polynomial type.
nterms = 0
# (c0ptr + 1) is the pointer to the first data point in COEFF.
c0ptr = ASI_COEFF(asi) - 1 + ASI_OFFSET(asi)
switch (ASI_TYPE(asi)) {
case II_NEAREST:
der[1] = COEFF(c0ptr + int(x[1] + 0.5))
return
case II_LINEAR:
nearx = x[1]
der[1] = (x[1] - nearx) * COEFF(c0ptr + nearx + 1) +
(nearx + 1 - x[1]) * COEFF(c0ptr + nearx)
if (nder > 1)
der[2] = COEFF(c0ptr + nearx + 1) - COEFF(c0ptr + nearx)
return
case II_SINC, II_LSINC:
call ii_sincder (x[1], der, nder,
COEFF(ASI_COEFF(asi) + ASI_OFFSET(asi)), ASI_NCOEFF(asi),
ASI_NSINC(asi), DX)
return
case II_DRIZZLE:
if (ASI_PIXFRAC(asi) >= 1.0)
call ii_driz1 (x, der[1], 1, COEFF(ASI_COEFF(asi) +
ASI_OFFSET(asi)), ASI_BADVAL(asi))
else
call ii_driz (x, der[1], 1, COEFF(ASI_COEFF(asi) +
ASI_OFFSET(asi)), ASI_PIXFRAC(asi), ASI_BADVAL(asi))
if (nder > 1) {
deltax = x[2] - x[1]
if (deltax == 0.0)
der[2] = 0.0
else {
tmpx[1] = x[1]
tmpx[2] = (x[1] + x[2]) / 2.0
if (ASI_PIXFRAC(asi) >= 1.0)
call ii_driz1 (tmpx, accum, 1, COEFF(ASI_COEFF(asi) +
ASI_OFFSET(asi)), ASI_BADVAL(asi))
else
call ii_driz (tmpx, accum, 1, COEFF(ASI_COEFF(asi) +
ASI_OFFSET(asi)), ASI_PIXFRAC(asi), ASI_BADVAL(asi))
tmpx[1] = tmpx[2]
tmpx[2] = x[2]
if (ASI_PIXFRAC(asi) >= 1.0)
call ii_driz1 (tmpx, der[2], 1, COEFF(ASI_COEFF(asi) +
ASI_OFFSET(asi)), ASI_BADVAL(asi))
else
call ii_driz (tmpx, der[2], 1, COEFF(ASI_COEFF(asi) +
ASI_OFFSET(asi)), ASI_PIXFRAC(asi), ASI_BADVAL(asi))
der[2] = 2.0 * (der[2] - accum) / deltax
}
}
return
case II_POLY3:
nterms = 4
case II_POLY5:
nterms = 6
case II_SPLINE3:
nterms = 4
default:
call error (0, "ASIDER: Unknown interpolant type")
}
# Routines falls through to this point if the interpolant is one of
# the higher order polynomial types or a third order spline.
nearx = x[1]
n0 = c0ptr + nearx
deltax = x[1] - nearx
# Compute the number of derivatives needed.
nd = nder
if (nder > nterms)
nd = nterms
# Generate the polynomial coefficients.
if (ASI_TYPE(asi) == II_SPLINE3) {
pcoeff[1] = COEFF(n0-1) + 4. * COEFF(n0) + COEFF(n0+1)
pcoeff[2] = 3. * (COEFF(n0+1) - COEFF(n0-1))
pcoeff[3] = 3. * (COEFF(n0-1) - 2. * COEFF(n0) + COEFF(n0+1))
pcoeff[4] = -COEFF(n0-1) + 3. * COEFF(n0) - 3. * COEFF(n0+1) +
COEFF(n0+2)
# Newton's form written in line to get polynomial from data
} else {
# Load data.
do i = 1, nterms
diff[i] = COEFF(n0 - nterms/2 + i)
# Generate difference table.
do k = 1, nterms - 1
do i = 1, nterms - k
diff[i] = (diff[i+1] - diff[i]) / k
# Shift to generate polynomial coefficients.
do k = nterms, 2, -1
do i = 2,k
diff[i] = diff[i] + diff[i-1] * (k - i - nterms/2)
do i = 1,nterms
pcoeff[i] = diff[nterms + 1 - i]
}
# Compute the derivatives. As the loop progresses pcoeff contains
# coefficients of higher and higher derivatives.
do k = 1, nd {
# Evaluate using nested multiplication.
accum = pcoeff[nterms - k + 1]
do j = nterms - k, 1, -1
accum = pcoeff[j] + deltax * accum
der[k] = accum
# Differentiate polynomial.
do j = 1, nterms - k
pcoeff[j] = j * pcoeff[j + 1]
}
end
|
