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
|
include "../nlfit/nlfit.h"
# PREP_DATA -- Prepare the x,y,err data to be fitted.
#
# This routine looks for problems in the x,y,err vectors as INDEF-valued
# data or errors, zeroed errors and such, and cleans up to the standards
# needed by the nlfit/curfit packages. This means setting the weight vector
# properly, as well as setting eventual x,y weird values to something safe.
# Routine also builds the weights to allow chi-square computation.
#
# Created 2/21/96 (I Busko).
procedure prep_data (x, y, err, npix, nl, wts)
real x[ARB] # io: x data
real y[ARB] # io: y data
real err[ARB] # i: err data
int npix # i: number of data points
pointer nl # i: NLFIT pointer
pointer wts # o: weight vector.
int i
real xsafe, ysafe
real sigma, epadu, rnoise, rn2
int nl_stati()
real nl_statr()
begin
sigma = nl_statr (nl, "sigma")
epadu = nl_statr (nl, "epadu")
rnoise = nl_statr (nl, "readnoise")
rn2 = rnoise * rnoise
do i = 1, npix {
# If any INDEF value, or zeroed error bar, set weight to zero
# once and for all. Otherwise, error is a valid one, thus build
# weight for chi-square computation.
if ( IS_INDEFR (x[i]) ||
IS_INDEFR (y[i]) ||
IS_INDEFR (err[i]) ||
(err[i] == 0.0) ) {
Memr[wts+i-1] = 0.0
} else {
Memr[wts+i-1] = 1.0 / err[i]
# If sigma parameter was provided by user, use it
# instead of any error bars present in the data.
if (!IS_INDEF(sigma))
Memr[wts+i-1] = 1.0 / sigma
# If ccd type was set, this takes preference.
if (nl_stati(nl,"errtyp") == CCD) {
if (y[i] <= 0.) {
if (rnoise > 0)
Memr[wts+i-1] = 1.0 / rnoise
else
Memr[wts+i-1] = 0.0
} else
Memr[wts+i-1] = 1.0 / sqrt (y[i]/epadu + rn2)
}
}
}
# Set all the zero weighted values to something safe.
do i = 1, npix {
if ( Memr[wts+i-1] != 0.0 ) {
xsafe = x[i]
ysafe = y[i]
}
}
do i = 1, npix {
if ( Memr[wts+i-1] == 0.0) {
x[i] = xsafe
y[i] = ysafe
}
}
end
|
