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
|
include <math/gsurfit.h>
include "epix.h"
# EP_GSFIT -- Fit the background annulus.
procedure ep_gsfit (ep, data, mask, x, y, w, nx, ny, gs)
pointer ep # EPIX structure
real data[nx,ny] # Data subraster
int mask[nx,ny] # Mask subraster
real x[nx,ny] # X positions
real y[nx,ny] # Y positions
real w[nx,ny] # Weights
int nx, ny # Subraster size
pointer gs # Surface pointer (returned)
int i, j, n, npts, xo, yo
pointer sp, work
real amedr()
begin
call smark (sp)
call salloc (work, nx * ny, TY_REAL)
gs = NULL
npts = nx * ny
if (EP_XORDER(ep) == 0 || EP_YORDER(ep) == 0) {
n = 0
do j = 1, ny {
do i = 1, nx {
if (mask[i,j] == 2) {
Memr[work+n] = data[i,j]
n = n + 1
}
}
}
call amovkr (amedr (Memr[work], n), Memr[work], npts)
xo = 1
yo = 1
} else {
call amovr (data, Memr[work], npts)
xo = EP_XORDER(ep)
yo = EP_YORDER(ep)
}
n = 0
do j = 1, ny {
do i = 1, nx {
x[i,j] = i
y[i,j] = j
if (mask[i,j] == 2) {
w[i,j] = 1.
n = n + 1
} else
w[i,j] = 0.
}
}
if (n > 7) {
repeat {
call gsinit (gs, GS_POLYNOMIAL, xo, yo, YES,
1., real (nx), 1., real (ny))
call gsfit (gs, x, y, Memr[work], w, npts, WTS_USER, n)
if (n == OK)
break
xo = max (1, xo - 1)
yo = max (1, yo - 1)
}
} else
call eprintf ("ERROR: Insufficient background points\n")
call sfree (sp)
end
|
