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
|
# FINDTHRESH - estimate the expected random error per pixel (in ADU) of
# the background, given the gain and read noise (in electrons) of a CCD.
#
# random error in 1 pixel = sqrt (sky*p(N) + r(N)**2) / p(N)
#
# r(N) is the effective read noise (electrons), corrected for N frames
# p(N) is the effective gain (electrons/ADU), corrected for N frames
#
# In our implementation, the `mean' used to estimate the sky may actually
# be any of `mean', `midpt', or `mode' as in the IMSTATISTICS task.
procedure findthresh (data)
real data {prompt="Sky level (ADU)"}
string images = "" {prompt="List of images"}
string section = "[*,*]" {prompt="Selected image section"}
string center = "mean" {prompt="Central statistical measure",
enum="mean|midpt|mode"}
real binwidth = 0.1 {prompt="Bin width of histogram in sigma\n"}
real gain {prompt="CCD gain in electrons/ADU"}
real readnoise {prompt="CCD read noise in electrons"}
int nframes = 1 {prompt="Number of coadded frames",
min=1}
string coaddtype = "average" {prompt="Type of coaddition",
enum="average|sum"}
bool verbose = yes {prompt="Verbose output?\n"}
string *list1
string *list2
begin
string img, tmpfile, statsfile
real reff, peff, mean, stddev, random
peff = gain
reff = readnoise
if (nframes > 1) {
reff *= sqrt (nframes)
if (coaddtype == "average")
peff *= nframes
if (verbose) {
print ("effective gain = ", peff, " (electrons/ADU)")
print ("effective readnoise = ", reff, " (electrons)\n")
}
}
if (images != "" && $nargs == 0) {
statsfile = mktemp ("tmp$junk")
tmpfile = mktemp ("tmp$junk")
sections (images, > tmpfile)
list1 = tmpfile
while (fscan (list1, img) != EOF) {
imstatistics (img//section, fields=center//",stddev",
lower=INDEF, upper=INDEF, binwidth=binwidth, format-,
> statsfile)
list2 = statsfile
if (fscan (list2, mean, stddev) != 2)
break
list2 = ""; delete (statsfile, ver-, >& "dev$null")
random = sqrt (mean*peff + reff**2) / peff
# round to three decimal places
stddev = real (nint (stddev * 1000.)) / 1000.
random = real (nint (random * 1000.)) / 1000.
if (verbose) {
print (" sigma (computed) = ", random, " (ADU)")
print (" (measured) = ", stddev, " (ADU)\n")
} else
print (random, "\t", stddev)
}
list1 = ""; delete (tmpfile, ver-, >& "dev$null")
list2 = ""; delete (statsfile, ver-, >& "dev$null")
} else {
mean = data
random = sqrt (mean*peff + reff**2) / peff
# round to three decimal places
random = real (nint (random * 1000.)) / 1000.
if (verbose)
print (" sigma (computed) = ", random, " (ADU)")
else
print (random)
}
end
|
