/**************************************************************************
 *                                                                        *
 * 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 the mean of data files

		MN source  [box=d] [NOZERO] [PIX=p] [circle=] [help]

	BOX       operates over box if specified
	CIRCLE    operates over circle if specified
	NOZERO    ignore zeroes
	PIX       use only every PIX'th pixel in calculations

	writes MEAN= and MEAN_source= and STD_DEV= into symbol table

	2/21/92 - exit(0) added -rrt
	9/92/92 - usage -rrt
	11/17/92 - removed call to exten -rrt
	9/22/94 - add circle -rrt
			- make output contain val=strings
	1/7/96  - added 'help' option
*/

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

int main(int, char **);

char *usage = "usage: mn file [box=] [pix=] [circle=] [nozero] [help]";

int 
main(argc, argv)
int argc;
char *argv[];
{
	int circle_flag=0;
	FITS_HANDLE fh;
	static int i, row;
	int circno,cr,cc,rad;
	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;
	double val;
	char *gotit, buf[NBUF], fname[NBUF];

	if (argc == 1) {
		error(-1, usage);
	}
	if (strcmp(argv[1], "help") == 0) {
		printf("%s\n", usage);
		exit(0);
	}

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

	for (i = 2; i < argc; i++) {
		if (keyword("help", argv[i])) {
			printf("%s\n", usage);
			exit(0);
		}
		if ((gotit = find("box", argv[i])) != NULL) {
			boxno = (int)(evaluate(gotit) + SMALL);
			if (getbox(boxno, &sr, &sc, &nr, &nc))
				error(-1, "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, "mn: 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;
		}
		fprintf(stderr, "mn: unknown argument '%s'\n",argv[i]);
		exit(1);
	}

	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];
				sum += val;		
				sumsq += (val*val);
				pixnum += 1.0;
			}
		}
	}

	if (pixnum == 0.0)
		error(-1, "mn: No pixels read");

	mean = sum/pixnum;
	diff = sumsq/pixnum - mean*mean;
	if (diff >= 0.0)
		std_dev = sqrt(diff);
	else {
		error(0, "mn: error in std. dev");	
		std_dev = -1.0;
	}
	printf("mean=%.2f rms=%.2f sum=%.0f npixels=%.0f\n",
			mean, std_dev, sum,pixnum);

	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);
	exit(0);
}