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 dcv_bcheb (x, npts, order, k1, k2, basis)
double x[npts] # array of data points
int npts # number of points
int order # order of polynomial, order = 1, constant
double k1, k2 # normalizing constants
double basis[ARB] # basis functions
int k, bptr
begin
bptr = 1
do k = 1, order {
if (k == 1)
call amovkd (double(1.0), basis, npts)
else if (k == 2)
call altad (x, basis[bptr], npts, k1, k2)
else {
call amuld (basis[1+npts], basis[bptr-npts], basis[bptr],
npts)
call amulkd (basis[bptr], double(2.0), basis[bptr], npts)
call asubd (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 dcv_bleg (x, npts, order, k1, k2, basis)
double x[npts] # number of data points
int npts # number of points
int order # order of polynomial, 1 is a constant
double k1, k2 # normalizing constants
double basis[ARB] # array of basis functions
int k, bptr
double ri, ri1, ri2
begin
bptr = 1
do k = 1, order {
if (k == 1)
call amovkd (double(1.0), basis, npts)
else if (k == 2)
call altad (x, basis[bptr], npts, k1, k2)
else {
ri = k
ri1 = (double(2.0) * ri - double(3.0)) / (ri - double(1.0))
ri2 = - (ri - double(2.0)) / (ri - double(1.0))
call amuld (basis[1+npts], basis[bptr-npts], basis[bptr],
npts)
call awsud (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 dcv_bspline1 (x, npts, npieces, k1, k2, basis, left)
double x[npts] # set of data points
int npts # number of points
int npieces # number of polynomial pieces minus 1
double k1, k2 # normalizing constants
double basis[ARB] # basis functions
int left[ARB] # indices of the appropriate spline functions
int k
begin
call altad (x, basis[1+npts], npts, k1, k2)
call achtdi (basis[1+npts], left, npts)
call aminki (left, npieces, left, npts)
do k = 1, npts {
basis[npts+k] = max (double(0.0), min (double(1.0),
basis[npts+k] - left[k]))
basis[k] = max (double(0.0), min (double(1.0), double(1.0) -
basis[npts+k]))
}
end
# CV_BSPLINE3 -- Procedure to evaluate all the non-zero basis functions
# for a cubic spline.
procedure dcv_bspline3 (x, npts, npieces, k1, k2, basis, left)
double x[npts] # array of data points
int npts # number of data points
int npieces # number of polynomial pieces minus 1
double k1, k2 # normalizing constants
double basis[ARB] # array of basis functions
int left[ARB] # array of indices for first non-zero spline
int i
pointer sp, sx, tx
double dsx, dtx
begin
# allocate space
call smark (sp)
call salloc (sx, npts, TY_DOUBLE)
call salloc (tx, npts, TY_DOUBLE)
# calculate the index of the first non-zero coeff
call altad (x, Memd[sx], npts, k1, k2)
call achtdi (Memd[sx], left, npts)
call aminki (left, npieces, left, npts)
do i = 1, npts {
Memd[sx+i-1] = max (double(0.0), min (double(1.0),
Memd[sx+i-1] - left[i]))
Memd[tx+i-1] = max (double(0.0), min (double(1.0), double(1.0) -
Memd[sx+i-1]))
}
# calculate the basis function
#call apowk$t (Mem$t[tx], 3, basis, npts)
do i = 1, npts {
dsx = Memd[sx+i-1]
dtx = Memd[tx+i-1]
basis[i] = dtx * dtx * dtx
basis[npts+i] = double(1.0) + dtx * (double(3.0) + dtx *
(double(3.0) - double(3.0) * dtx))
basis[2*npts+i] = double(1.0) + dsx * (double(3.0) + dsx *
(double(3.0) - double(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
|
