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
245
246
247
248
249
250
251
252
253
254
255
256
|
include <gset.h>
include <pkg/gtools.h>
include "apertures.h"
# AP_NOISE -- Model residuals.
procedure ap_noise (ap, gain, dbuf, nc, nl, c1, l1, sbuf, spec, profile, nx, ny,
xs, ys, sum2, sum4, nsum, nbin, dmin, dmax)
pointer ap # Aperture structure
real gain # Gain
pointer dbuf # Data buffer
int nc, nl # Size of data buffer
int c1, l1 # Origin of data buffer
pointer sbuf # Sky buffer (NULL if no sky)
real spec[ny] # Normalization spectrum
real profile[ny,nx] # Profile
int nx, ny # Size of profile array
int xs[ny], ys # Start of spectrum in image
real sum2[nbin] # Sum of residuals squared in bin
real sum4[nbin] # Sum of residuals squared in bin
int nsum[nbin] # Number of values in bin
int nbin # Number of bins
real dmin, dmax # Data limits of bins
int i, ix, iy, ix1, ix2
real dstep, low, high, s, x1, x2, model, data, ap_cveval()
pointer cv, sptr, dptr
begin
dstep = (dmax - dmin) / nbin
i = AP_AXIS(ap)
low = AP_CEN(ap,i) + AP_LOW(ap,i)
high = AP_CEN(ap,i) + AP_HIGH(ap,i)
cv = AP_CV(ap)
do iy = 1, ny {
i = iy + ys - 1
s = ap_cveval (cv, real (i))
x1 = max (0.5, low + s)
x2 = min (c1 + nc - 0.49, high + s)
if (x1 > x2)
next
ix1 = nint (x1) - xs[iy] + 1
ix2 = nint (x2) - xs[iy] + 1
s = spec[iy]
if (sbuf != NULL)
sptr = sbuf + (iy - 1) * nx - 1
dptr = dbuf + (i - l1) * nc + nint(x1) - c1
do ix = ix1, ix2 {
if (sbuf != NULL) {
model = (s * profile[iy,ix] + Memr[sptr]) / gain
sptr = sptr + 1
} else
model = (s * profile[iy,ix]) / gain
data = Memr[dptr] / gain
dptr = dptr + 1
if (model < dmin || model >= dmax)
next
i = (model - dmin) / dstep + 1
sum2[i] = sum2[i] + (data - model) ** 2
sum4[i] = sum4[i] + (data - model) ** 4
nsum[i] = nsum[i] + 1
}
}
end
define HELP "noao$twodspec/apextract/apnoise.key"
define PROMPT "apextract options"
# AP_NPLOT -- Plot and examine noise characteristics.
procedure ap_nplot (image, im, sigma, sigerr, npts, dmin, dmax)
char image[SZ_FNAME] # Image
pointer im # Image pointer
real sigma[npts] # Sigma values
real sigerr[npts] # Sigma errors
int npts # Number of sigma values
real dmin, dmax # Data min and max
real rdnoise # Read noise
real gain # Gain
int i, newgraph, wcs, key
real wx, wy, x, x1, x2, dx, y, ymin, ymax
pointer sp, cmd, gp, gt
int gt_gcur()
real apgimr()
#int apgwrd()
#bool ap_answer()
pointer gt_init()
errchk ap_gopen
begin
# Query user.
call smark (sp)
call salloc (cmd, SZ_LINE, TY_CHAR)
#call sprintf (Memc[cmd], SZ_LINE, "Edit apertures for %s?")
# call pargstr (image)
#if (!ap_answer ("ansedit", Memc[cmd])) {
# call sfree (sp)
# return
#}
gain = apgimr ("gain", im)
rdnoise = apgimr ("readnoise", im)
dx = (dmax - dmin) / npts
x1 = dmin + dx / 2
x2 = dmax - dx / 2
ymin = sigma[1] - sigerr[1]
ymax = sigma[1] + sigerr[1]
do i = 2, npts {
ymin = min (ymin, sigma[i] - sigerr[i])
ymax = max (ymax, sigma[i] + sigerr[i])
}
# Set up the graphics.
call sprintf (Memc[cmd], SZ_LINE, "Noise characteristics of image %s")
call pargstr (image)
call ap_gopen (gp)
gt = gt_init()
call gt_sets (gt, GTTITLE, Memc[cmd])
call gt_sets (gt, GTXLABEL, "Data value")
call gt_sets (gt, GTYLABEL, "Sigma")
call gt_sets (gt, GTTYPE, "mark")
call gt_sets (gt, GTMARK, "plus")
# Enter cursor loop.
key = 'r'
repeat {
switch (key) {
case '?': # Print help text.
call gpagefile (gp, HELP, PROMPT)
case ':': # Colon commands.
if (Memc[cmd] == '/')
call gt_colon (Memc[cmd], gp, gt, newgraph)
else
call ap_ncolon (Memc[cmd], rdnoise, gain, newgraph)
case 'q':
break
case 'r': # Redraw the graph.
newgraph = YES
case 'w': # Window graph
call gt_window (gt, gp, "gcur", newgraph)
case 'I': # Interrupt
call fatal (0, "Interrupt")
default: # Ring bell for unrecognized commands.
call printf ("\007")
}
# Update the graph if needed.
if (newgraph == YES) {
call sprintf (Memc[cmd], SZ_LINE,
"Read noise = %g e-, Gain = %g e-/DN")
call pargr (rdnoise)
call pargr (gain)
call gt_sets (gt, GTPARAMS, Memc[cmd])
call gclear (gp)
y = sqrt ((rdnoise/gain)**2 + dmax/gain)
call gswind (gp, dmin, dmax, ymin, max (ymax, y))
call gt_swind (gp, gt)
call gt_labax (gp, gt)
do i = 1, npts {
if (sigma[i] > 0) {
x = x1 + (i-1) * dx
call gmark (gp, x, sigma[i], GM_VEBAR+GM_HLINE, -dx/2,
-sigerr[i])
}
}
do i = 1, npts {
x = x1 + (i-1) * dx
y = sqrt ((rdnoise/gain)**2 + x/gain)
if (i == 1)
call gamove (gp, x, y)
else
call gadraw (gp, x, y)
}
newgraph = NO
}
} until (gt_gcur ("gcur", wx, wy, wcs, key, Memc[cmd], SZ_LINE) == EOF)
# Free memory.
call gt_free (gt)
call sfree (sp)
end
# List of colon commands.
define CMDS "|readnoise|gain|"
define RDNOISE 1 # Read noise
define GAIN 2 # Gain
# AP_NCOLON -- Respond to colon command from ap_nplot.
procedure ap_ncolon (command, rdnoise, gain, newgraph)
char command[ARB] # Colon command
real rdnoise # Readout noise
real gain # Gain
int newgraph # New graph?
real rval
int ncmd, nscan(), strdic()
pointer sp, cmd
begin
# Scan the command string and get the first word.
call smark (sp)
call salloc (cmd, SZ_LINE, TY_CHAR)
call sscan (command)
call gargwrd (cmd, SZ_LINE)
ncmd = strdic (cmd, cmd, SZ_LINE, CMDS)
switch (ncmd) {
case RDNOISE:
call gargr (rval)
if (nscan() == 2) {
rdnoise = rval
newgraph = YES
} else {
call printf ("rdnoise %g\n")
call pargr (rdnoise)
}
case GAIN:
call gargr (rval)
if (nscan() == 2) {
gain = rval
newgraph = YES
} else {
call printf ("gain %g\n")
call pargr (gain)
}
default:
call printf ("Unrecognized or ambiguous command\007")
}
call sfree (sp)
end
|
