procedure qdb5( pxm, nx, ny, x, y, bkg, bkgn)

# determine iterative median sky background level and number of pixels used

int	nx
int	ny
real	pxm[nx,ny]
int	x
int	y
double  mean
int     mean_n
double  bkg
int     bkgn

define	SIZE 56
real    bv[SIZE]
real    tmpv[SIZE]
int 	xxx[SIZE]
data  	xxx/       3, 4, 5, 6, 7, 8, 9,
                   3, 4, 5, 6, 7, 8, 9,
             1, 2,                     10,11,
             1, 2,                     10,11,
             1, 2,                     10,11,
             1, 2,                     10,11,
             1, 2,                     10,11,
             1, 2,                     10,11,
             1, 2,                     10,11,
                   3, 4, 5, 6, 7, 8, 9,
                   3, 4, 5, 6, 7, 8, 9/
int 	yyy[SIZE]
data	yyy/       1, 1, 1, 1, 1, 1, 1,
                   2, 2, 2, 2, 2, 2, 2,
	     3, 3,                      3, 3,
             4, 4,                      4, 4,
             5, 5,                      5, 5, 
             6, 6,                      6, 6, 
             7, 7,                      7, 7,
             8, 8,                      8, 8,
             9, 9,                      9, 9,
                  10,10,10,10,10,10,10,
                  11,11,11,11,11,11,11/
int	xx[11]
int	yy[11]
real	px
int	k
double  mdn
double  adv
int     tmp_n

begin

 xx[1] = x - 5
 xx[2] = x - 4
 xx[3] = x - 3
 xx[4] = x - 2
 xx[5] = x - 1
 xx[6] = x 
 xx[7] = x + 1
 xx[8] = x + 2
 xx[9] = x + 3
xx[10] = x + 4
xx[11] = x + 5

 yy[1] = y - 5
 yy[2] = y - 4
 yy[3] = y - 3
 yy[4] = y - 2
 yy[5] = y - 1
 yy[6] = y 
 yy[7] = y + 1
 yy[8] = y + 2
 yy[9] = y + 3
yy[10] = y + 4
yy[11] = y + 5

mean = 0
do k = 1, SIZE {
   px = pxm[xx[xxx[k]],yy[yyy[k]]]
   bv[k] = px
   mean = mean + px
}

call qdimdn( bv, SIZE, tmpv, tmp_n, mdn, adv)

mean_n = SIZE
mean = mean / mean_n

bkg = mdn
bkgn = tmp_n

end
#