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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <math/gsurfit.h>
include "gsurfitdef.h"
# 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.
procedure gsrefit (sf, x, y, z, w, ier)
pointer sf # surface descriptor
real x[ARB] # array of x values
real y[ARB] # array of y values
real z[ARB] # data array
real 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
real adotr()
errchk smark, salloc, sfree
begin
# clear accumulator
call aclrr (VECTOR(GS_VECTOR(sf)), GS_NCOEFF(sf))
# if first call to gsefit calculate basis functions
if (GS_XBASIS(sf) == NULL || GS_YBASIS(sf) == NULL) {
call malloc (GS_WZ(sf), GS_NPTS(sf), TY_REAL)
switch (GS_TYPE(sf)) {
case GS_LEGENDRE:
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)
call rgs_bleg (x, GS_NPTS(sf), GS_XORDER(sf), GS_XMAXMIN(sf),
GS_XRANGE(sf), XBASIS(GS_XBASIS(sf)))
call rgs_bleg (y, GS_NPTS(sf), GS_YORDER(sf), GS_YMAXMIN(sf),
GS_YRANGE(sf), YBASIS(GS_YBASIS(sf)))
case GS_CHEBYSHEV:
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)
call rgs_bcheb (x, GS_NPTS(sf), GS_XORDER(sf), GS_XMAXMIN(sf),
GS_XRANGE(sf), XBASIS(GS_XBASIS(sf)))
call rgs_bcheb (y, GS_NPTS(sf), GS_YORDER(sf), GS_YMAXMIN(sf),
GS_YRANGE(sf), YBASIS(GS_YBASIS(sf)))
case GS_POLYNOMIAL:
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)
call rgs_bpol (x, GS_NPTS(sf), GS_XORDER(sf), GS_XMAXMIN(sf),
GS_XRANGE(sf), XBASIS(GS_XBASIS(sf)))
call rgs_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)
call salloc (bwz, GS_NPTS(sf), TY_REAL)
# index the pointers
vzptr = GS_VECTOR(sf) - 1
byptr = GS_YBASIS(sf)
switch (GS_TYPE(sf)) {
case GS_LEGENDRE, GS_CHEBYSHEV, GS_POLYNOMIAL:
call amulr (w, z, Memr[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 amulr (Memr[GS_WZ(sf)], YBASIS(byptr), Memr[bwz],
GS_NPTS(sf))
bxptr = GS_XBASIS(sf)
do k = 1, xorder {
vindex = vzptr + k
VECTOR(vindex) = VECTOR(vindex) + adotr (Memr[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 rgschoslv (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
|
