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
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
|
include <mach.h>
include <math/gsurfit.h>
define ECFTYPES "|chebyshev|legendre|" # Fit types
# ECF_SOLVE -- Fit
#
# f(x, slope*y+offset) = (y+slope*offset)*z
#
# with offset minimizing the RMS.
procedure ecf_solve (ecf, x, y, z, w, r, npts, fixedorder)
pointer ecf # GSURFIT pointer
double x[npts] # X points
double y[npts] # Y points
double z[npts] # Z points
double w[npts] # Weights
double r[npts] # Residuals
int npts # Number of points
int fixedorder # Fixed order?
int i, j, k, err
double ya, yb, newrms, ecf_rms()
pointer sp, y1, ecf1
errchk ecf_solve1
include "ecffit.com"
define fit_ 99
begin
if (fixedorder == YES) {
call ecf_solve1 (ecf, x, y, z, w, r, npts)
return
}
call smark (sp)
call salloc (y1, npts, TY_DOUBLE)
# Determine if the orders are reversed.
j = 1
k = 1
do i = 1, npts {
if (z[i] < z[j])
j = i
if (z[i] > z[k])
k = i
}
if (y[j] >= y[k]) {
slope = 1
offset = max (offset, int(1. - ymin))
} else {
slope = -1
offset = max (offset, int(1. + ymax))
}
call dgsfree (ecf)
shift = 0.
rms = MAX_DOUBLE
j = 1
k = 0
for (i=offset;;i=i+j) {
if (slope == 1) {
ya = i + ymin
yb = i + ymax
} else {
ya = i - ymax
yb = i - ymin
}
if (ya < 1.)
break
call altmd (y, Memd[y1], npts, double(slope), double(i))
call amuld (Memd[y1], z, r, npts)
fit_ call dgsinit (ecf1, gstype, xorder, yorder, YES, xmin, xmax, ya, yb)
call dgsfit (ecf1, x, Memd[y1], r, w, npts, WTS_USER, err)
if (err != OK) {
if (xorder > 2 || yorder > 2) {
call dgsfree (ecf)
xorder = max (2, xorder - 1)
yorder = max (2, yorder - 1)
goto fit_
}
switch (err) {
case SINGULAR:
call dgsfree (ecf)
ecf = ecf1
call eprintf ("Singular solution\n")
case NO_DEG_FREEDOM:
call sfree (sp)
call error (0, "No degrees of freedom")
}
}
call dgsvector (ecf1, x, Memd[y1], r, npts)
call adivd (r, Memd[y1], r, npts)
call asubd (z, r, r, npts)
newrms = ecf_rms (r, w, npts)
k = k + 1
if (newrms / rms < 0.999) {
call dgsfree (ecf)
ecf = ecf1
offset = i
rms = newrms
} else {
call dgsfree (ecf1)
if (k > 2)
break
i = offset
j = -j
}
}
call altmd (y, Memd[y1], npts, double(slope), double(offset))
call dgsvector (ecf, x, Memd[y1], r, npts)
call adivd (r, Memd[y1], r, npts)
call asubd (z, r, r, npts)
call sfree (sp)
end
# ECF_SOLVE1 -- Fit f(x, y+offset) = (y+offset)*z with offset fixed.
procedure ecf_solve1 (ecf, x, y, z, w, r, npts)
pointer ecf # GSURFIT pointer
double x[npts] # X points
double y[npts] # Y points
double z[npts] # Z points
double w[npts] # Weights
double r[npts] # Residuals
int npts # Number of points
int err
pointer sp, y1
double ya, yb, ecf_rms()
include "ecffit.com"
define fit_ 99
begin
call smark (sp)
call salloc (y1, npts, TY_DOUBLE)
call dgsfree (ecf)
shift = 0.
if (slope == 1) {
offset = max (offset, int(1. - ymin))
ya = offset + ymin
yb = offset + ymax
} else {
offset = max (offset, int(1. + ymax))
ya = offset - ymax
yb = offset - ymin
}
call altmd (y, Memd[y1], npts, double (slope), double (offset))
call amuld (Memd[y1], z, r, npts)
fit_ call dgsinit (ecf, gstype, xorder, yorder, YES, xmin, xmax,
min (ya, yb), max (ya, yb))
call dgsfit (ecf, x, Memd[y1], r, w, npts, WTS_USER, err)
if (err != OK) {
if (xorder > 2 || yorder > 2) {
call dgsfree (ecf)
xorder = max (2, xorder - 1)
yorder = max (2, yorder - 1)
goto fit_
}
switch (err) {
case SINGULAR:
call eprintf ("Singular solution\n")
case NO_DEG_FREEDOM:
call sfree (sp)
call error (0, "No degrees of freedom")
}
}
call dgsvector (ecf, x, Memd[y1], r, npts)
call adivd (r, Memd[y1], r, npts)
call asubd (z, r, r, npts)
rms = ecf_rms (r, w, npts)
call sfree (sp)
end
|
