#include	"focas1.h"

/*  CONINDEX -- Compute mean surface brightness and concentration index.
    This program computes quantities similar to those described in
    "Morphological Classification of Galaxies with Simple Photometric
    Parameters", Doi, Fukugita, and Okamura.  The main difference is
    in the use of ellipitical apertures instead of circular apertures.
    Each object has a set of intensity weighted moments computed by
    FOCAS using pixels within a detection isophote.  The first order
    intensity weighted x and y moments define a centroid.  This is
    used for the aperture center.  The second order central intensity
    weighted moments, ixx, iyy, and ixy, define an elliptical light
    distribution.  The elliptical "isophotes" satisfy the equation:

	(1)  r**2 = iyy * x**2 + ixx * y**2 - 2 * ixy * x * y

    where x and y are the positions relative to the intensity centroid, and
    r**2 is a normalized radius which is constant on the elliptical
    isophote.  The moments are normalized in such a way that the radius is
    1 when the area of the isophote equals the isophotal area computed by
    FOCAS; i.e.	 the number of pixels within the detection isophote
    (which has some ragged shaped).

    The flux inside an ellipical isophote is obtained by specifying
    a normalized radius and summing the sky subtracted pixel values
    with values of equation (1) less than or equal to the specified
    radius.  In order to provided a finer measurement the pixels are
    divided into a set of equal valued subpixels (the pixel value
    divided by the number of subpixels).  Then if the center of a
    subpixel falls within the isophote it is included.	The current
    subsampling is 100 pixels or a 0.1 pixel boxes.  The area is
    the sum of the areas of the subpixels used.

    In this program two normalized radii are specified.	 One is usually
    1 (the elliptical isophote with area equal to the detection area)
    though some other number could be used.  The second is a fraction
    of this; for example 0.3.  The photometry results are entered
    into the catalog in the following fields:

	xx - the flux within the first elliptical isophote
	yy - the flux within the second elliptical isophote
	xy - the area of the first isophote
	cx - mean surface brightness of the first isophote = xx/xy
	cy - concentration index = yy/xx
*/

#define		SUBPIX		.1


static char *help[] = {
    "Compute mean surface brightness and concentration index.",
    "  usage: conindex catalog r1 r2 [filter]",
    "  Arguments: catalog - catalog of objects",
    "		  r1 - first normalized radius (usually 1.)",
    "		  r2 - second normalized radius (usually ~0.3)",
    "		  [filter] - filter options",
    "  Output: modified catalog fields cx, cy, xx, yy, and xy",
    "  Additional files required:  field file, area file",
    "",
    "  This should be applied only to brighter, larger objects say",
    "  with a filter such as A 100 10000.",
    0
};

main (argc, argv)
int	argc;
char   *argv[];
{
    double	R1, R2;

    struct objrec ob;
    struct areas  ar;
    struct image  pic;
    PIXEL	  *p;

    int		i;
    long	szrec, szobjrec();
    double	s, a, b, c, e, r, r1, r2, x, y, x1, y1;
    double	data1, data2, flux1, flux2, area, subarea, subpix;

    if ((argc < 4) || (argv[1][0] == '^')) {
	for (i = 0; help[i] != 0; i++)
	    fprintf (stderr, "%s\n", help[i]);
	exit (0);
    }

    /* Open the data files */
    cfd = catopen (argv[1], 2);
    pfd = fndfld (sp.ptfl, 0);
    afd = fndar (sp.arfl, 0);
    /* cmmnt (argc, argv); */
    wtcathdr (cfd, 0);
    szrec = szobjrec ();

    /* Get the isophote radii parameters */
    sscanf (argv[2], "%lf", &R1);
    sscanf (argv[3], "%lf", &R2);
    R1 *= R1;
    R2 *= R2;

    /* Initialize */
    subpix = 0.5 - SUBPIX / 2;
    subarea = SUBPIX * SUBPIX;
    setflt (argc - 4, &argv[4]);

    /* Loop through each object in the catalog */
    for (;;) {
	if (rdcatob (cfd, 0L, &ob))
	    break;
	if (ffilter (&ob) != 1)
	    continue;
	if (rdarea (afd, ob.entnum, ob.subent, ob.arpos, &ar))
	    continue;
	if (rdimg (pfd, ar.rct, &pic) < 0)
	    continue;
	p = pic.pic;

	/* Set the isophote parameters from the moments. */
	s = ob.ixx + ob.iyy;
	a = ob.iyy / s;
	b = ob.ixx / s;
	c = ob.ixy / s;
	e = hypot (2 * c, b - a);
	s = 2 * M_PI / ob.area / sqrt (1 - e * e);
	a *= s;
	b *= s;
	c *= -2 * s;

	/* Do the photometry */
	flux1 = 0.;
	flux2 = 0.;
	area = 0.;
	for (y = pic.y0-ob.icy; y < pic.y0+pic.dy-ob.icy; y++) {
	    for (x = pic.x0-ob.icx; x < pic.x0+pic.dx-ob.icx; x++) {
		data1 = (*p++) - ob.sbr;
		x1 = fabs (x) - subpix;
		y1 = fabs (y) - subpix;
		r1 = a * x1 * x1 + b * y1 * y1 + c * x1 * y1;
		if (r1 > R1)
		    continue;
		x1 = fabs (x) + subpix;
		y1 = fabs (y) + subpix;
		r2 = a * x1 * x1 + b * y1 * y1 + c * x1 * y1;
		if (r2 < R1) {
		    flux1 += data1;
		    area += 1.;
		} else {
		    data2 = data1 * subarea;
		    for (x1 = x-subpix; x1 <= x+.5; x1 += SUBPIX) {
			for (y1 = y-subpix; y1 <= y+.5; y1 += SUBPIX) {
			    r = a * x1 * x1 + b * y1 * y1 + c * x1 * y1;
			    if (r < R1) {
				flux1 += data2;
				area += subarea;
			    }
			}
		    }
		}

		if (r1 > R2)
		    continue;
		if (r2 < R2) {
		    flux2 += data1;
		} else {
		    data2 = data1 * subarea;
		    for (x1 = x-subpix; x1 <= x+.5; x1 += SUBPIX) {
			for (y1 = y-subpix; y1 <= y+.5; y1 += SUBPIX) {
			    r = a * x1 * x1 + b * y1 * y1 + c * x1 * y1;
			    if (r < R2)
				flux2 += data2;
			}
		    }
		}
	    }
	}

	/* Save the results */
	ob.xx = flux1;
	ob.yy = flux2;
	ob.xy = area;
	ob.cx = flux1 / area;
	ob.cy = flux2 / flux1;
	lseek (cfd, -szrec, 1);
	wtcat (&ob, cfd);
    }

    /* Finish up */
    imclose (pfd);
    fclose (afd);
    fclose (cfd);
}