/* MX LIBRARY FUNCTION:  mxqdphot0_qdpsap_vf11
 *
 * PURPOSE:  peaker "square" aperture photometry
 *
 * AUTHOR:  Kenneth J. Mighell (mighell@noao.edu)
 *
 * LANGUAGE:  ANSI C
 *
 * DOCUMENTATION:  http://www.noao.edu/staff/mighell/mx
 *
 * DATE:  20000316
 *
 * MOD:  2
 * 
 * Copyleft (L) 2000 Kenneth J. Mighell
 */

#include <mx.h>
#include "mxqdphot.h"

void mxqdphot0_qdpsap_vf11 (
double **buff1,
double **buff2,
int nx,
int ny,
int dx,
int dy,
double peakmin,
double peakmax,
double *outv,
int szoutv,
int *nstr_
){

  char mxfunc[] = "mxqdphot0_qdpsap_vf11";
  int status;

  int nstr;
  int x;
  int y;
  int xx;
  int yy;
  double pixel;
  double neighbor;
  int bingo;
  int count;
  int 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;
  int y0;
  int x0;
  int yy0;
  int xx0;
  int true = MX_TRUE;
  int false = MX_FALSE;
  int mxnstr;

  nstr = 0;
  mxnstr = szoutv/(2*11);

for (y=6; y<=(ny-5); ++y) {
y0 = y - 1;
for (x=6; x<=(nx-5); ++x) {
x0 = x - 1;
pixel = buff1[y0][x0];
if ( (pixel>=peakmin) && (pixel<peakmax) ) {

  bingo = true;
  for (yy=y-1; yy<=(y+1); ++yy) {
  yy0 = yy - 1;
  for (xx=x-1; xx<=(x+1); ++xx) {
    xx0 = xx - 1;
    if (bingo) {
      neighbor = buff1[yy0][xx0];
      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;
	  }
	} 
      } 
    } 
  } /* for xx */
  } /* for yy */

  if (bingo) {

/* 1st aperture ============================================================ */

    /* signal */
    mxqdphot0_qds2_vf5 (buff1, x, y, &sgnl1, &sgnl1_n);
    /* 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 */
    mxqdphot0_qdb5_vf5 (buff1, x, y, &bkg1, &bkg1_n);
    /* 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 */
    mxqdphot0_qds2_vf5 (buff2, x, y, &sgnl2, &sgnl2_n);
    /* 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 */
    mxqdphot0_qdb5_vf5 (buff2, x, y, &bkg2, &bkg2_n);
    /* 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 ================================================================= */


    /* accept star only if OK in both observations *and* nstr<=mxnstr */
    ok = (mg1<cmg1) && (mg2<cmg2) && (nstr<=mxnstr);
    if (ok) { 
      nstr = nstr + 1;
      xd = x + dx;
      yd = y + dy;
      index = ((nstr-1)*(2*11)) - 1;
      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;

    }



  } /* if bingo */

} /* if */
} /* for x */
} /* for y */

ok:
  *nstr_ = nstr; 
  return;

}
/* end-of-file */