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
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <fset.h>
define POLY_INSTRUMENTAL 1
define POLY_UNIFORM 2
define POLY_STATISTICAL 3
define MAX_ITEMS 1000 # maximum number of data elements
# POLYFIT -- Fit a polynomial to the list of input pairs (x,y, sigmay)
# A polynomial fit of user specifiable order is made to the data.
#
# y = a0 + a1*x + a2*x**2 + ...
#
# The values for the coeficients a0 - aN are printed on the
# first line of the standard output. The uncertainties in the
# coeficients are then printed on the next line.
#
# Optionally (verbose = yes) the values for chi-square, ftest,
# and correlation coefficient are printed along with the calculated
# values from the fit for the dependent variable.
#
# If listdata = yes, then the only output will be pairs of
# X,Yc values. Yc is the value of dependent variable as Calculated
# from the fit. This option allows piping to the GRAPH task.
#
# The data are taken from STDIN, a file, or a list of files.
# In the latter case, each data file results in an independent
# set of results.
# The routines REGRES and MATINV from Bevington are used to perfrom
# the fit.
procedure t_polyfit()
char fname[SZ_FNAME], weights[SZ_FNAME]
bool verbose, listdata
int filelist, order, weighting
int clpopni(), clgfil(), clgeti(), clgwrd()
bool clgetb()
define exit_ 91
begin
# Input can come from the standard input, a file, or a list of files.
# The following procedure makes both cases look like a list of files.
filelist = clpopni ("input")
order = clgeti ("order")
weighting = clgwrd ("weighting", weights, SZ_FNAME,
",instrumental,uniform,statistical,")
verbose = clgetb ("verbose")
listdata = clgetb ("listdata")
while (clgfil (filelist, fname, SZ_FNAME) != EOF)
call pf_fitdatalist (fname, order, weighting, verbose, listdata)
call clpcls (filelist)
end
# PF_FITDATALIST -- Perform polynomial fit to data from named file.
procedure pf_fitdatalist (listfile, order, weighting, verbose, listdata)
char listfile[SZ_FNAME] # input file
int order # polynomial order
int weighting # mode of weighting fit
bool verbose, listdata # printout options
int i, in, item, m[50], mode, nterms, line_number
real x[MAX_ITEMS], y[MAX_ITEMS], yfit[MAX_ITEMS], sigy[MAX_ITEMS]
real a0, a[50], siga0, siga[50], r[50], rmul, chisqr, ftest
real stdev
extern pf_fctn()
bool fp_equalr()
int fscan(), nscan(), open()
errchk open, fscan, printf
begin
# Set term selection array for REGRES.
for (i=1; i <= 50; i=i+1)
m[i] = i
in = open (listfile, READ_ONLY, TEXT_FILE)
# Read successive X,Y, SIGMAY triples from the standard input,
# accumulating the values in the arrays X, Y, weight. Skip list
# elements containing less than two numbers. The maximum number of
# elements that can be read is fixed.
item = 1
line_number = 0
while ((fscan (in) != EOF) && (item < MAX_ITEMS+1)) {
line_number = line_number + 1
call gargr (x[item])
call gargr (y[item])
# There must be two items per entry for the x,y pair.
if (nscan() < 2) {
call eprintf ("Bad entry in list on line:%d - item ignored\n")
call pargi (line_number)
next
}
# Set undefined errors to 0.0.
if (weighting == POLY_INSTRUMENTAL) {
call gargr (sigy[item])
if (nscan() < 3) {
call eprintf ("Undefined sigmay on line:%d - item ignored\n")
call pargi (line_number)
next
} else if (fp_equalr (sigy[item], 0.0)) {
call eprintf ("Zero-valued sigmay on line:%d - item ignored\n")
call pargi (line_number)
next
}
} else
sigy[item] = 0
item = item + 1
}
item = item - 1
if (item > MAX_ITEMS) {
call printf ("Number of data elements exceeded - max=%d\n")
call pargi (MAX_ITEMS)
}
if (item <= order) {
call eprintf ("Not enough data for fit: order=%d, items=%d\n")
call pargi (order)
call pargi (item)
goto exit_
}
# It is necessary to scale the dependent variable values
# to 1.0 on average to minimize the dynamic range during
# matrix inversion.
nterms = order
switch (weighting) {
case POLY_INSTRUMENTAL:
mode = 1
case POLY_UNIFORM:
mode = 0
case POLY_STATISTICAL:
mode = -1
default:
mode = 0
}
call pf_regres (x, y, sigy, item, nterms, m, mode, yfit,
a0, a, siga0, siga, r, rmul, chisqr, ftest, pf_fctn)
# Compute standard deviation of residuals from reduced chi-sqr.
switch (weighting) {
case POLY_STATISTICAL, POLY_INSTRUMENTAL:
stdev = 0.0
do i = 1, item
stdev = stdev + (y[i] - yfit[i]) ** 2
stdev = sqrt (stdev / (item - 1))
case POLY_UNIFORM:
stdev = sqrt ((item - nterms - 1) * chisqr / (item - 1))
default:
stdev = sqrt ((item - nterms - 1) * chisqr / (item - 1))
}
# Print coefficients scaled back to input Y levels
if (!listdata) {
call printf ("%12.7g")
call pargr (a0)
for (i=1; i <= order; i=i+1) {
call printf (" %12.7g")
call pargr (a[i])
}
call printf ("\n")
# Print sigmas.
call printf ("%12.7g")
call pargr (siga0)
for (i=1; i <= order; i=i+1) {
call printf (" %12.7g")
call pargr (siga[i])
}
call printf ("\n")
# If verbose option specified, also print chi-square, ftest
# correlation coefficient, standard deviation of residuals,
# number of input pairs, and calculated y-values.
# ***25Nov85,SeH - ftest is undefined when correlation = 1.
if (verbose) {
if (fp_equalr (ftest, -99999.)) {
call printf ("\nchi sqr: %7g ftest: UNDEF ")
call pargr (chisqr)
if (fp_equalr (rmul, -99999.)) {
call printf (" correlation: UNDEF\n")
call pargr (rmul)
} else {
call printf (" correlation: %7g\n")
call pargr (rmul)
}
} else {
call printf ("\nchi sqr: %7g ftest: %7g ")
call pargr (chisqr)
call pargr (ftest)
if (fp_equalr (rmul, -99999.)) {
call printf (" correlation: UNDEF\n")
call pargr (rmul)
} else {
call printf ("correlation: %7g\n")
call pargr (rmul)
}
}
call printf (" nr pts: %7g std dev res: %0.6g\n")
call pargi (item)
call pargr (stdev)
call printf ("\nx(data) y(calc) y(data) sigy(data)\n")
}
}
if (verbose || listdata) {
for (i=1; i <= item; i=i+1) {
call printf ("%7g %7g %7g %7g\n")
call pargr (x[i])
call pargr (yfit[i])
call pargr (y[i])
call pargr (sigy[i])
}
}
exit_
call close (in)
end
|
