aboutsummaryrefslogtreecommitdiff
path: root/noao/nproto/ace/edgewts.xNEW
blob: fdd8d8ddc69ca0134e575ab30256e4ec3d6e62dd (plain) (blame)
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
task	test

procedure test ()

double	dx, dy, r[11], w[11], clgetd()
int	i, nr

begin
	dx = clgetd ("dx")
	dy = clgetd ("dy")
	nr = 11

	call edgewts (dx, dy, r, w, nr)
	do i = 1, nr {
	    call eprintf ("%.2f %.4g\n")
		call pargd (r[i])
		call pargd (w[i])
	}
end

procedure edgewts (dx, dy, r, w, nr)

double	dx			#I Distance from aperture center to pixel center
double	dy			#I Distance from aperture center to pixel center
double	r[nr]			#O Aperture radii
double	w[nr]			#O Weights
int	nr			#O Number of aperture radius points

int	i, j, k, n
double	r2, rmin, rmax, dr, a, d, rap2, y2

begin
	rmin = sqrt ((max(0.,dx-0.6))**2+(max(0.,dy-0.6))**2)
	rmax = sqrt ((dx+0.6)**2+(dy+0.6)**2)
	dr = (rmax - rmin) / nr
	rmin = rmin + dr / 2

	n = 100
	d = 1.0D0 / (2 * n + 1)
	a = d * d

	do k = 1, nr {
	    rap2 = (rmin + (k - 1) * dr) ** 2
	    r[k] = sqrt (rap2)
	    w[k] = 0.0D0
	    do j = -n, n {
		y2 = (dy + j * d) ** 2
		do i = -n, n {
		    r2 = y2 + (dx + i * d) ** 2
		    if (r2 > rap2)
			break
		    w[k] = w[k] + a
		}
	    }
	}
end