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
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
|
include <imhdr.h>
include <tbset.h>
include "../lib/daophotdef.h"
include "../lib/apseldef.h"
include "../lib/peakdef.h"
# DP_PEAKPHOT -- Fit the PSF to a single star.
procedure dp_peakphot (dao, im, tp, tpout, tprej, ap_text)
pointer dao # pointer to the daophot structure
pointer im # the input image descriptor
int tp # the input photometry file descriptor
int tpout # the ouput photometry file descriptor
int tprej # the rejections file descriptor
bool ap_text # which style of photometry
real rel_bright, xold, yold, x, y, dx, dy, mag, sky, errmag
real chi, sharp, radius, itx, ity, otx, oty
pointer psffit, key, sp, subim, colpoint, indices, fields, perror
int id, in_nrow, instar, lowx, lowy, nxpix, nypix, niter, out_nrow
int rout_nrow, nterm, ier, plen
int tbpsta(), dp_rrphot(), dp_pkfit(), dp_gpkerr()
pointer dp_gsubrast()
begin
# Get the daophot pointers.
psffit = DP_PSFFIT (dao)
# Store the original fitting radius.
radius = DP_FITRAD(dao)
# Check that the fitting radius is less than the psf radius.
DP_FITRAD(dao) = min (DP_FITRAD(dao), DP_PSFRAD(dao))
DP_SFITRAD(dao) = DP_FITRAD(dao) * DP_SCALE(dao)
# Allocate working space.
call smark (sp)
call salloc (indices, NAPPAR, TY_INT)
call salloc (fields, SZ_LINE, TY_CHAR)
call salloc (colpoint, PK_NOUTCOL, TY_INT)
call salloc (perror, SZ_FNAME, TY_CHAR)
# Initialze the output table.
if (DP_TEXT(dao) == YES) {
call dp_xnewpeak (dao, tpout)
if (tprej != NULL)
call dp_xnewpeak (dao, tprej)
} else {
call dp_tnewpeak (dao, tpout, Memi[colpoint])
if (tprej != NULL)
call dp_tnewpeak (dao, tprej, Memi[colpoint])
}
# Intialize the input table.
if (ap_text) {
call pt_kyinit (key)
Memi[indices] = DP_PAPID
Memi[indices+1] = DP_PAPXCEN
Memi[indices+2] = DP_PAPYCEN
Memi[indices+3] = DP_PAPMAG1
Memi[indices+4] = DP_PAPSKY
call dp_gappsf (Memi[indices], Memc[fields], NAPRESULT)
in_nrow = 0
} else {
call dp_tpkinit (tp, Memi[indices])
in_nrow = tbpsta (tp, TBL_NROWS)
}
# Initialize the photometry file reading code.
instar = 0
# Initialize the fitting code.
if (DP_RECENTER(dao) == YES)
nterm = 3
else
nterm = 1
if (DP_FITSKY(dao) == YES)
nterm = nterm + 1
call dp_mempk (dao, nterm)
out_nrow = 0
rout_nrow = 0
repeat {
# Read in the photometry for a single star.
if (dp_rrphot (tp, key, Memc[fields], Memi[indices], id, itx,
ity, sky, mag, instar, in_nrow) == EOF)
break
# Convert to and from logical coordinates.
call dp_win (dao, im, itx, ity, x, y, 1)
call dp_wout (dao, im, x, y, otx, oty, 1)
if (DP_VERBOSE(dao) == YES) {
call printf (
"Star: %5d X: %8.2f Y: %8.2f Mag: %8.2f Sky: %8.2f\n")
call pargi (id)
call pargr (otx)
call pargr (oty)
call pargr (mag)
call pargr (sky)
}
# Check that the center is defined.
if (IS_INDEFR(x) || IS_INDEFR(y)) {
if (DP_VERBOSE(dao) == YES) {
call printf (
"\tWarning: X and/or Y for star %d are undefined\n")
call pargi (id)
}
next
}
# Read in the subraster.
subim = dp_gsubrast (im, x, y, DP_FITRAD(dao), lowx, lowy,
nxpix, nypix)
if (subim == NULL) {
if (DP_VERBOSE(dao) == YES) {
call printf (
"\tWarning: Cannot read in image data for star %d\n")
call pargi (id)
}
next
}
# Save the old x and y values for use with the variable psf
# option.
xold = x
yold = y
call dp_wpsf (dao, im, xold, yold, xold, yold, 1)
# Compute the relative centers and the relative brightness and
# fit the star.
if (IS_INDEFR(sky)) {
ier = PKERR_INDEFSKY
} else {
x = x - lowx + 1.0
y = y - lowy + 1.0
dx = (xold - 1.0) / DP_PSFX(psffit) - 1.0
dy = (yold - 1.0) / DP_PSFY(psffit) - 1.0
if (IS_INDEFR(mag))
mag = DP_PSFMAG (psffit) + DELTA_MAG
rel_bright = DAO_RELBRIGHT (psffit, mag)
ier = dp_pkfit (dao, Memr[subim], nxpix, nypix, DP_FITRAD(dao),
x, y, dx, dy, rel_bright, sky, errmag, chi, sharp, niter)
x = x + lowx - 1.0
y = y + lowy - 1.0
}
call dp_wout (dao, im, x, y, otx, oty, 1)
if (ier != PKERR_OK) {
# Set fitting parameters to INDEF.
mag = INDEFR
niter = 0
errmag = INDEFR
chi = INDEFR
sharp = INDEFR
if (DP_VERBOSE(dao) == YES) {
switch (ier) {
case PKERR_INDEFSKY:
call printf (
"\tWarning: The sky value for star %d is undefined\n")
call pargi (id)
case PKERR_NOPIX:
call printf (
"\tWarning: Too few pixels to fit star %d\n")
call pargi (id)
case PKERR_SINGULAR:
call printf (
"\tWarning: Singular matrix computed for star %d\n")
call pargi (id)
case PKERR_FAINT:
call printf ("\tWarning: Star %d too faint\n")
call pargi (id)
default:
call printf (
"\tWarning: Unknown error detected for star %d\n")
call pargi (id)
}
}
} else {
# Compute the results.
mag = DP_PSFMAG (psffit) - 2.5 * log10 (rel_bright)
errmag = 1.085736 * errmag / rel_bright
if (errmag >= 2.0)
errmag = INDEFR
if (DP_VERBOSE (dao) == YES) {
call printf (
"\tFIT: Star: %5d X: %8.2f Y: %8.2f Mag: %8.2f Sky =%8.2f\n")
call pargi (id)
call pargr (otx)
call pargr (oty)
call pargr (mag)
call pargr (sky)
}
}
# Get the error code.
plen = dp_gpkerr (ier, Memc[perror], SZ_FNAME)
# Now write the results to the output photometry or rejections
# file.
if (DP_TEXT(dao) == YES) {
if ((tprej != NULL) && (ier != PKERR_OK))
call dp_xpkwrite (tprej, id, otx, oty, mag, errmag, sky,
niter, chi, sharp, ier, Memc[perror], plen)
else
call dp_xpkwrite (tpout, id, otx, oty, mag, errmag, sky,
niter, chi, sharp, ier, Memc[perror], plen)
} else {
if ((tprej != NULL) && (ier != PKERR_OK)) {
rout_nrow = rout_nrow + 1
call dp_tpkwrite (tprej, Memi[colpoint], id, otx, oty, mag,
errmag, sky, niter, chi, sharp, ier,
Memc[perror], plen, rout_nrow)
} else {
out_nrow = out_nrow + 1
call dp_tpkwrite (tpout, Memi[colpoint], id, otx, oty, mag,
errmag, sky, niter, chi, sharp, ier, Memc[perror],
plen, out_nrow)
}
}
}
# Free the text file descriptor.
if (ap_text)
call pt_kyfree (key)
# Restore the original fitting radius.
DP_FITRAD(dao) = radius
DP_SFITRAD(dao) = DP_FITRAD(dao) * DP_SCALE(dao)
call sfree (sp)
end
# DP_GPKERR -- Set the PEAK task error code string.
int procedure dp_gpkerr (ier, perror, maxch)
int ier # the integer error code
char perror # the output error code string
int maxch # the maximum size of the error code string
int plen
int gstrcpy()
begin
switch (ier) {
case PKERR_OK:
plen = gstrcpy ("No_error", perror, maxch)
case PKERR_INDEFSKY:
plen = gstrcpy ("Bad_sky", perror, maxch)
case PKERR_NOPIX:
plen = gstrcpy ("Npix_too_few", perror, maxch)
case PKERR_SINGULAR:
plen = gstrcpy ("Singular", perror, maxch)
case PKERR_FAINT:
plen = gstrcpy ("Too_faint", perror, maxch)
default:
plen = gstrcpy ("No_error", perror, maxch)
}
return (plen)
end
|
