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
|
# RAVERAGE -- Running average, standard deviation, and enevelope of a list.
procedure raverage (input, nwin)
file input {prompt="Input file"}
int nwin {prompt="Number of points in window", min=1}
bool sort = no {prompt="Sort input?"}
real nsig = 0 {prompt="Number of sigma for envelope"}
struct *fd1, *fd2
begin
file in, tmp
int i, j, n, nin
real x, y, z, x1, y1, sum1, sum2, sum3, xavg, yavg, sig
real ylow, yhigh
# Get query parameters.
in = input
nin = nwin
if (nwin < 1)
error (1, "Window size must be greater than zero")
# Buffer the standard input and sort if requested.
tmp = in
if (in == "STDIN") {
tmp = mktemp ("tmp$iraf")
i = 0
while (YES) {
j = scan(x,y)
if (j != 2) {
if (j == EOF)
break
else if (j < 1)
next
else if (j < 2) {
y = x
x = i + 1
}
}
i += 1
print (x, y, >> tmp)
}
if (sort) {
rename (tmp, tmp//"a")
sort (tmp//"a", col=1, num+, rev-, > tmp)
delete (tmp//"a", verify-)
}
} else if (sort) {
tmp = mktemp ("tmp$iraf")
sort (in,, col=1, num+, rev-, > tmp)
}
# Initialize.
i = 0; n = 0; sum1 = 0; sum2 = 0; sum3 = 0
# Accumulate the first window.
fd1 = tmp
while (n<nin) {
j = fscan (fd1, x, y)
if (j != 2) {
if (j == EOF)
break
else if (j < 1)
next
else if (j < 2) {
y = x
x = i + 1
}
}
i += 1
n += 1
sum1 += x
sum2 += y
sum3 += y*y
}
# Read the rest of the list adding and subtracting.
fd2 = tmp
while (YES) {
j = fscan (fd1, x, y)
if (j != 2) {
if (j == EOF)
break
else if (j == 0)
next
else if (j == 1) {
y = x
x = i + 1
}
}
while (YES) {
j = fscan(fd2, x1, y1)
if (j != 2) {
if (j == 0)
next
else if (j == 1) {
y1 = x1
x1 = i - n + 1
}
}
break
}
xavg = sum1 / n
yavg = sum2 / n
sig = sqrt (max (0., (n * sum3 - sum2 * sum2) / (n * max(1,n-1))))
if (nsig > 0) {
ylow = yavg - nsig * sig
yhigh = yavg + nsig * sig
printf ("%g %g %g %d %g %g\n", xavg, yavg, sig, n, ylow, yhigh)
} else
printf ("%g %g %g %d\n", xavg, yavg, sig, n)
i += 1
sum1 += x - x1
sum2 += y - y1
sum3 += y*y - y1*y1
}
fd1 = ""; fd2 = ""
# Compute the final values.
if (n <= 0) {
xavg = INDEF
yavg = INDEF
sig = INDEF
ylow = INDEF
yhigh = INDEF
} else {
xavg = sum1 / n
yavg = sum2 / n
sig = sqrt (max (0., (n * sum3 - sum2 * sum2) / (n * max(1,n-1))))
ylow = yavg - nsig * sig
yhigh = yavg + nsig * sig
}
if (nsig > 0)
printf ("%g %g %g %d %g %g\n", xavg, yavg, sig, n, ylow, yhigh)
else
printf ("%g %g %g %d\n", xavg, yavg, sig, n)
# Delete temporary files.
if (tmp != in)
delete (tmp, verify-)
end
|
