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
|
include "epix.h"
# EP_STATISTICS -- Compute and print statistics for the input aperture.
procedure ep_statistics (ep, ap, xa, ya, xb, yb, box)
pointer ep # EPIX structure
int ap # Aperture type
int xa, ya, xb, yb # Aperture coordinates
int box # Print box?
int i, x1, x2, y1, y2
pointer mask, x, y, w, gs
begin
i = max (5., abs (EP_SEARCH(ep))+EP_BUFFER(ep)+EP_WIDTH(ep)+1)
x1 = min (xa, xb) - i
x2 = max (xa, xb) + i
y1 = min (ya, xb) - i
y2 = max (ya, yb) + i
EP_OUTDATA(ep) = NULL
call ep_gindata (ep, x1, x2, y1, y2)
if (EP_INDATA(ep) != NULL) {
call malloc (mask, EP_NPTS(ep), TY_INT)
call malloc (x, EP_NPTS(ep), TY_REAL)
call malloc (y, EP_NPTS(ep), TY_REAL)
call malloc (w, EP_NPTS(ep), TY_REAL)
call ep_search (ep, Memr[EP_INDATA(ep)], EP_NX(ep),
EP_NY(ep), ap, xa, ya, xb, yb)
call ep_mask (ep, mask, ap, xa, ya, xb, yb)
call ep_gsfit (ep, Memr[EP_INDATA(ep)], Memi[mask],
Memr[x], Memr[y], Memr[w], EP_NX(ep), EP_NY(ep), gs)
call ep_statistics1 (Memr[EP_INDATA(ep)], Memi[mask],
EP_NX(ep), EP_NY(ep), EP_X1(ep), EP_Y1(ep),
(xa+xb)/2, (ya+yb)/2, gs)
if (box == YES)
call ep_box (Memr[EP_INDATA(ep)], EP_NX(ep), EP_NY(ep),
EP_X1(ep), EP_Y1(ep), xa, ya, xb, yb)
call mfree (mask, TY_INT)
call mfree (x, TY_REAL)
call mfree (y, TY_REAL)
call mfree (w, TY_REAL)
call gsfree (gs)
}
end
# EP_STATISTICS1 -- Compute and print statistics.
procedure ep_statistics1 (data, mask, nx, ny, x1, y1, x, y, gs)
real data[nx,ny] # Input data subraster
int mask[nx,ny] # Mask subraster
int nx, ny # Size of subraster
int x1, y1 # Origin of subraster
int x, y # Center of object
pointer gs # GSURFIT pointer
int i, j, area, nsky
real flux, sky, sigma, d, gseval()
begin
flux = 0.
area = 0
sky = 0.
sigma = 0.
nsky = 0
do j = 1, ny {
do i = 1, nx {
if (mask[i,j] == 1) {
d = data[i,j]
if (gs != NULL)
d = d - gseval (gs, real (i), real (j))
flux = flux + d
area = area + 1
} else if (mask[i,j] == 2) {
d = data[i,j] - gseval (gs, real (i), real (j))
sky = sky + data[i,j]
sigma = sigma + d * d
nsky = nsky + 1
}
}
}
call printf ("x=%d y=%d z=%d mean=%g area=%d")
call pargi (x)
call pargi (y)
call pargr (data[x-x1+1,y-y1+1])
call pargr (flux / area)
call pargi (area)
if (nsky > 0) {
call printf (" sky=%g sigma=%g nsky=%d")
call pargr (sky / nsky)
call pargr (sqrt (sigma / nsky))
call pargi (nsky)
}
call printf ("\n")
end
# EP_BOX -- Print box of pixel values.
procedure ep_box (data, nx, ny, xo, yo, xa, ya, xb, yb)
real data[nx,ny] # Input data subraster
int nx, ny # Size of subraster
int xo, yo # Origin of subraster
int xa, ya, xb, yb # Aperture
int i, j, x1, x2, y1, y2, x, y
begin
x1 = min (xa, xb)
x2 = max (xa, xb)
y1 = min (ya, yb)
y2 = max (ya, yb)
if (x2 - x1 + 1 <= 10) {
x1 = max (xo, x1 - 1)
x2 = min (xo + nx - 1, x2 + 1)
}
y1 = max (yo, y1 - 1)
y2 = min (yo + ny - 1, y2 + 1)
call printf ("%4w")
do x = x1, x2 {
call printf (" %4d ")
call pargi (x)
}
call printf ("\n")
do y = y2, y1, -1 {
call printf ("%4d")
call pargi (y)
j = y - yo + 1
do x = x1, x2 {
i = x - xo + 1
call printf (" %5g")
call pargr (data[i,j])
}
call printf ("\n")
}
end
|
