/* @(#)cjmagn.c	17.1.1.1 (ESO-DMD) 01/25/02 17:40:32 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
  MA 02139, USA.
 
  Correspondence concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.COPYRIGHT:  Copyright (c) 1987 European Southern Observatory,
                                         all rights reserved
.IDENTIFIER  Cjmagn
.LANGUAGE    C
.AUTHOR      K. Banse		ESO - Garching
.KEYWORDS    
.PURPOSE     calculate the magnitude of a star using method JMETH
.ALGORITHM   
             
.IN/OUTPUT
  call as    stat = Cjmagn( jmeth, arr, npix, ni, nb, xycen,
                            mag, dmag, sky, dsky, apix, flux )
  input:
          int   jmeth		method of magnitude calculation
                                 1: sum of 9 central pixels
                                 2: sum of all pixels in "flux" area (DEFAULT)
                                 3: circular aperture
          float *arr		data around the star
          int   npix[2]		dimension of ARR
          int   ni              number of pixels between "flux" area 
				and background
          int   nb		width of background along edge
          float xycen[2]	center of star in pixels in array "arr" (x,y)
  output:
          float *mag		magnitude of star
          float *dmag		uncertainty in magnitude determination
          float *sky		mean sky intensity
          float *dsky		uncertainty in sky determination
          float *nrpix          number of pixels used for magn. determination
          float *flux		total flux of the star

.RETURNS     status	 1: put star in center of image
                         0: ok
		        -1: array "arr" too small
			-2: center not in array "arr"
.ENVIRONment MIDAS
             #include <midas_def.h>      Prototypes for MIDAS interfaces
 
.VERSIONS    1.00       940411  from JMAGN.FOR    RvH
------------------------------------------------------------*/

/* Define _POSIX_SOURCE to indicate that this is a POSIX program */

#define  _POSIX_SOURCE 1


#include <math.h>
 
/*

*/

/*++++++++++++++++++++++++++++++
.IDENTIFIER  static function pixfrac
.PURPOSE     calculate the fraction of a pixel which is inside a circle
             with radius r around (xc,yc)
.IN/OUTPUT
  call as    static function
  input: 
          double f		pixel value (flux)
          double df[4]		neg. partial derivative of f with respect to x
          			pos. partial derivative of f with respect to x
          			neg. partial derivative of f with respect to y
          			pos. partial derivative of f with respect to y
          float  xdelt		xdelta (x - ycenter)
          float  ydelt		ydelta (y - ycenter)
          double r		radius of circle
  output:
          double *fi		flux inside circle
          double *pi		area inside circle
------------------------------*/

#ifdef __STDC__
      static void PFRAC(double f, double df[4], float xdelt, float ydelt,
                        double r, double *fi, double *pi)
#else
      static void PFRAC(f, df, xdelt, ydelt, r, fi, pi)
      float  xdelt, ydelt;
      double f, df[4], r, *fi, *pi;
#endif

{
double fx, fy, x, y, ai;
double v, vy, rx, ry, rp;
double dx = 0.1, dy = 0.1, da = 0.01;


*fi = *pi = 0.0;
ai = f - 0.5 * (df[0] - df[1] + df[2] - df[3]);
y  = -0.5 + (0.5 * dy);


/* Subdivide pixel and check */

sect_100:
if (y > 0.0)
   fy = df[2];
else
   fy = df[3];
x = -0.5 + (0.5 * dx);
vy = fy * y;
ry = (ydelt + y);
ry *= ry;

sect_200:
if (x > 0.0) 
   fx = df[0];
else
   fx = df[1];
v = vy + (fx * x);
rx = (xdelt + x);
rx *= rx; 
rp = sqrt(rx + ry);


/* check if inside */

if ((r - rp) >= 0.0)
   {
   *fi += v;
   *pi += da;
   }

x += dx;
if (x < 0.5) goto sect_200;

y += dy;
if (y < 0.5) goto sect_100;


*fi = (ai * *pi) + (*fi * da);
}

/*

*/

int Cjmagn( jmeth, arr, npix, ni, nb, pfac, xycen,
            mag, dmag, sky, dsky, nrpix, flux )
int   jmeth, *npix, ni, nb;
float *arr, *pfac, *xycen, *mag, *dmag, *sky, *dsky, *nrpix, *flux;

{
register int ii, ix, iy;
int   nn, ninb, nxb[2], nyb[2];
int   ierr, xmin, ymin, nrx, nry;

float  rn, rns, dflux;
float  *pntra, *pntr, a, xrval, yrval, fac;

double s, sk, sq, akap;


fac = *pfac;
*mag  = -9999;					/* set return values */
*dmag = *sky = *dsky = *flux = 0.0;
*nrpix = 0.0;
ierr = 0;
ninb = ni + nb;
nn = (2 * ninb) + 3;


/* check dimensions and position */

if ((npix[0] <nn) || (npix[1] < nn)) return (-3);	 /* array too small */

if ((*xycen < (ninb-1)) || (xycen[0] > (float)(npix[0]-ninb-1)))
   {
   if (jmeth == 1) return (-2);			/* center not in array */
   ierr = 1;					/* put star in center */
   xycen[0] = (npix[0]-1) * 0.5;
   }

if ((xycen[1] < (ninb-1)) || (xycen[1] > (float)(npix[1]-ninb-1))) 
   {
   if (jmeth == 1) return (-2);			/* center not in array */
   ierr = 1;					/* put star in center */
   xycen[1] = (npix[1]-1) * 0.5;
   }
  
 
/* for all methods we determine the sky intensity by clipping
   with 2.0 * sigma */

akap = 1.0e+30;				/* huge Kappa */
nxb[0] = nb - 1;
nxb[1] = npix[0] - nb;
nyb[0] = nb - 1;
nyb[1] = npix[1] - nb;


/* calculate background sky, if so required */

if (nb > 0)     
   {
   xrval = 0.0;				/* xrval = 0.0 (initial *sky) */

   for (ii=0; ii<10; ii++)			/* 10 iterations... */
      {
      s = 0.0;
      sq = 0.0;
      nn = 0;
      pntr = arr;

      for (ix=0; ix<npix[0]; ix++)
         {
         for (iy=0; iy<npix[1]; iy++)
            {
            if ((ix <= nxb[0]) || (ix >= nxb[1]) ||
                (iy <= nyb[0]) || (iy >= nyb[1]))
               {
               a = *pntr;
               if ( fabs((double) (a - xrval) ) <= akap )
                  {
                  nn++;
                  s  += a;
                  sq += (a*a);
                  }
               }
            pntr++;
            }
         }

      if (nn > 0)
         {
         rns = (float) (nn);
         xrval  = (float) (s / rns);
         a = sq/rns - (xrval*xrval);
         if (a > 0.0)
            {
            *dsky = (float) sqrt((double) a); 
            akap  = fac * (*dsky);
            }
         else
            {
            akap = *dsky = 0.0; 
            }
         }
      }

   *sky  = xrval;
   }


/* calculate the magnitude for a given method */

if (jmeth != 3) 
   {

   /* method #1: sum of 9 central pixels */
   
   if (jmeth == 1)
      {
      xmin = (int) floor((double) xycen[0] );	/* C indexing, so remove */
      nrx = xmin + 2;			/* subtraction of 1.0 from xycen */
      ymin = (int) floor((double) xycen[1] );
      nry = ymin + 2;
      }

   /* method #2: sum of all pixels within background */

   else   
      {
      xmin = ymin = ninb + 1;
      nrx = npix[0] - ninb;
      nry = npix[1] - ninb;
      }

   /* sum flux within area */

   nn = 0;
   s = 0.0;

   if ((xmin > nrx) || (ymin > nry)) return (-3);

   pntra = arr + (xmin-1 + (npix[0]*(ymin-1)));		/* point to start */
   for (iy=0; iy<(nry-ymin+1); iy++)
      {
      pntr = pntra;
      for (ix=0; ix<(nrx-xmin+1); ix++) 
         {
         nn ++;
         s += *pntr++;
         }
      pntra += npix[0];
      }

   rn = (float)(nn);
   }


/* method #3: circular aperture */

else
   {
   int      npxm1, npym1;
   float    ry;
   double   fi, r, dr, ro, ri, pi, val, dval[4];

   ro = npix[0] - xycen[0] - 1.0 - nb;
   ri = ro - ni;
   dr = sqrt((double)0.5);
   sk = 0.0;
   sq = 0.0;
   s = 0.0;
   rns = 0.0;
   rn = 0.0;
   npxm1 = npix[0] - 1;
   npym1 = npix[1] - 1;

   pntra = arr;
   for (iy=0; iy<npix[1]; iy++)
      {
      pntr = pntra;
      yrval = (float) (iy - xycen[1]);
      ry = yrval * yrval;
      for (ix=0; ix<npix[0]; ix++)
         {
         val = *pntr;
         xrval = (float) (ix-xycen[0]);
         r = sqrt((double) ((xrval*xrval) + ry));
         if ( (nb > 0) && (r >= ro) && (fabs(val - *sky) <= akap) )     
            {
            rns += 1.0;
            sk += val;
            sq += (val*val);
            }
               
         if ((r - dr) < ri) 
            {
            if ((r + dr) <= ri) 
               {
               s += val;
               rn += 1.0;
               }
            else
               {
               dval[0] = *(pntr+1) - val;
               a = val - *(pntr-1);
               if (ix == npxm1) dval[0] = a;
               dval[1] = a;
               if (ix == 0) dval[1] = dval[0];

               dval[2] = *(pntr+npix[0]) - val;
               a = val - *(pntr-npix[0]);
               if (iy == npym1) dval[2] = a;
               dval[3] = a;
               if (iy == 0) dval[3] = dval[2];

               PFRAC(val,dval,xrval, yrval,ri,&fi,&pi);
               s  += fi;
               rn += pi;
               }
            }
         pntr++;
         }
             
      pntra += npix[0];
      }
      
   if (rns > 1.0)
      { 
      xrval = (float) (sk / rns);
      *sky = xrval;
      a = (float) (sq/rns - (xrval*xrval));
      if (a > 0.0)
         *dsky = (float) sqrt((double) a);
      else
         *dsky = 0.0; 
      }
   }
      

/* For all methods: compute magnitude and uncertainty */

*nrpix = rn;

if (nb > 0)
   {
   *flux = s - (rn * *sky);
   dflux = rn * (*dsky) * (float) sqrt((double) (1.0/rn + 1.0/rns) );
   }
else
   {
   *flux = s;
   dflux = 0.;
   }

if (*flux >= (0.1*dflux))     
   {
   *mag  = -2.5 * log10( (double) *flux );
   *dmag = 1.0857362 * dflux / *flux;
   }
else
   ierr = 2;

return ierr;
}