# include "data.h"

double get_aperture_stats (Matrix *matrix, int X, int Y, int Npix, int Nborder, double max) {

  double *ring;
  double x, y, x2, y2, I, sky, FWHMx, FWHMy, value, mag;
  int i, j, n, Npix2, Nring, Nmax;
  double Npts, gain, dsky2, dmag;
  char *string;
  
  string = get_variable ("GAIN");
  if (string == (char *) NULL) {
    fprintf (stderr, "assuming a value of 1.0\n");
    gain = 1.0;
  } else {
    gain = atof (string);
  }
  Nborder = MAX (1, Nborder);
  Nborder = MIN (1000, Nborder);
  
  Npix2 = (int)(0.5*Npix);
  Npix = 2 * Npix2 + 1;
  Nring = 4*Nborder*(Nborder + Npix);
  ALLOCATE (ring, double, Nring);
  bzero (ring, sizeof(double)*Nring);

  n = 0;  
  for (j = 0; j < Nborder; j++) {
    for (i = X - Npix2 - Nborder; i < X + Npix2 + Nborder + 1; i++, n+=2) {
      ring[n]   = fits_get_matrix_value (matrix, i, (int)(Y - Npix2 - j));
      ring[n+1] = fits_get_matrix_value (matrix, i, (int)(Y + Npix2 + j));
    }
    for (i = Y - Npix2; i < Y + Npix2 + 1; i++, n+=2) {
      ring[n]   = fits_get_matrix_value (matrix, (int)(X - Npix2 - j), i);
      ring[n+1] = fits_get_matrix_value (matrix, (int)(X + Npix2 + j), i);
    }
  }
  sort (ring, Nring);
  for (Npts = sky = dsky2 = 0, i = 0.25*Nring; i < 0.75*Nring; i++, Npts += 1.0) {
    sky += ring[i];
    dsky2 += ring[i]*ring[i];
  }
  sky = sky / Npts;
  dsky2 = dsky2 / Npts - sky*sky;
  free (ring);

  Npts = Nmax = 0;
  x = y = x2 = y2 = I = 0;
  for (i = X - Npix2; i < X + Npix2 + 1; i++) {
    for (j = Y - Npix2; j < Y + Npix2 + 1; j++) {
      value = fits_get_matrix_value (matrix, i, j);
      x  += i*(value - sky);
      y  += j*(value - sky);
      x2 += i*i*(value - sky);
      y2 += j*j*(value - sky);
      I  += (value - sky);
      Npts ++;
      if (value > max) {
	Nmax ++;
      }
    }
  }

  x = x / I;
  y = y / I;
  FWHMx = 2.355*sqrt (fabs(x2 / I - x*x));
  FWHMy = 2.355*sqrt (fabs(y2 / I - y*y));
  mag = -2.5*log10(I);
  dmag = sqrt (fabs(1.0 / (gain*I) + Npts*dsky2 / (I*I)));
  x = x + 1;
  y = y + 1;
  
  set_variable ("Xg", x);
  set_variable ("Yg", y);
  set_variable ("SXg", FWHMx);
  set_variable ("SYg", FWHMy);
  set_variable ("Sg", sky);
  set_variable ("dSg", sqrt (fabs (dsky2)));
  set_variable ("Zg", mag);
  set_variable ("dZg", dmag);
  set_variable ("Zcg", I);
  set_int_variable ("Nsat", Nmax);
  
  fprintf (GetOutfile(), "%f %f %f %f %f %f %f %f\n", x, y, FWHMx, FWHMy, sky, I, mag, dmag);

  return (mag);

}

