/**************************************************************************
 *                                                                        *
 * Copyright (c) 1996 Michael Richmond and Richard Treffers               *
 *                                                                        *
 *                    This software may be copied and distributed for     *
 *                    educational, research and not for profit services   *
 *                    provided that this copyright and statement are      *
 *                    included in all such copies.                        *
 *                                                                        *
 **************************************************************************/


/*	
	computes some statistics on a data file

		ABX source  [box=d] [NOZERO] [PIX=p] [sky=s]


	BOX       operates over box if specified
	NOZERO    ignore zeroes
	PIX       use only every PIX'th pixel in calculations
	SKY       subtract the supplied value from all pixels

	9/22/92 - call to exit(0) and minor cosmetics -rrt
	10/25/92- fixed a bug in column position of min, max vals when 
	              the ABX command is used on a box of an image.  MWR

	9/22/94 - circle added
	2/20/2000 - now write MEAN, MN_name and STD_DEV to symbol table.  MWR
*/

#include "pcvista.h"
#include "fits.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <unistd.h>


int 
main(argc, argv)
int argc;
char *argv[];
{
	FITS_HANDLE fh;
	char fname[256];
	char buf[256];
	int circle_flag=0;
	int circno,cr,cc,rad;
	static int i, row, first = 0, rmax, rmin, cmax, cmin;
	static int sc, sr, nr, nc, nozero_flag, nskip = 1, nrow, ncol, boxno;
	static int16 data[NMAX];
	static double sum, pixnum, mean, sumsq, diff, std_dev, sky;
	static double min, max;
	double val;
	char *gotit;

	if (argc == 1)
		error(-1, "usage: abx file [box=] [circle=] [pix=] [nozero] [sky=]");

	strcpy(fname, argv[1]);
	fh = fits_open(argv[1], "r", &nrow, &ncol);
	nr = nrow;
	nc = ncol;
	sky = 0.0;

	for (i = 2; i < argc; i++) {
		if ((gotit = find("box", argv[i])) != NULL) {
			boxno = (int)(evaluate(gotit) + SMALL);
			if (getbox(boxno, &sr, &sc, &nr, &nc))
				error(-1, "abx: box not found");
			check_box(boxno, nrow, ncol);
			continue;
		}
		if ((gotit = find("circle", argv[i])) != NULL) {
			circno = (int)(evaluate(gotit) + SMALL);
			if (getcircle(circno, &cr, &cc, &rad))
				error(-1, "abx: circle not found");
			check_circle(circno, nrow, ncol);
			circle_flag=1;
			continue;
		}
		if ((gotit = find("pix", argv[i])) != NULL) {
			nskip = (int)(evaluate(gotit) + SMALL);
			continue;
		}
		if (keyword("nozero", argv[i])) {
			nozero_flag = 1;
			continue;
		}
		if ((gotit = find("sky", argv[i])) != NULL) {
			sky = evaluate(gotit) + SMALL;
			continue;
		}
	}

	if(circle_flag){
		sr=cr-rad;
		nr=2*rad;
	}

	for (row = sr; row < (sr + nr); row++) {
		if(circle_flag){
			int chord;
			chord=find_chord(row-cr,rad);
			nc=2*chord;
			sc=cc-chord;
		}
		fits_get_data(fh, row, sc, data, nc);
		for (i = 0; i < nc; i += nskip) {
			if (data[i] || !nozero_flag) {
				val = (double) (data[i] - sky);
				if (first == 0) {
					min = val;
					max = val;
					first = 1;
				}
				else {
					if (val < min) {
						min = val;
						rmin = row;
						cmin = i + sc;
					}
					else if (val > max) {
						max = val;
						rmax = row;
						cmax = i + sc;
					}
				}
				sum += val;		
				sumsq += (val*val);
				pixnum += 1.0;
			}
		}
	}

	if (pixnum == 0.0)
		error(-1, "abx: no pixels read");

	mean = sum/pixnum;
	diff = sumsq/pixnum - mean*mean;
	if (diff >= 0.0)
		std_dev = sqrt(diff);
	else {
		error(0, "error in std. dev");	
		std_dev = -1.0;
	}

	/* write MEAN, STD_DEV, and MN_name into symbol table */
	sprintf(buf, "mn_%s=%9.2f", fname, mean);	/* update symbol table */
	putsym(buf);
	sprintf(buf, "mean=%9.2f", mean);
	putsym(buf);
	sprintf(buf, "std_dev=%9.2f", std_dev);
	putsym(buf);

	printf("tot. adus  %10.2f num pixels   %-10.0f\n", sum, pixnum);
	printf("mean adus  %10.2f std_dev adus %-10.2f\n", mean, std_dev);
	printf("max. adus  %10.2f at row= %4d  col= %4d \n", max, rmax, cmax);
	printf("min. adus  %10.2f at row= %4d  col= %4d \n", min, rmin, cmin);
	exit(0);
}