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
|
include <mach.h>
include "../lib/apphotdef.h"
include "../lib/apphot.h"
include "../lib/fitpsfdef.h"
include "../lib/noisedef.h"
include "../lib/fitpsf.h"
# APSFFIT -- Procedure to fit an analytic function to the PSF.
int procedure apsffit (ap, im, wx, wy)
pointer ap # pointer to the apphot structure
pointer im # pointer to the IRAF image
real wx, wy # object coordinates
int ier, fier
pointer psf, nse
real datamin, datamax, dmin, dmax, threshold
int apfbuf(), apsfradgauss(), apsfelgauss(), apsfmoments()
begin
# Initialize.
psf = AP_PPSF(ap)
nse = AP_NOISE(ap)
AP_PFXCUR(psf) = wx
AP_PFYCUR(psf) = wy
if (IS_INDEFR(wx) || IS_INDEFR(wy)) {
AP_OPFXCUR(psf) = INDEFR
AP_OPFYCUR(psf) = INDEFR
} else {
switch (AP_WCSOUT(ap)) {
case WCS_WORLD, WCS_PHYSICAL:
call ap_ltoo (ap, wx, wy, AP_OPFXCUR(psf), AP_OPFYCUR(psf), 1)
case WCS_TV:
call ap_ltov (im, wx, wy, AP_OPFXCUR(psf), AP_OPFYCUR(psf), 1)
default:
AP_OPFXCUR(psf) = wx
AP_OPFYCUR(psf) = wy
}
}
call amovkr (INDEFR, Memr[AP_PPARS(psf)], AP_MAXNPARS(psf))
call amovkr (INDEFR, Memr[AP_PPERRS(psf)], AP_MAXNPARS(psf))
# Fetch the buffer of pixels.
ier = apfbuf (ap, im, wx, wy)
if (ier == AP_NOPSFAREA)
return (AP_NOPSFAREA)
# Compute the min and max of the data subraster.
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)
switch (AP_PSFUNCTION(psf)) {
case AP_RADGAUSS:
fier = apsfradgauss (Memr[AP_PSFPIX(psf)], AP_PNX(psf), AP_PNY(psf),
AP_POSITIVE(ap), AP_FWHMPSF(ap) * AP_SCALE(ap), datamin,
datamax, AP_NOISEFUNCTION(nse), AP_EPADU(nse),
AP_READNOISE(nse) / AP_EPADU(nse), AP_PMAXITER(psf),
AP_PK2(psf), AP_PNREJECT(psf), Memr[AP_PPARS(psf)],
Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf))
Memr[AP_PPARS(psf)+1] = Memr[AP_PPARS(psf)+1] + wx - AP_PXC(psf)
Memr[AP_PPARS(psf)+2] = Memr[AP_PPARS(psf)+2] + wy - AP_PYC(psf)
case AP_ELLGAUSS:
fier = apsfelgauss (Memr[AP_PSFPIX(psf)], AP_PNX(psf), AP_PNY(psf),
AP_POSITIVE(ap), AP_FWHMPSF(ap) * AP_SCALE(ap), datamin,
datamax, AP_NOISEFUNCTION(nse), AP_EPADU(nse),
AP_READNOISE(nse) / AP_EPADU(nse), AP_PMAXITER(psf),
AP_PK2(psf), AP_PNREJECT(psf), Memr[AP_PPARS(psf)],
Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf))
Memr[AP_PPARS(psf)+1] = Memr[AP_PPARS(psf)+1] + wx - AP_PXC(psf)
Memr[AP_PPARS(psf)+2] = Memr[AP_PPARS(psf)+2] + wy - AP_PYC(psf)
case AP_MOMENTS:
call alimr (Memr[AP_PSFPIX(psf)], AP_PNX(psf) * AP_PNY(psf),
dmin, dmax)
dmin = max (dmin, datamin)
dmax = min (dmax, datamax)
threshold = 0.0
if (AP_POSITIVE(ap) == YES)
fier = apsfmoments (Memr[AP_PSFPIX(psf)], AP_PNX(psf),
AP_PNY(psf), dmin + threshold, dmax,
AP_POSITIVE(ap), Memr[AP_PPARS(psf)], Memr[AP_PPERRS(psf)],
AP_PSFNPARS(psf))
else
fier = apsfmoments (Memr[AP_PSFPIX(psf)], AP_PNX(psf),
AP_PNY(psf), dmax - threshold, dmin,
AP_POSITIVE(ap), Memr[AP_PPARS(psf)],
Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf))
Memr[AP_PPARS(psf)+1] = Memr[AP_PPARS(psf)+1] + wx - AP_PXC(psf)
Memr[AP_PPARS(psf)+2] = Memr[AP_PPARS(psf)+2] + wy - AP_PYC(psf)
default:
# do nothing gracefully
}
switch (AP_WCSOUT(ap)) {
case WCS_WORLD, WCS_PHYSICAL:
call ap_ltoo (ap, Memr[AP_PPARS(psf)+1], Memr[AP_PPARS(psf)+2],
Memr[AP_PPARS(psf)+1], Memr[AP_PPARS(psf)+2], 1)
case WCS_TV:
call ap_ltov (im, Memr[AP_PPARS(psf)+1], Memr[AP_PPARS(psf)+2],
Memr[AP_PPARS(psf)+1], Memr[AP_PPARS(psf)+2], 1)
default:
;
}
# Return the appropriate error code.
if (fier == AP_OK) {
if (ier == AP_PSF_OUTOFBOUNDS)
return (AP_PSF_OUTOFBOUNDS)
else
return (AP_OK)
} else if (fier == AP_NPSF_TOO_SMALL) {
call amovkr (INDEFR, Memr[AP_PPARS(psf)], AP_PSFNPARS(psf))
call amovkr (INDEFR, Memr[AP_PPERRS(psf)], AP_PSFNPARS(psf))
return (AP_NPSF_TOO_SMALL)
} else
return (fier)
end
|
