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 "dgsurfitdef.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 dgsrefit (sf, x, y, z, w, ier)
pointer sf # surface descriptor
double x[ARB] # array of x values
double y[ARB] # array of y values
double z[ARB] # data array
double 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
double adotd()
errchk smark, salloc, sfree
begin
# clear accumulator
call aclrd (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_DOUBLE)
switch (GS_TYPE(sf)) {
case GS_LEGENDRE:
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)
call dgs_bleg (x, GS_NPTS(sf), GS_XORDER(sf), GS_XMAXMIN(sf),
GS_XRANGE(sf), XBASIS(GS_XBASIS(sf)))
call dgs_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_DOUBLE)
call malloc (GS_YBASIS(sf), GS_NPTS(sf) * GS_YORDER(sf),
TY_DOUBLE)
call dgs_bcheb (x, GS_NPTS(sf), GS_XORDER(sf), GS_XMAXMIN(sf),
GS_XRANGE(sf), XBASIS(GS_XBASIS(sf)))
call dgs_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_DOUBLE)
call malloc (GS_YBASIS(sf), GS_NPTS(sf) * GS_YORDER(sf),
TY_DOUBLE)
call dgs_bpol (x, GS_NPTS(sf), GS_XORDER(sf), GS_XMAXMIN(sf),
GS_XRANGE(sf), XBASIS(GS_XBASIS(sf)))
call dgs_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_DOUBLE)
# index the pointers
vzptr = GS_VECTOR(sf) - 1
byptr = GS_YBASIS(sf)
switch (GS_TYPE(sf)) {
case GS_LEGENDRE, GS_CHEBYSHEV, GS_POLYNOMIAL:
call amuld (w, z, Memd[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 amuld (Memd[GS_WZ(sf)], YBASIS(byptr), Memd[bwz],
GS_NPTS(sf))
bxptr = GS_XBASIS(sf)
do k = 1, xorder {
vindex = vzptr + k
VECTOR(vindex) = VECTOR(vindex) + adotd (Memd[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 dgschoslv (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
|
