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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
# CV_BCHEB -- Procedure to evaluate all the non-zero Chebyshev functions for
# a set of points and given order.
procedure $tcv_bcheb (x, npts, order, k1, k2, basis)
PIXEL x[npts] # array of data points
int npts # number of points
int order # order of polynomial, order = 1, constant
PIXEL k1, k2 # normalizing constants
PIXEL basis[ARB] # basis functions
int k, bptr
begin
bptr = 1
do k = 1, order {
if (k == 1)
call amovk$t (PIXEL(1.0), basis, npts)
else if (k == 2)
call alta$t (x, basis[bptr], npts, k1, k2)
else {
call amul$t (basis[1+npts], basis[bptr-npts], basis[bptr],
npts)
call amulk$t (basis[bptr], PIXEL(2.0), basis[bptr], npts)
call asub$t (basis[bptr], basis[bptr-2*npts], basis[bptr], npts)
}
bptr = bptr + npts
}
end
# CV_BLEG -- Procedure to evaluate all the non zero Legendre function
# for a given order and set of points.
procedure $tcv_bleg (x, npts, order, k1, k2, basis)
PIXEL x[npts] # number of data points
int npts # number of points
int order # order of polynomial, 1 is a constant
PIXEL k1, k2 # normalizing constants
PIXEL basis[ARB] # array of basis functions
int k, bptr
PIXEL ri, ri1, ri2
begin
bptr = 1
do k = 1, order {
if (k == 1)
call amovk$t (PIXEL(1.0), basis, npts)
else if (k == 2)
call alta$t (x, basis[bptr], npts, k1, k2)
else {
ri = k
ri1 = (PIXEL(2.0) * ri - PIXEL(3.0)) / (ri - PIXEL(1.0))
ri2 = - (ri - PIXEL(2.0)) / (ri - PIXEL(1.0))
call amul$t (basis[1+npts], basis[bptr-npts], basis[bptr],
npts)
call awsu$t (basis[bptr], basis[bptr-2*npts],
basis[bptr], npts, ri1, ri2)
}
bptr = bptr + npts
}
end
# CV_BSPLINE1 -- Evaluate all the non-zero spline1 functions for a set
# of points.
procedure $tcv_bspline1 (x, npts, npieces, k1, k2, basis, left)
PIXEL x[npts] # set of data points
int npts # number of points
int npieces # number of polynomial pieces minus 1
PIXEL k1, k2 # normalizing constants
PIXEL basis[ARB] # basis functions
int left[ARB] # indices of the appropriate spline functions
int k
begin
call alta$t (x, basis[1+npts], npts, k1, k2)
call acht$ti (basis[1+npts], left, npts)
call aminki (left, npieces, left, npts)
do k = 1, npts {
basis[npts+k] = max (PIXEL(0.0), min (PIXEL(1.0),
basis[npts+k] - left[k]))
basis[k] = max (PIXEL(0.0), min (PIXEL(1.0), PIXEL(1.0) -
basis[npts+k]))
}
end
# CV_BSPLINE3 -- Procedure to evaluate all the non-zero basis functions
# for a cubic spline.
procedure $tcv_bspline3 (x, npts, npieces, k1, k2, basis, left)
PIXEL x[npts] # array of data points
int npts # number of data points
int npieces # number of polynomial pieces minus 1
PIXEL k1, k2 # normalizing constants
PIXEL basis[ARB] # array of basis functions
int left[ARB] # array of indices for first non-zero spline
int i
pointer sp, sx, tx
PIXEL dsx, dtx
begin
# allocate space
call smark (sp)
call salloc (sx, npts, TY_PIXEL)
call salloc (tx, npts, TY_PIXEL)
# calculate the index of the first non-zero coeff
call alta$t (x, Mem$t[sx], npts, k1, k2)
call acht$ti (Mem$t[sx], left, npts)
call aminki (left, npieces, left, npts)
do i = 1, npts {
Mem$t[sx+i-1] = max (PIXEL(0.0), min (PIXEL(1.0),
Mem$t[sx+i-1] - left[i]))
Mem$t[tx+i-1] = max (PIXEL(0.0), min (PIXEL(1.0), PIXEL(1.0) -
Mem$t[sx+i-1]))
}
# calculate the basis function
#call apowk$t (Mem$t[tx], 3, basis, npts)
do i = 1, npts {
dsx = Mem$t[sx+i-1]
dtx = Mem$t[tx+i-1]
basis[i] = dtx * dtx * dtx
basis[npts+i] = PIXEL(1.0) + dtx * (PIXEL(3.0) + dtx *
(PIXEL(3.0) - PIXEL(3.0) * dtx))
basis[2*npts+i] = PIXEL(1.0) + dsx * (PIXEL(3.0) + dsx *
(PIXEL(3.0) - PIXEL(3.0) * dsx))
basis[3*npts+i] = dsx * dsx * dsx
}
#call apowk$t (Mem$t[sx], 3, basis[1+3*npts], npts)
# release space
call sfree (sp)
end
|
