/*      FILE:  phot_signl.c
 *   PURPOSE:  Calculate the signal within 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_CalcSgnl_f2 (
  qdphot5_ImageS *ImageS,
  qdphot5_PhotS *PhotS
){

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

  double sgnl;
#define SIZE 21
  int xxx[SIZE] = {
      2,3,4,
    1,2,3,4,5,
    1,2,3,4,5,
    1,2,3,4,5,
      2,3,4
  };
  int yyy[SIZE] = {
      1,1,1,
    2,2,2,2,2,
    3,3,3,3,3,
    4,4,4,4,4,
      5,5,5
  };
  int xx[5+1]; /* SPP->C: 0th element added to simplify conversion */
  int yy[5+1]; /* SPP->C: 0th element added to simplify conversion */
  int k;
  double pa=-180.;
  double m0;
  double mx;
  double my;
  double mxx;
  double myy;
  double mxy;
  double m00;
  double I[SIZE];
  double X[SIZE];
  double Y[SIZE];
  double xc;
  double yc;
  double e = -1.0;
  double top;
  double bottom;
  double **buff;
  double sharp=-1.;
  double peakpx;
  double ratio;
  double kappa;
  int x0ini;
  int y0ini;
  int x0wki;
  int y0wki;
  double sky;
  double x0wkd;
  double y0wkd;
  double x1outd;
  double y1outd;
  int delta_x0;
  int delta_y0;
  int big_drift;
  int drift;
  int pass;

  x0ini = PhotS->x0ini;
  y0ini = PhotS->y0ini;
  sky = PhotS->bkg;
  buff = ImageS->matrixd;
  kappa = -88.8888;

  x0wki = x0ini;
  y0wki = y0ini;
  pass = 0;

/*
 * ============================================================================
 * === STAGE ONE ==============================================================
 * ============================================================================
 */
stage_one:

  /* increment pass counter 
   */
  pass++; 

  /* put the zero-offset x coordinates in xx
   */
  xx[1] = x0wki - 2;
  xx[2] = x0wki - 1;
  xx[3] = x0wki;
  xx[4] = x0wki + 1;
  xx[5] = x0wki + 2;

  /* put the zero-offset y coordinates in yy
   */
  yy[1] = y0wki - 2;
  yy[2] = y0wki - 1;
  yy[3] = y0wki;
  yy[4] = y0wki + 1;
  yy[5] = y0wki + 2;

  peakpx = ((buff[y0wki][x0wki])/(ImageS->gain)) + (ImageS->bias);

  /* compute 1st order moments
   */
  m0 = 0.;
  mx = 0.;
  my = 0.;
  m00 = 0.;
  for (k=0; k<SIZE; ++k) { 
    X[k] = xx[xxx[k]] - x0wki;
    Y[k] = yy[yyy[k]] - y0wki;
    I[k] = buff[yy[yyy[k]]][xx[xxx[k]]] - sky;
    if (I[k]<0) I[k]=0.;
    m0 += I[k];
    mx += X[k]*I[k];
    my += Y[k]*I[k];
    m00 += I[k]*I[k];
  }
  sgnl = m0;
  /* 
   * abort if no signal!
   */
  if (sgnl<=0.) goto abort; 

  /* compute offsets in x and y
   */
  xc = mx / m0;
  yc = my / m0;
  /* 
   * add them to our working values
   */
  x0wkd = x0wki + (ImageS->llx0i) + xc;
  y0wkd = y0wki + (ImageS->lly0i) + yc;
  /*
   * compute indices
   */
  x0wki = x0wkd - (ImageS->llx0i) + 0.5;
  y0wki = y0wkd - (ImageS->lly0i) + 0.5;

  /* compute the shift of the central pixel of the aperture [integer pixels]
   */
  delta_x0 = x0wki - x0ini;
  delta_y0 = y0wki - y0ini;

  /* do we need to move the center of the aperture?
   */
  drift = ( (delta_x0!=0) || (delta_y0!=0) );
  /*
   * if not, goto stage_two
   */
  if (!drift) goto stage_two;

  /* do we need to move the center aperture by more than 1 px in x or y?
   */
  big_drift = ( (fabs(delta_x0)>1) || (fabs(delta_y0)>1) );
  /*
   * if so, abort
   */
  if (big_drift) goto abort;

  if (pass<2) goto stage_one;

/*
 * ============================================================================
 * === STAGE TWO ==============================================================
 * ============================================================================
 */
stage_two: 

  /* compute one-offset output coordinates */
  x1outd = x0wkd + 1.0;
  y1outd = y0wkd + 1.0;

  /* compute 2nd order moments */
  mxx = 0.;
  myy = 0.;
  mxy = 0.;
  for (k=0; k<SIZE; ++k) { 
    X[k] -= xc;
    Y[k] -= yc;
    mxx += (X[k]*X[k])*I[k];
    myy += (Y[k]*Y[k])*I[k];
    mxy += (X[k]*Y[k])*I[k];
  }
  mxx = mxx/m0;
  myy = myy/m0;
  mxy = mxy/m0;

  /* 
   * compute SHAPE (eccentricity) 
   *     and ANGLE (position angle)
   */
  bottom = (mxx+myy);
  if (bottom==0) goto bye;
  top = (mxx-myy)*(mxx-myy) + 4.0*(mxy*mxy);
  e = sqrt( top ) / bottom;                  /* eccentricity */
  pa = 0.5 * atan2 ( (2*mxy), (mxx - myy) ); /* position angle */
  pa *= (45.0/atan(1.0));                    /* convert to degrees */

  /* compute SHARP [ 0 (flat) <= SHARP <= 1 (spike)] */
  sharp = m00/(m0*m0);

  /* compute KAPPA [ ~0 for a critically-sampled normal distribution ] */
  top = I[10];               /* intensity of central pixel */
  bottom = (m0-I[10])/20.0;  /* average of all other pixels */
  ratio = (top/bottom)/3.36; /* ratio~1 for a normal distribution */
  if (ratio<=0) goto abort;
  kappa = log10(ratio);  /* kappa~0 for a normal distribution */
  if (kappa>9.999) kappa=9.999;
  if (kappa<-9.999) kappa=-9.999;

  goto bye;

abort:
  sgnl = 0.;

bye:
  PhotS->round = e;
  PhotS->angle = pa;
  PhotS->sharp = sharp;
  PhotS->kappa = kappa;
  PhotS->peakpx = peakpx;
  PhotS->sgnl = sgnl;
  PhotS->area = SIZE;
  PhotS->x1outd = x1outd;
  PhotS->y1outd = y1outd;

  return( status );

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