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
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <math/gsurfit.h>
$if (datatype == r)
include "gsurfitdef.h"
$else
include "dgsurfitdef.h"
$endif
# GSREFIT -- Procedure to refit the surface assuming that the x, y and w
# values and the matrices MATRIX and CHOFAC have remained unchanged. It
# is necessary only to accumulate a new VECTOR. The new coefficients
# are calculated by forward and back subsitution and stored in COEFF.
$if (datatype == r)
procedure gsrefit (sf, x, y, z, w, ier)
$else
procedure dgsrefit (sf, x, y, z, w, ier)
$endif
pointer sf # surface descriptor
PIXEL x[ARB] # array of x values
PIXEL y[ARB] # array of y values
PIXEL z[ARB] # data array
PIXEL w[ARB] # array of weights
int ier # ier = OK, everything OK
# ier = SINGULAR, matrix is singular, 1 or more
# coefficients are 0.
# ier = NO_DEG_FREEDOM, too few points to solve matrix
int k, l
int xorder, nfree, maxorder
pointer sp, vzptr, vindex, bxptr, byptr, bwz
PIXEL adot$t()
errchk smark, salloc, sfree
begin
# clear accumulator
call aclr$t (VECTOR(GS_VECTOR(sf)), GS_NCOEFF(sf))
# if first call to gsefit calculate basis functions
if (GS_XBASIS(sf) == NULL || GS_YBASIS(sf) == NULL) {
$if (datatype == r)
call malloc (GS_WZ(sf), GS_NPTS(sf), TY_REAL)
$else
call malloc (GS_WZ(sf), GS_NPTS(sf), TY_DOUBLE)
$endif
switch (GS_TYPE(sf)) {
case GS_LEGENDRE:
$if (datatype == r)
call malloc (GS_XBASIS(sf), GS_NPTS(sf) * GS_XORDER(sf),
TY_REAL)
call malloc (GS_YBASIS(sf), GS_NPTS(sf) * GS_YORDER(sf),
TY_REAL)
$else
call malloc (GS_XBASIS(sf), GS_NPTS(sf) * GS_XORDER(sf),
TY_DOUBLE)
call malloc (GS_YBASIS(sf), GS_NPTS(sf) * GS_YORDER(sf),
TY_DOUBLE)
$endif
call $tgs_bleg (x, GS_NPTS(sf), GS_XORDER(sf), GS_XMAXMIN(sf),
GS_XRANGE(sf), XBASIS(GS_XBASIS(sf)))
call $tgs_bleg (y, GS_NPTS(sf), GS_YORDER(sf), GS_YMAXMIN(sf),
GS_YRANGE(sf), YBASIS(GS_YBASIS(sf)))
case GS_CHEBYSHEV:
$if (datatype == r)
call malloc (GS_XBASIS(sf), GS_NPTS(sf) * GS_XORDER(sf),
TY_REAL)
call malloc (GS_YBASIS(sf), GS_NPTS(sf) * GS_YORDER(sf),
TY_REAL)
$else
call malloc (GS_XBASIS(sf), GS_NPTS(sf) * GS_XORDER(sf),
TY_DOUBLE)
call malloc (GS_YBASIS(sf), GS_NPTS(sf) * GS_YORDER(sf),
TY_DOUBLE)
$endif
call $tgs_bcheb (x, GS_NPTS(sf), GS_XORDER(sf), GS_XMAXMIN(sf),
GS_XRANGE(sf), XBASIS(GS_XBASIS(sf)))
call $tgs_bcheb (y, GS_NPTS(sf), GS_YORDER(sf), GS_YMAXMIN(sf),
GS_YRANGE(sf), YBASIS(GS_YBASIS(sf)))
case GS_POLYNOMIAL:
$if (datatype == r)
call malloc (GS_XBASIS(sf), GS_NPTS(sf) * GS_XORDER(sf),
TY_REAL)
call malloc (GS_YBASIS(sf), GS_NPTS(sf) * GS_YORDER(sf),
TY_REAL)
$else
call malloc (GS_XBASIS(sf), GS_NPTS(sf) * GS_XORDER(sf),
TY_DOUBLE)
call malloc (GS_YBASIS(sf), GS_NPTS(sf) * GS_YORDER(sf),
TY_DOUBLE)
$endif
call $tgs_bpol (x, GS_NPTS(sf), GS_XORDER(sf), GS_XMAXMIN(sf),
GS_XRANGE(sf), XBASIS(GS_XBASIS(sf)))
call $tgs_bpol (y, GS_NPTS(sf), GS_YORDER(sf), GS_YMAXMIN(sf),
GS_YRANGE(sf), YBASIS(GS_YBASIS(sf)))
default:
call error (0, "GSREFIT: Unknown curve type.")
}
}
call smark (sp)
$if (datatype == r)
call salloc (bwz, GS_NPTS(sf), TY_REAL)
$else
call salloc (bwz, GS_NPTS(sf), TY_DOUBLE)
$endif
# index the pointers
vzptr = GS_VECTOR(sf) - 1
byptr = GS_YBASIS(sf)
switch (GS_TYPE(sf)) {
case GS_LEGENDRE, GS_CHEBYSHEV, GS_POLYNOMIAL:
call amul$t (w, z, Mem$t[GS_WZ(sf)], GS_NPTS(sf))
xorder = GS_XORDER(sf)
maxorder = max (GS_XORDER(sf) + 1, GS_YORDER(sf) + 1)
do l = 1, GS_YORDER(sf) {
call amul$t (Mem$t[GS_WZ(sf)], YBASIS(byptr), Mem$t[bwz],
GS_NPTS(sf))
bxptr = GS_XBASIS(sf)
do k = 1, xorder {
vindex = vzptr + k
VECTOR(vindex) = VECTOR(vindex) + adot$t (Mem$t[bwz],
XBASIS(bxptr), GS_NPTS(sf))
bxptr = bxptr + GS_NPTS(sf)
}
vzptr = vzptr + xorder
switch (GS_XTERMS(sf)) {
case GS_XNONE:
xorder = 1
case GS_XHALF:
if ((l + GS_XORDER(sf) + 1) > maxorder)
xorder = xorder - 1
default:
;
}
byptr = byptr + GS_NPTS(sf)
}
default:
call error (0, "GSACCUM: Unknown curve type.")
}
# test for number of degrees of freedom
ier = OK
nfree = GS_NPTS(sf) - GS_NCOEFF(sf)
if (nfree < 0) {
ier = NO_DEG_FREEDOM
return
}
# calculate the values of the coefficients
switch (GS_TYPE(sf)) {
case GS_LEGENDRE, GS_CHEBYSHEV, GS_POLYNOMIAL:
# solve for the coefficients by forward and back substitution
call $tgschoslv (CHOFAC(GS_CHOFAC(sf)), GS_NCOEFF(sf),
GS_NCOEFF(sf), VECTOR(GS_VECTOR(sf)), COEFF(GS_COEFF(sf)))
default:
call error (0, "GSSOLVE: Illegal surface type.")
}
# release the space
call sfree (sp)
end
|
