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
|
include <math/gsurfit.h>
# IGS_SOLVE1 -- Fit z = f(x, y).
define SFTYPES "|chebyshev|legendre|" # Surface types
procedure igs_solve1 (sf, x, y, z, w, npts)
pointer sf # GSURFIT pointer
real x[npts] # X points
real y[npts] # Y points
real z[npts] # Z points
real w[npts] # Weights
int npts # Number of points
int i, nfunc, ix, iy
pointer sf1, sf2, resids
int strdic()
include "igsfit.com"
begin
# Determine the function type.
nfunc = strdic (function, function, SZ_LINE, SFTYPES)
# Fit the first surface.
ix = min (2, xorder)
iy = min (2, yorder)
call xgsinit (sf1, nfunc, ix, iy, NO, xmin, xmax, ymin, ymax)
call xgsfit (sf1, x, y, z, w, npts, WTS_USER, i)
switch (i) {
case SINGULAR:
call eprintf ("Singular solution\n")
case NO_DEG_FREEDOM:
call error (0, "No degrees of freedom")
}
# Evaluate the first surface and fit the residuals.
call malloc (resids, npts, TY_REAL)
call xgsvector (sf1, x, y, Memr[resids], npts)
call asubr (z, Memr[resids], Memr[resids], npts)
call xgsinit (sf2, nfunc, xorder, yorder, YES, xmin,xmax,ymin,ymax)
call xgsfit (sf2, x, y, Memr[resids], w, npts, WTS_USER, i)
switch (i) {
case SINGULAR:
call eprintf ("Singular solution\n")
case NO_DEG_FREEDOM:
call error (0, "No degrees of freedom")
}
# Add the two surfaces and free memory.
call xgsadd (sf1, sf2, sf)
call xgsfree (sf1)
call xgsfree (sf2)
call mfree (resids, TY_REAL)
end
# IGS_SOLVE2 -- Fit z = x + f(y).
procedure igs_solve2 (sf, x, y, z, w, npts)
pointer sf # GSURFIT pointer
real x[npts] # X points
real y[npts] # Y points
real z[npts] # Z points
real w[npts] # Weights
int npts # Number of points
int i, nfunc
real a
pointer sf1
int strdic()
real xgsgcoeff()
include "igsfit.com"
begin
nfunc = strdic (function, function, SZ_LINE, SFTYPES)
call xgsinit (sf1, nfunc, 1, yorder, NO, xmin, xmax, ymin, ymax)
call asubr (z, x, z, npts)
call xgsfit (sf1, x, y, z, w, npts, WTS_USER, i)
call aaddr (z, x, z, npts)
switch (i) {
case SINGULAR:
call eprintf ("Singular solution\n")
case NO_DEG_FREEDOM:
call error (0, "No degrees of freedom")
}
call xgsfree (sf)
call xgsinit (sf, nfunc, 2, yorder, NO, xmin, xmax, ymin, ymax)
a = xgsgcoeff (sf1, 1, 1)
a = a + (xmin + xmax) / 2
call xgsscoeff (sf, 1, 1, a)
a = (xmax - xmin) / 2
call xgsscoeff (sf, 2, 1, a)
do i = 2, yorder {
a = xgsgcoeff (sf1, 1, i)
call xgsscoeff (sf, 1, i, a)
}
call xgsfree (sf1)
end
# IGS_SOLVE3 -- Fit z = y + f(x).
procedure igs_solve3 (sf, x, y, z, w, npts)
pointer sf # GSURFIT pointer
real x[npts] # X points
real y[npts] # Y points
real z[npts] # Z points
real w[npts] # Weights
int npts # Number of points
int i, nfunc
real a
pointer sf1
int strdic()
real xgsgcoeff()
include "igsfit.com"
begin
nfunc = strdic (function, function, SZ_LINE, SFTYPES)
call xgsinit (sf1, nfunc, xorder, 1, NO, xmin, xmax, ymin, ymax)
call asubr (z, y, z, npts)
call xgsfit (sf1, x, y, z, w, npts, WTS_USER, i)
call aaddr (z, y, z, npts)
switch (i) {
case SINGULAR:
call eprintf ("Singular solution\n")
case NO_DEG_FREEDOM:
call error (0, "No degrees of freedom")
}
call xgsfree (sf)
call xgsinit (sf, nfunc, xorder, 2, NO, xmin, xmax, ymin, ymax)
a = xgsgcoeff (sf1, 1, 1)
a = a + (ymin + ymax) / 2
call xgsscoeff (sf, 1, 1, a)
a = (ymax - ymin) / 2
call xgsscoeff (sf, 1, 2, a)
do i = 2, xorder {
a = xgsgcoeff (sf1, i, 1)
call xgsscoeff (sf, i, 1, a)
}
call xgsfree (sf1)
end
|
