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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
# SF_BCHEB -- Procedure to evaluate all the non-zero Chebyshev functions for
# a set of points and given order.
procedure sf_bcheb (x, npts, order, k1, k2, basis)
real x[npts] # array of data points
int npts # number of points
int order # order of polynomial, order = 1, constant
real k1, k2 # normalizing constants
real basis[ARB] # basis functions
int k, bptr
begin
bptr = 1
do k = 1, order {
if (k == 1)
call amovkr (1., basis, npts)
else if (k == 2)
call altar (x, basis[bptr], npts, k1, k2)
else {
call amulr (basis[1+npts], basis[bptr-npts], basis[bptr],
npts)
call amulkr (basis[bptr], 2., basis[bptr], npts)
call asubr (basis[bptr], basis[bptr-2*npts], basis[bptr], npts)
}
bptr = bptr + npts
}
end
# SF_BLEG -- Procedure to evaluate all the non zero Legendre function
# for a given order and set of points.
procedure sf_bleg (x, npts, order, k1, k2, basis)
real x[npts] # number of data points
int npts # number of points
int order # order of polynomial, 1 is a constant
real k1, k2 # normalizing constants
real basis[ARB] # array of basis functions
int k, bptr
real ri, ri1, ri2
begin
bptr = 1
do k = 1, order {
if (k == 1)
call amovkr (1., basis, npts)
else if (k == 2)
call altar (x, basis[bptr], npts, k1, k2)
else {
ri = k
ri1 = (2. * ri - 3.) / (ri - 1.)
ri2 = - (ri - 2.) / (ri - 1.)
call amulr (basis[1+npts], basis[bptr-npts], basis[bptr],
npts)
call awsur (basis[bptr], basis[bptr-2*npts],
basis[bptr], npts, ri1, ri2)
}
bptr = bptr + npts
}
end
# SF_BSPLINE1 -- Evaluate all the non-zero spline1 functions for a set
# of points.
procedure sf_bspline1 (x, npts, npieces, k1, k2, basis, left)
real x[npts] # set of data points
int npts # number of points
int npieces # number of polynomial pieces minus 1
real k1, k2 # normalizing constants
real basis[ARB] # basis functions
int left[ARB] # indices of the appropriate spline functions
int k
begin
call altar (x, basis[1+npts], npts, k1, k2)
call achtri (basis[1+npts], left, npts)
call aminki (left, npieces, left, npts)
do k = 1, npts {
basis[npts+k] = basis[npts+k] - left[k]
basis[k] = 1. - basis[npts+k]
}
end
# SF_BSPLINE3 -- Procedure to evaluate all the non-zero basis functions
# for a cubic spline.
procedure sf_bspline3 (x, npts, npieces, k1, k2, basis, left)
real x[npts] # array of data points
int npts # number of data points
int npieces # number of polynomial pieces minus 1
real k1, k2 # normalizing constants
real basis[ARB] # array of basis functions
int left[ARB] # array of indices for first non-zero spline
int i
pointer sp, sx, tx
begin
# allocate space
call smark (sp)
call salloc (sx, npts, TY_REAL)
call salloc (tx, npts, TY_REAL)
# calculate the index of the first non-zero coeff
call altar (x, Memr[sx], npts, k1, k2)
call achtri (Memr[sx], left, npts)
call aminki (left, npieces, left, npts)
# normalize x to 0 to 1
do i = 1, npts {
Memr[sx+i-1] = Memr[sx+i-1] - left[i]
Memr[tx+i-1] = 1. - Memr[sx+i-1]
}
# calculate the basis function
call apowkr (Memr[tx], 3, basis, npts)
do i = 1, npts {
basis[npts+i] = 1. + Memr[tx+i-1] * (3. + Memr[tx+i-1] * (3. -
3. * Memr[tx+i-1]))
basis[2*npts+i] = 1. + Memr[sx+i-1] * (3. + Memr[sx+i-1] * (3. -
3. * Memr[sx+i-1]))
}
call apowkr (Memr[sx], 3, basis[1+3*npts], npts)
# release space
call sfree (sp)
end
|
