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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
# AWVG -- Compute the mean and standard deviation (sigma) of a sample. Pixels
# whose value lies outside the specified lower and upper limits are not used.
# If the upper and lower limits have the same value (e.g., zero), no limit
# checking is performed. The number of pixels in the sample is returned as the
# function value.
int procedure awvg$t (a, npix, mean, sigma, lcut, hcut)
PIXEL a[ARB]
$if (datatype == dl)
double mean, sigma, lcut, hcut
$else
real mean, sigma, lcut, hcut
$endif
double sum, sumsq, value, temp
int npix, i, ngpix
begin
sum = 0.0
sumsq = 0.0
ngpix = 0
# Accumulate sum, sum of squares. The test to disable limit checking
# requires numerical equality of two floating point numbers; this should
# be ok since they are used as flags not as numbers (they are not used
# in computations).
if (hcut == lcut) {
do i = 1, npix {
$if (datatype == x)
value = abs (a[i])
$else
value = a[i]
$endif
sum = sum + value
sumsq = sumsq + value ** 2
}
ngpix = npix
} else {
do i = 1, npix {
$if (datatype == x)
value = abs (a[i])
$else
value = a[i]
$endif
if (value >= lcut && value <= hcut) {
ngpix = ngpix + 1
sum = sum + value
sumsq = sumsq + value ** 2
}
}
}
switch (ngpix) { # compute mean and sigma
case 0:
$if (datatype == dl)
mean = $INDEFD
sigma = $INDEFD
$else
mean = $INDEFR
sigma = $INDEFR
$endif
case 1:
mean = sum
$if (datatype == dl)
sigma = $INDEFD
$else
sigma = $INDEFR
$endif
default:
mean = sum / ngpix
temp = (sumsq - (sum/ngpix) * sum) / (ngpix - 1)
if (temp < 0) # possible with roundoff error
sigma = 0.0
else
sigma = sqrt (temp)
}
return (ngpix)
end
|
