/**************************************************************************
 *                                                                        *
 * 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 SKY value in an image file by inspecting the
	histogram near the max value

                SKY file [box=] [circle=] [bin=] [pix=] [quiet] [help]
	
	puts SKY and SKYSIG symbols into the PCVISTA table.

	1/30/91 - added 'floor' to "sky" expression, and <math.h>   - MWR

	2/21/92 - added exit(0) -rrt

	6/2/92  - Installed new method of calculating sky sigma -rrt 
	6/4/92  - allow negative values of sky -rrt
	9/22/92 - removed OLD_SKY trash -rrt
	10/22/92- removed small bug from 'fsky()' - we were including the
	             'topbin' value returned by the histogram routines as valid,
	             but it contains the number of all pixels with values
	             > OR = the top bin - thus, it can easily overwhelm the
	             true sky peak if there is a large tail of bright pixels. MWR
	10/29/92 - sky used without arguments gives usage -rrt
	12/15/92 - rewrote entirely to use gaussian fit to histogram in order
	             to determine the sky sigma; much faster than old method.
	             Also, with new, much larger NHIST, no need for 'bin'
	             parameter, so I've deleted it. MWR 
    12/17/92 - skysig is printed with fewer digits -rrt
	9/28/94  - circle added
			 - output now has keyword=value strings
			 - error messages try to contain program name -rrt
	1/8/96   - added 'help' option. MWR
	2/19/96  - now use fit to gaussian to determine sky value, as
	             well as sky-sigma value.  MWR
	7/13/96  - added "bin" option, allowing user to set bin size in
	             histogram.  MWR
   2/23/03  - moved "find_sky" routine (and all support functions) 
	             to new library file "findsky.c",
	             so it is now available to other programs.  MWR
*/

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

#ifndef DEBUG
#undef DEBUG					/* define this to get some diag output */
#endif

#ifdef PROTO
int main(int, char **);
#else
#endif

#define BINSIZE   1        /* default binsize in sky histogram */
#define FRACTION  0.40     /* use only bins w/ at least this fraction */
#define MINBIN    5	   /* use at least this many bins around the peak */
                   	   /*   on either side of the peak, when fitting */
                   	   /*   a gaussian to find sky sigma             */


static int circle_flag;
static int cr,cc,rad;

#define USAGE "usage: sky file [box=] [circle=] [bin=] [pix=] [quiet] [help]"

int
main(argc, argv)
int argc;
char *argv[];
{
	int i, boxno=-1, nskip = 1, quiet_flag = 0, nrow, ncol,circno;
	int binsize = BINSIZE;
	int16 sky;
	FITS_HANDLE fh;
	static int nr, nc, sr, sc;
	char *gotit, buf[NBUF],fname[NBUF];
	float skysig;
	int gotfname=0;

	if (argc == 1) {
		error(-1, USAGE);
	}

	/* placate compiler */
	circno = 0;

	for (i = 1; i < argc; i++) {
		if (keyword("help", argv[i])) {
			fprintf(stdout,"%s\n",USAGE);
			exit(0);
		}
		if ((gotit = find("box", argv[i])) != NULL) {
			boxno = (int) (evaluate(gotit) + 0.5);
			continue;
		}
		if ((gotit = find("circle", argv[i])) != NULL) {
			circno = (int)(evaluate(gotit) + SMALL);
			if (getcircle(circno, &cr, &cc, &rad)){
				fprintf(stderr, "sky: circle %d not found\n",circno);
				exit(1);
			}
			circle_flag=1;
			continue;
		}
		if ((gotit = find("bin", argv[i])) != NULL) {
			binsize = (int) (evaluate(gotit) + 0.5);
			if (binsize < 1)
				error(-1, "bin parameter must be greater than or equal to 1");
			continue;
		}
		if ((gotit = find("pix", argv[i])) != NULL) {
			nskip = (int) (evaluate(gotit) + 0.5);
			if (nskip < 1)
				error(-1, "pix parameter must be greater than or equal to 1");
			continue;
		}
		if (keyword("quiet", argv[i])) {
			quiet_flag = 1;
			continue;
		}
		if(!gotfname){
			strcpy(fname,argv[i]);
			gotfname=1;
			continue;
		}
		fprintf(stderr,"sky: unknown argument '%s'\n",argv[i]);
		exit(1);
	}

	if(!gotfname)
		error(1,"no filename specified");

	fh = fits_open(fname, "r", &nrow, &ncol);
	nr = nrow;	
	nc = ncol;
	if(boxno !=-1){
		if(circle_flag)
			error(1,"sky: can't have circles and boxes");
		if (getbox(boxno, &sr, &sc, &nr, &nc))
			error(-1, "sky: box not found");
		check_box(boxno, nrow, ncol);
	}	
	if(circle_flag)
		check_circle(circno, nrow, ncol);
	
	if (find_sky(fh, sr, sc, nr, nc, nskip, binsize, &sky, &skysig) != 0) {
		error(0, "warning: sky or skysig may not be correct");
	}

	if (!quiet_flag) {
		printf("sky=%d skysig=%.2f\n", sky, skysig);
	}
	sprintf(buf, "sky=%d", sky);
	putsym(buf);
	sprintf(buf, "skysig=%.2f", skysig);
	putsym(buf);

	fits_close(fh);

	exit(0);
}