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
|
include <mach.h>
include "../lib/apphotdef.h"
include "../lib/apphot.h"
include "../lib/noisedef.h"
include "../lib/photdef.h"
include "../lib/phot.h"
define CONVERSION 2.772588
# AP_WMAG -- Procedure to compute the magnitudes inside a set of apertures for
# a single of object.
int procedure ap_wmag (ap, im, wx, wy, positive, skyval, skysig, nsky)
pointer ap # pointer to the apphot structure
pointer im # pointer to the IRAF image
real wx, wy # object coordinates
int positive # emission or absorption features
real skyval # sky value
real skysig # sky sigma
int nsky # number of sky pixels
int c1, c2, l1, l2, ier, nap
pointer sp, nse, phot, temp
real datamin, datamax, zmag, wsigsq, wvarsky
int apmagbuf()
begin
# Initalize.
phot = AP_PPHOT(ap)
nse = AP_NOISE(ap)
AP_PXCUR(phot) = wx
AP_PYCUR(phot) = wy
if (IS_INDEFR(wx) || IS_INDEFR(wy)) {
AP_OPXCUR(phot) = INDEFR
AP_OPYCUR(phot) = INDEFR
} else {
switch (AP_WCSOUT(ap)) {
case WCS_WORLD, WCS_PHYSICAL:
call ap_ltoo (ap, wx, wy, AP_OPXCUR(phot), AP_OPYCUR(phot), 1)
case WCS_TV:
call ap_ltov (im, wx, wy, AP_OPXCUR(phot), AP_OPYCUR(phot), 1)
default:
AP_OPXCUR(phot) = wx
AP_OPYCUR(phot) = wy
}
}
call amovkd (0.0d0, Memd[AP_SUMS(phot)], AP_NAPERTS(phot))
call amovkd (0.0d0, Memd[AP_AREA(phot)], AP_NAPERTS(phot))
call amovkr (INDEFR, Memr[AP_MAGS(phot)], AP_NAPERTS(phot))
call amovkr (INDEFR, Memr[AP_MAGERRS(phot)], AP_NAPERTS(phot))
# Make sure the center is defined.
if (IS_INDEFR(wx) || IS_INDEFR(wy))
return (AP_APERT_NOAPERT)
# Fetch the aperture pixels.
ier = apmagbuf (ap, im, wx, wy, c1, c2, l1, l2)
if (ier == AP_APERT_NOAPERT)
return (AP_APERT_NOAPERT)
call smark (sp)
call salloc (temp, AP_NAPERTS(phot), TY_REAL)
# Compute the min and max.
if (IS_INDEFR(AP_DATAMIN(ap)))
datamin = -MAX_REAL
else
datamin = AP_DATAMIN(ap)
if (IS_INDEFR(AP_DATAMAX(ap)))
datamax = MAX_REAL
else
datamax = AP_DATAMAX(ap)
# Do photometry for all the apertures.
AP_NMINAP(phot) = AP_NMAXAP(phot) + 1
call amulkr (Memr[AP_APERTS(phot)], AP_SCALE(ap), Memr[temp],
AP_NAPERTS(phot))
switch (AP_PWEIGHTS(phot)) {
case AP_PWCONSTANT:
if (IS_INDEFR(AP_DATAMIN(ap)) && IS_INDEFR(AP_DATAMAX(ap)))
call apmeasure (im, wx, wy, c1, c2, l1, l2, Memr[temp],
Memd[AP_SUMS(phot)], Memd[AP_AREA(phot)], AP_NMAXAP(phot))
else
call apbmeasure (im, wx, wy, c1, c2, l1, l2, datamin,
datamax, Memr[temp], Memd[AP_SUMS(phot)],
Memd[AP_AREA(phot)], AP_NMAXAP(phot), AP_NMINAP(phot))
case AP_PWCONE:
wsigsq = AP_FWHMPSF(ap) * AP_SCALE(ap)
wvarsky = (AP_READNOISE(nse) / AP_EPADU(nse)) ** 2
if (IS_INDEFR(AP_DATAMIN(ap)) && IS_INDEFR(AP_DATAMAX(ap)))
call ap_tmeasure (im, wx, wy, c1, c2, l1, l2, Memr[temp],
Memd[AP_SUMS(phot)], Memd[AP_AREA(phot)], AP_NMAXAP(phot),
wsigsq, AP_EPADU(nse), wvarsky)
else
call ap_btmeasure (im, wx, wy, c1, c2, l1, l2, datamin,
datamax, Memr[temp], Memd[AP_SUMS(phot)],
Memd[AP_AREA(phot)], AP_NMAXAP(phot), AP_NMINAP(phot),
wsigsq, AP_EPADU(nse), wvarsky)
case AP_PWGAUSS:
wsigsq = (AP_FWHMPSF(ap) * AP_SCALE(ap)) ** 2 / CONVERSION
wvarsky = (AP_READNOISE(nse) / AP_EPADU(nse)) ** 2
if (IS_INDEFR(AP_DATAMIN(ap)) && IS_INDEFR(AP_DATAMAX(ap)))
call ap_gmeasure (im, wx, wy, c1, c2, l1, l2, Memr[temp],
Memd[AP_SUMS(phot)], Memd[AP_AREA(phot)], AP_NMAXAP(phot),
wsigsq, AP_EPADU(nse), wvarsky)
else
call ap_bgmeasure (im, wx, wy, c1, c2, l1, l2,
datamin, datamax, Memr[temp], Memd[AP_SUMS(phot)],
Memd[AP_AREA(phot)], AP_NMAXAP(phot), AP_NMINAP(phot),
wsigsq, AP_EPADU(nse), wvarsky)
default:
if (IS_INDEFR(AP_DATAMIN(ap)) && IS_INDEFR(AP_DATAMAX(ap)))
call apmeasure (im, wx, wy, c1, c2, l1, l2, Memr[temp],
Memd[AP_SUMS(phot)], Memd[AP_AREA(phot)], AP_NMAXAP(phot))
else
call apbmeasure (im, wx, wy, c1, c2, l1, l2, datamin,
datamax, Memr[temp], Memd[AP_SUMS(phot)],
Memd[AP_AREA(phot)], AP_NMAXAP(phot), AP_NMINAP(phot))
}
# Make sure that the sky value has been defined.
if (IS_INDEFR(skyval))
ier = AP_APERT_NOSKYMODE
else {
# Check for bad pixels.
if ((ier == AP_OK) && (AP_NMINAP(phot) <= AP_NMAXAP(phot)))
ier = AP_APERT_BADDATA
nap = min (AP_NMINAP(phot) - 1, AP_NMAXAP(phot))
# Compute the magnitudes and errors.
if (positive == YES)
call apcopmags (Memd[AP_SUMS(phot)], Memd[AP_AREA(phot)],
Memr[AP_MAGS(phot)], Memr[AP_MAGERRS(phot)], nap, skyval,
skysig, nsky, AP_ZMAG(phot), AP_NOISEFUNCTION(nse),
AP_EPADU(nse))
else
call apconmags (Memd[AP_SUMS(phot)], Memd[AP_AREA(phot)],
Memr[AP_MAGS(phot)], Memr[AP_MAGERRS(phot)], nap, skyval,
skysig, nsky, AP_ZMAG(phot), AP_NOISEFUNCTION(nse),
AP_EPADU(nse), AP_READNOISE(nse))
# Compute correction for itime.
zmag = 2.5 * log10 (AP_ITIME(ap))
call aaddkr (Memr[AP_MAGS(phot)], zmag, Memr[AP_MAGS(phot)], nap)
}
call sfree (sp)
if (ier != AP_OK)
return (ier)
else
return (AP_OK)
end
|
