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
|
include <math.h>
include "../lib/center.h"
define TOL 0.001 # tolerance for fitting algorithm
# AP_LGCTR1D -- Procedure to compute the center from the 1D marginal
# distributions using a simplified version of the optimal filtering
# technique and addopting a Gaussian model for the fit. The method
# is streamlined by replacing the Gaussian with a simple triangle
# following L.Goad.
int procedure ap_lgctr1d (ctrpix, nx, ny, cx, cy, sigma, maxiter, norm,
skysigma, xc, yc, xerr, yerr)
real ctrpix[nx, ny] # object to be centered
int nx, ny # dimensions of subarray
real cx, cy # center in subraster coordinates
real sigma # sigma of PSF
int maxiter # maximum number of iterations
real norm # the normalization factor
real skysigma # standard deviation of the pixels
real xc, yc # computed centers
real xerr, yerr # first guess at errors
int nxiter, nyiter
pointer sp, xm, ym
real ratio, constant
int aptopt()
real asumr()
begin
# Allocate working space.
call smark (sp)
call salloc (xm, nx, TY_REAL)
call salloc (ym, ny, TY_REAL)
# Compute the marginal distributions.
call ap_mkmarg (ctrpix, Memr[xm], Memr[ym], nx, ny)
xerr = asumr (Memr[xm], nx)
yerr = asumr (Memr[ym], ny)
call adivkr (Memr[xm], real (nx), Memr[xm], nx)
call adivkr (Memr[ym], real (ny), Memr[ym], ny)
# Compute the x center and error.
xc = cx
nxiter = aptopt (Memr[xm], nx, xc, sigma, TOL, maxiter, YES)
if (xerr <= 0.0)
xerr = INDEFR
else {
if (IS_INDEFR(skysigma))
constant = 0.0
else
constant = 4.0 * SQRTOFPI * sigma * skysigma ** 2
ratio = constant / xerr
xerr = sigma ** 2 / (xerr * norm)
xerr = sqrt (max (xerr, ratio * xerr))
if (xerr > real (nx))
xerr = INDEFR
}
# Compute the y center and error.
yc = cy
nyiter = aptopt (Memr[ym], ny, yc, sigma, TOL, maxiter, YES)
if (yerr <= 0.0)
yerr = INDEFR
else {
if (IS_INDEFR(skysigma))
constant = 0.0
else
constant = 4.0 * SQRTOFPI * sigma * skysigma ** 2
ratio = constant / yerr
yerr = sigma ** 2 / (yerr * norm)
yerr = sqrt (max (yerr, ratio * yerr))
if (yerr > real (ny))
yerr = INDEFR
}
# Return appropriate error code.
call sfree (sp)
if (nxiter < 0 || nyiter < 0)
return (AP_CTR_SINGULAR)
else if (nxiter > maxiter || nyiter > maxiter)
return (AP_CTR_NOCONVERGE)
else
return (AP_OK)
end
|
