/*      FILE:  phot_bkg.c
 *   PURPOSE: 
 * 
 *  Calculate the robust median sky background level near the aperture.
 *
 *    AUTHOR:  Kenneth J. Mighell (mighell@noao.edu)
 *  LANGUAGE:  ANSI C
 *      DATE:  2001JUL03
 * COPYRIGHT:  (C) 2001 Assoc. of Universities for Research in Astronomy Inc.
 */

#include "mx.h"
#include "qdphot.h"

int qdphot5_PhotS_CalcBkg_f2 (
  qdphot5_ImageS *Image,
  qdphot5_PhotS *PhotS
){

  char mxfunc[] = "qdphot5_PhotS_CalcBkg_f2";
  int status = 0;

  double bkg = 0.;
  double ebkg = 0.;
  int bkgn;
#define	SIZE 56
  double bv[SIZE];
  double tmpv[SIZE];
  int xxx[SIZE] = {
          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] = {
          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+1]; /* SPP->C: 0th element added to simplify conversion */
  int yy[11+1]; /* SPP->C: 0th element added to simplify conversion */
  int k;
  double mdn;
  double adv;
  int tmp_n;
  double top;
  double bottom;
  double ron;
  int x0ini;
  int y0ini;

  x0ini = PhotS->x0ini;
  xx[1] = x0ini - 5;
  xx[2] = x0ini - 4;
  xx[3] = x0ini - 3;
  xx[4] = x0ini - 2;
  xx[5] = x0ini - 1;
  xx[6] = x0ini;
  xx[7] = x0ini + 1;
  xx[8] = x0ini + 2;
  xx[9] = x0ini + 3;
  xx[10] = x0ini + 4;
  xx[11] = x0ini + 5;

  y0ini = PhotS->y0ini;
  yy[1] = y0ini - 5;
  yy[2] = y0ini - 4;
  yy[3] = y0ini - 3;
  yy[4] = y0ini - 2;
  yy[5] = y0ini - 1;
  yy[6] = y0ini;
  yy[7] = y0ini + 1;
  yy[8] = y0ini + 2;
  yy[9] = y0ini + 3;
  yy[10] = y0ini + 4;
  yy[11] = y0ini + 5;

  for (k=0; k<SIZE; ++k) {
    bv[k] = Image->matrixd[yy[yyy[k]]][xx[xxx[k]]];
  }

  qdphot5_VecdRobustMdn_vf6 (bv, SIZE, tmpv, &tmp_n, &mdn, &adv); 

  bkgn = tmp_n;
  bkg = mdn;
  if (bkgn==0) {
    status = -1;
    goto bye;
  }

  /* 
   * error estimate (~1 sigma) of background 
   * 
   * Nota bene: 
   *
   *   average deviation = sqrt(2/pi) 
   *                     ~ 0.7989 standard deviations (1 sigma)
   *   with *no* sigma clipping
   *
   */
  ron = Image->ron;
  bottom = sqrt(bkgn);
  top = adv / 0.8;
  if (ron>top) top = ron;
  ebkg = top/bottom;
 
bye:
  PhotS->bkg_n = bkgn;
  PhotS->bkg = bkg;
  PhotS->ebkg = ebkg;
  return( status );

#undef SIZE
}
/* end-of-file */