/*      FILE: /mxtools/src/psfstar/g2d.c
 *   PURPOSE: Compute a 2-d Gaussian
 *    AUTHOR: Kenneth J. Mighell (mighell@noao.edu)
 *  LANGUAGE: ANSI C
 *      DATE: 2001NOV06
 * COPYRIGHT: (C) 2001 Assoc. of Universities for Research in Astronomy Inc.
 *
 * COMMENT:
 *  
 *  All X and Y coordinates in units of physical pixels, defined as follows:
 *
 *  The X coordinate of the   left edge of the lower left pixel is 0.0
 *  The X coordinate of the  right edge of the lower left pixel is 1.0
 *  The Y coordinate of the bottom edge of the lower left pixel is 0.0
 *  The Y coordinate of the    top edge of the lower left pixel is 1.0
 *
 *  In general, subtract 0.5 pixels from IRAF images coordinates:
 *    x = x_iraf - 0.5
 *    y = y_iraf - 0.5
 * 
 */

#include "mx.h"
#include "inc.h"

void Gaussian2d_v8( 
  double ig,      /*                volume of 2-d Gaussian */
  double xgp,     /* X physical coordinate of 2-d Gaussian */
  double ygp,     /* Y physical coordinate of 2-d Gaussian */
  double sg,      /*                 sigma of 2-d Gaussian */
  int    res_p,   /* resolution = sqrt(total number of subpixels) */
  double xp,      /* X physical coordinate of pixel */
  double yp,      /* Y physical coordinate of pixel */
  double *z_p     /* value of 2-d Gaussian at Z = Gaussian2d(X,Y) */
) {

  char func[] = "Gaussian2d_v8";

  double z;   

  double twopi;
  double c;

  double p;
  double dx;
  double dy;
  double r;
  double sg2;
  double five = 5.0;
  double offset;
  double delta;
  int resolution;
  int resolution2;
  double X;
  double Y;
  double dX;
  double dY;
  int ix;
  int iy;
  double R2;
  int res = res_p;

  twopi = 8.0*atan(1.0); 
  c = 1.0/twopi;
  sg2 = sg*sg;

  /* initialize output to zero */
  z = 0.0;

  /* negative sigmas? */
  if (sg<=0.0) goto bye;  /* abort */

  /* Are we more than 8 sigma away from the center of the 2-d Gaussian ? */
  dx = xp - xgp;
  dy = yp - ygp;
  r = sqrt( (dx*dx) + (dy*dy) );
  if ((r/sg)>five) goto bye; /* abort */

  res = 1;
  if (r<6.0) res++;
  if (r<4.5) res++;
  if (r<3.0) res++;
  if (r<1.5) res++;
  if (sg<0.8) res++;
  if (sg<1.2) res++;
  if (sg>2) res--;
  if (sg>4) res--;
  if (res<1) res=1;

  resolution = 1;
  if (res>resolution) resolution = res;
  if (resolution>100) resolution = 100; /* maximum is 10^4 subpixels */
  resolution2 = resolution * resolution;
  delta = 1.0/resolution;
  offset = -(resolution+1.0)/(2.0*resolution) ;
  for (iy=1; iy<=resolution; iy++) {
    Y = yp + offset + (iy*delta);
    dY = Y - ygp;
    for (ix=1; ix<=resolution; ix++) {
      X = xp + offset + (ix*delta);
      dX = X - xgp;
      R2 = (dX*dX) + (dY*dY);
      p = c * (ig/sg2) * exp(-R2/(2.0*sg2)) / resolution2;
      z += p;
    }
  }

bye:
  *z_p = z;
  return;

}
/* end-of-file */