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
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
|
include <mach.h>
include <imhdr.h>
include "../lib/daophotdef.h"
include "../lib/apseldef.h"
include "../lib/allstardef.h"
# DP_ASTAR -- Begin doing the photometry.
procedure dp_astar (dao, im, subim, allfd, rejfd, cache, savesub, version)
pointer dao # pointer to the daophot structure
pointer im # pointer to the input image
pointer subim # pointer to the subtracted image
int allfd # file descriptor for the output photometry file
int rejfd # file descriptor for rejections files
int cache # cache the data ?
int savesub # save the subtracted image ?
int version # version number
bool clip
int nstar, row1_number, row2_number, niter, itsky, x1, x2, y1, y2, istar
int lstar, ldummy, newsky
pointer apsel, psffit, allstar, sp, indices, colpoints
real fradius, sepcrt, sepmin, clmpmx, wcrit, radius, xmin, xmax, ymin, ymax
real csharp, rdummy
int dp_alphot(), dp_alnstar()
pointer dp_gst(), dp_gwt(), dp_gdc()
real dp_usepsf()
begin
# Define some daophot structure pointers.
apsel = DP_APSEL(dao)
psffit = DP_PSFFIT (dao)
allstar = DP_ALLSTAR (dao)
# Store the original fitting radius
fradius = DP_FITRAD(dao)
DP_FITRAD(dao) = min (DP_FITRAD(dao), DP_PSFRAD(dao))
DP_SFITRAD(dao) = DP_FITRAD(dao) * DP_SCALE(dao)
# Get some working memory.
call smark (sp)
call salloc (indices, NAPPAR, TY_INT)
call salloc (colpoints, ALL_NOUTCOLUMN, TY_INT)
# Define some daophot constants.
nstar = DP_APNUM(apsel)
csharp = 1.4427 * Memr[DP_PSFPARS(psffit)] *
Memr[DP_PSFPARS(psffit)+1] / dp_usepsf (DP_PSFUNCTION(psffit), 0.0,
0.0, DP_PSFHEIGHT(psffit), Memr[DP_PSFPARS(psffit)],
Memr[DP_PSFLUT(psffit)], DP_PSFSIZE(psffit), DP_NVLTABLE(psffit),
DP_NFEXTABLE(psffit), 0.0, 0.0, rdummy, rdummy)
if (IS_INDEFR(DP_MERGERAD(dao)))
sepcrt = 2.0 * (Memr[DP_PSFPARS(psffit)] ** 2 +
Memr[DP_PSFPARS(psffit)+1] ** 2)
else
sepcrt = DP_MERGERAD(dao) ** 2
sepmin = min (1.0, FRACTION_MINSEP * sepcrt)
itsky = 3
# Initialize the output photometry file for writing.
row1_number = 0
row2_number = 0
if (DP_TEXT(dao) == YES) {
call dp_xnewal (dao, allfd)
if (rejfd != NULL)
call dp_xnewal (dao, rejfd)
} else {
call dp_tnewal (dao, allfd, Memi[colpoints])
if (rejfd != NULL)
call dp_tnewal (dao, rejfd, Memi[colpoints])
}
# Allocate memory for the input, output and fitting arrays.
call dp_gnindices (Memi[indices])
call dp_rmemapsel (dao, Memi[indices], NAPPAR, nstar)
call dp_almemstar (dao, nstar, DP_MAXGROUP(dao))
call dp_cache (dao, im, subim, cache)
# Read in the input data and assign the initial weights.
call dp_setwt (dao, im)
# Convert the magnitudes to relative brightnesses.
call dp_alzero (dao, Memr[DP_APMAG(apsel)], nstar)
# Initialize all the fit/nofit indices and the error codes.
call amovki (NO, Memi[DP_ASKIP(allstar)], nstar)
call amovki (ALLERR_OK, Memi[DP_AIER(allstar)], nstar)
# Remove stars that have INDEF centers, are off the image altogether,
# or are too close to another star before beginning to fit.
call dp_strip (Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
Memr[DP_APMSKY(apsel)], Memi[DP_ASKIP(allstar)],
Memi[DP_AIER(allstar)], nstar, sepmin, int(IM_LEN(im,1)),
int(IM_LEN(im,2)), DP_FITRAD(dao), DP_VERBOSE(dao))
# Write out results for the rejected stars.
call dp_wout (dao, im, Memr[DP_APXCEN(apsel)+nstar],
Memr[DP_APYCEN(apsel)+nstar], Memr[DP_APXCEN(apsel)+nstar],
Memr[DP_APYCEN(apsel)+nstar], DP_APNUM(apsel) - nstar)
if (DP_TEXT(dao) == YES) {
call dp_xalwrite (allfd, rejfd, Memi[DP_APID(apsel)],
Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)],
Memr[DP_APMSKY(apsel)], Memr[DP_APCHI(apsel)],
Memr[DP_ANUMER(allstar)], Memr[DP_ADENOM(allstar)],
Memi[DP_ASKIP(allstar)], Memi[DP_AIER(allstar)], niter,
nstar + 1, DP_APNUM(apsel), DP_PSFMAG(psffit), csharp)
} else {
call dp_talwrite (allfd, rejfd, Memi[colpoints],
Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
Memr[DP_APERR(apsel)], Memr[DP_APMSKY(apsel)],
Memr[DP_APCHI(apsel)], Memr[DP_ANUMER(allstar)],
Memr[DP_ADENOM(allstar)], Memi[DP_ASKIP(allstar)],
Memi[DP_AIER(allstar)], niter, nstar + 1, DP_APNUM(apsel),
row1_number, row2_number, DP_PSFMAG(psffit), csharp)
}
# Do some initialization for the fit.
clmpmx = CLAMP_FRACTION * DP_FITRAD(dao)
call aclrr (Memr[DP_AXOLD(allstar)], nstar)
call aclrr (Memr[DP_AYOLD(allstar)], nstar)
call amovkr (clmpmx, Memr[DP_AXCLAMP(allstar)], nstar)
call amovkr (clmpmx, Memr[DP_AYCLAMP(allstar)], nstar)
call amovkr (1.0, Memr[DP_ASUMWT(allstar)], nstar)
# Begin iterating.
for (niter = 1; (nstar > 0) && (niter <= DP_MAXITER (dao));
niter = niter + 1) {
if (DP_VERBOSE(dao) == YES) {
call printf ("NITER = %d\n")
call pargi (niter)
}
if ((DP_CLIPEXP(dao) != 0) && (niter >= 4))
clip = true
else
clip = false
# Define the critical errors.
if (niter >= 15)
wcrit = WCRIT15
else if (niter >= 10)
wcrit = WCRIT10
else if (niter >= 5)
wcrit = WCRIT5
else
wcrit = WCRIT0
# Sort the data on y.
call dp_alsort (Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
Memr[DP_APMSKY(apsel)], Memr[DP_ASUMWT(allstar)],
Memr[DP_AXOLD(allstar)], Memr[DP_AYOLD(allstar)],
Memr[DP_AXCLAMP(allstar)], Memr[DP_AYCLAMP(allstar)], nstar)
# Compute the working radius. This is either the fitting radius
# or the outer sky radius if this is an iteration during which
# sky is to be determined whichever is larger.
if ((DP_FITSKY(dao) == YES) && (mod (niter, itsky) == 0)) {
newsky = YES
if ((DP_ANNULUS(dao) + DP_DANNULUS(dao)) > DP_FITRAD(dao))
radius = DP_ANNULUS(dao) + DP_DANNULUS(dao)
else
radius = DP_FITRAD(dao)
} else {
newsky = NO
radius = DP_FITRAD(dao)
}
# Define region of the image used to fit the remaining stars.
call alimr (Memr[DP_APXCEN(apsel)], nstar, xmin, xmax)
call alimr (Memr[DP_APYCEN(apsel)], nstar, ymin, ymax)
x1 = max (1, min (IM_LEN(im,1), int (xmin-radius)+1))
x2 = max (1, min (IM_LEN(im,1), int (xmax+radius)))
y1 = max (1, min (IM_LEN(im,2), int (ymin-radius)+1))
y2 = max (1, min (IM_LEN (im,2), int (ymax+radius)))
# Reinitialize the weight and scratch arrays / images .
call dp_wstinit (dao, im, Memr[DP_APXCEN(apsel)],
Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
nstar, radius, x1, x2, y1, y2)
# Recompute the initial sky estimates if that switch is enabled.
if (newsky == YES)
call dp_alsky (dao, im, Memr[DP_APXCEN(apsel)],
Memr[DP_APYCEN(apsel)], Memr[DP_APMSKY(apsel)], nstar,
x1, x2, y1, y2, DP_ANNULUS(dao), DP_ANNULUS(dao) +
DP_DANNULUS(dao), 100, -MAX_REAL)
# Group the remaining stars.
call dp_regroup (Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
Memr[DP_APMSKY(apsel)], Memr[DP_ASUMWT(allstar)],
Memr[DP_AXOLD(allstar)], Memr[DP_AYOLD(allstar)],
Memr[DP_AXCLAMP(allstar)], Memr[DP_AYCLAMP(allstar)], nstar,
DP_FITRAD(dao), Memi[DP_ALAST(allstar)])
# Reset the error codes.
call amovki (ALLERR_OK, Memi[DP_AIER(allstar)], nstar)
# Do the serious fitting one group at a time.
for (istar = 1; istar <= nstar; istar = lstar + 1) {
# Do the fitting a group at a time.
lstar = dp_alphot (dao, im, nstar, istar, niter, sepcrt,
sepmin, wcrit, clip, clmpmx, version)
# Write the results.
if (DP_TEXT(dao) == YES) {
call dp_xalwrite (allfd, rejfd, Memi[DP_APID(apsel)],
Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)],
Memr[DP_APMSKY(apsel)], Memr[DP_APCHI(apsel)],
Memr[DP_ANUMER(allstar)], Memr[DP_ADENOM(allstar)],
Memi[DP_ASKIP(allstar)], Memi[DP_AIER(allstar)],
niter, istar, lstar, DP_PSFMAG(psffit), csharp)
} else {
call dp_talwrite (allfd, rejfd, Memi[colpoints],
Memi[DP_APID(apsel)], Memr[DP_APXCEN(apsel)],
Memr[DP_APYCEN(apsel)], Memr[DP_APMAG(apsel)],
Memr[DP_APERR(apsel)], Memr[DP_APMSKY(apsel)],
Memr[DP_APCHI(apsel)], Memr[DP_ANUMER(allstar)],
Memr[DP_ADENOM(allstar)], Memi[DP_ASKIP(allstar)],
Memi[DP_AIER(allstar)], niter, istar, lstar, row1_number,
row2_number, DP_PSFMAG(psffit), csharp)
}
}
# Find the last star in the list that still needs more work.
nstar = dp_alnstar (Memi[DP_ASKIP(allstar)], nstar)
# Remove the fitted stars from the list.
for (istar = 1; (istar < nstar) && nstar > 0; istar = istar + 1) {
call dp_alswitch (Memi[DP_APID(apsel)],
Memr[DP_APXCEN(apsel)], Memr[DP_APYCEN(apsel)],
Memr[DP_APMAG(apsel)], Memr[DP_APERR(apsel)],
Memr[DP_APMSKY(apsel)], Memr[DP_ASUMWT(allstar)],
Memr[DP_AXOLD(allstar)], Memr[DP_AYOLD(allstar)],
Memr[DP_AXCLAMP(allstar)], Memr[DP_AYCLAMP(allstar)],
Memi[DP_ASKIP(allstar)], istar, nstar)
}
# Flush the output buffers.
if (dp_gwt (dao, im, ldummy, ldummy, READ_WRITE, YES) == NULL)
;
if (dp_gst (dao, im, ldummy, ldummy, READ_ONLY, YES) == NULL)
;
if (dp_gdc (dao, im, ldummy, ldummy, READ_WRITE, YES) == NULL)
;
}
# Release cached memory.
call dp_uncache (dao, subim, savesub)
call sfree (sp)
# Restore the original fitting radius
DP_FITRAD(dao) = fradius
DP_SFITRAD(dao) = DP_FITRAD(dao) * DP_SCALE(dao)
end
# DP_ALNSTAR -- Compute the number of stars left
int procedure dp_alnstar (skip, nstar)
int skip[ARB] # skip array
int nstar # number of stars
begin
while (nstar > 0) {
if (skip[nstar] == NO)
break
nstar = nstar - 1
}
return (nstar)
end
# DP_ALSWITCH -- Switch array elements.
procedure dp_alswitch (id, xcen, ycen, mag, magerr, sky, chi, xold, yold,
xclamp, yclamp, skip, istar, nstar)
int id[ARB] # array of ids
real xcen[ARB] # array of x centers
real ycen[ARB] # array of y centers
real mag[ARB] # array of magnitudes
real magerr[ARB] # array of magnitude errors
real sky[ARB] # array of sky values
real chi[ARB] # array of chi values
real xold[ARB] # old x array
real yold[ARB] # old y array
real xclamp[ARB] # array of x clamps
real yclamp[ARB] # array of y clamps
int skip[ARB] # skip array
int istar # next star
int nstar # total number of stars
int dp_alnstar()
begin
if (skip[istar] == YES) {
id[istar] = id[nstar]
xcen[istar] = xcen[nstar]
ycen[istar] = ycen[nstar]
mag[istar] = mag[nstar]
magerr[istar] = magerr[nstar]
sky[istar] = sky[nstar]
chi[istar] = chi[nstar]
xold[istar] = xold[nstar]
yold[istar] = yold[nstar]
xclamp[istar] = xclamp[nstar]
yclamp[istar] = yclamp[nstar]
skip[istar] = NO
nstar = nstar - 1
nstar = dp_alnstar (skip, nstar)
}
end
|
