procedure qdpsap( buff1, buff2, nx, ny, dx, dy, peakmin, peakmax, 
   outv, szoutv, nstr)

int	nx
int	ny
real	buff1[nx,ny]
real	buff2[nx,ny]
int	dx
int	dy
real    peakmin
real    peakmax
int     szoutv
double  outv[szoutv]
int     nstr
int     mxnstr

#
int	x
int	y
int     xx
int     yy
real    pixel
real    neighbor
bool    bingo
int     count
bool    ok

double  sgnl1
int 	sgnl1_n
double  esgnl1
double  bkg1
int     bkg1_n
double  ebkg1
double  area1
double  flux1
double  mg1
double  emg1
double  cmg1

double  sgnl2
int 	sgnl2_n
double  esgnl2
double  bkg2
int     bkg2_n
double  ebkg2
double  area2
double  flux2
double  mg2
double  emg2
double  cmg2

double  xd
double  yd
int     fltrnm

int     index

#
begin

nstr = 0
mxnstr = szoutv/22
do y = 6, (ny-5) {
do x = 6, (nx-5) {
pixel = buff1[x,y]
if ( (pixel>=peakmin) && (pixel<peakmax) ) {

   bingo = true
   do yy = (y-1),(y+1) {
   do xx = (x-1),(x+1) {
      if (bingo) {
         neighbor = buff1[xx,yy]
         if ( neighbor > pixel ) {
            bingo = false
         } else if (neighbor == pixel) {
            if ( (xx!=x) || (yy!=y) ) {
               if (((xx<=x)&&(yy<y))||((xx>x)&&(yy<y))) bingo = false
            }
         }
      }
    }
    }
    if (bingo) {

### 1st aperture ##############################################################

# signal
       call qds2( buff1, nx, ny, x, y, sgnl1, sgnl1_n)
       if (sgnl1 < 0 ) {
          sgnl1 = 0
          esgnl1 = 0
       } else {
          esgnl1 = sqrt(sgnl1) # assumes *only* photon noise
       }

# background
       call qdb5( buff1, nx, ny, x, y, bkg1, bkg1_n)
       if (bkg1 < 0 ) {
          bkg1 = 0
          ebkg1 = 0
       } else {
          ebkg1 = sqrt(bkg1) # assumes *only* photon noise
       }

# magnitude
       area1 = sgnl1_n
       flux1 = sgnl1 - (bkg1*area1)
       mg1 = 99.999
       emg1 = 99.999
       if (flux1 > 0) {
          mg1 = -2.5*log10(flux1)
          emg1 = 1.0857*(esgnl1/sgnl1)
       }

# check magnitude
       flux1 = area1*(bkg1+ebkg1)
       cmg1 = 99.999
       if (flux1 > 0) cmg1 = -2.5*log10(flux1)

### 2nd aperture ##############################################################

# signal
       call qds2( buff2, nx, ny, x, y, sgnl2, sgnl2_n)
       if (sgnl2 < 0 ) {
          sgnl2 = 0
          esgnl2 = 0
       } else {
          esgnl2 = sqrt(sgnl2) # assumes *only* photon noise
       }

# background
       call qdb5( buff2, nx, ny, x, y, bkg2, bkg2_n)
       if (bkg2 < 0 ) {
          bkg2 = 0
          ebkg2 = 0
       } else {
          ebkg2 = sqrt(bkg2) # assumes *only* photon noise
       }

# magnitude
       area2 = sgnl2_n
       flux2 = sgnl2 - (bkg2*area2)
       mg2 = 99.999
       emg2 = 99.999
       if (flux2 > 0) {
          mg2 = -2.5*log10(flux2)
          emg2 = 1.0857*(esgnl2/sgnl2)
       }

# check magnitude
       flux2 = area2*(bkg2+ebkg2)
       cmg2 = 99.999

       if (flux2 > 0) cmg2 = -2.5*log10(flux2)

### results ###################################################################

       ok = (mg1<cmg1) && (mg2<cmg2) && (nstr<=mxnstr)
       if (ok) { # accept star only if OK in both observations

	  nstr = nstr + 1

	  xd = x + dx
	  yd = y + dy

          index = (nstr-1)*2*11

	  fltrnm = 1
	  outv[index+1] = xd
	  outv[index+2] = yd
	  outv[index+3] = nstr
	  outv[index+4] = fltrnm
	  outv[index+5] = mg1
	  outv[index+6] = emg1
	  outv[index+7] = sgnl1
	  outv[index+8] = area1
	  outv[index+9] = bkg1
	  outv[index+10]= ebkg1
	  outv[index+11]= bkg1_n

	  fltrnm = 2
	  outv[index+12]= xd
	  outv[index+13]= yd
	  outv[index+14]= nstr
	  outv[index+15]= fltrnm
	  outv[index+16]= mg2
	  outv[index+17]= emg2
	  outv[index+18]= sgnl2
	  outv[index+19]= area2
	  outv[index+20]= bkg2
	  outv[index+21]= ebkg2
	  outv[index+22]= bkg2_n

       }

    }
}
}
}

end
#