/**************************************************************************
 *                                                                        *
 * 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.                        *
 *                                                                        *
 **************************************************************************/


/*	DIVIDES data files
	
		DIV filea [fileb] [const=c] [flat] [mul=m] [dark=] [outfile=

	divides by constant if 'const' is specified
	multiplies by 'mul' if specified
	if 'flat' then mean of second data file is multiplied by at end
	BAD divisors are set to zero

	10/3/89 -div by 0 bug in case 9 fixed - RRT
	1/29/91 added 'floor' to expression for "const", <math.h>  - MWR
	9/9/92 - fixed bug that allowed values > MAX_DATA_VAL to be assigned
	         to output pixels (oops).  New calculation done entirely in
	         long, then checked at end to see if it's between the MIN
	         and MAX data values before being assigned.  Slows down
	         by about 10% on a Sun 1+, but we can live with it.  Besides,
	         without bug fix, peaks of stars, etc., get NEGATIVE values!
	         That's unacceptable!  - MWR

	9/15/92 - exit(0) added -rrt
	11/20/92 -added outfile option -rrt
	4/27/93  - bug in illegal option -rrt
*/

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

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

#define USAGE \
  "usage: div file [fileb] [const=] [flat] [mul=] [box=] [dark=] [outfile=] [help]"

int 
main(argc, argv)
int argc;
char *argv[];
{
	int i, row, boxno, dark = 0, op_flag;
	long constant = 1L, mul = 1L, lresult;
	int nrow, ncol, other_nrow, other_ncol;
	int sr = 0, sc = 0, nr, nc;
	FITS_HANDLE fh, fother = 0,fout;
	static int16 data[NMAX], other[NMAX];
	char *gotit, buf[NBUF], *meanstring, other_file[NBUF];
	char outfile[NBUF];
	int new_file=0;
	double temp;

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

	fh = fits_open(argv[1], "r", &nrow, &ncol);
	strcpy(outfile,argv[1]);
	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, "can't find box");
			check_box(boxno, nrow, ncol);			
			continue;
		}
		if ((gotit = find("const", argv[i])) != NULL) {
			constant = (long)(evaluate(gotit) + SMALL);
			if (constant == 0L)
				error(-1, "can't divide by zero");
			continue;
		}
		if ((gotit = find("outfile", argv[i])) != NULL) {
			strcpy(outfile,gotit);
			new_file=1;
			continue;
		}
		if ((gotit = find("dark", argv[i])) != NULL) {
			dark = (int)(evaluate(gotit) + SMALL);
			continue;
		}
		if ((gotit = find("mul", argv[i])) != NULL) {
			mul = (long)(evaluate(gotit) + SMALL);
			continue;
		}
		if (keyword("flat", argv[i])) {
			sprintf(buf, "mn_%s", other_file);
			if ((meanstring = getsym(buf)) == NULL) {
				sprintf(buf, "can't find mean of '%s'", other_file);
				error(-1, buf);
			}
			sscanf(meanstring, "%lf", &temp);
			mul = (long)floor(temp + 0.5);
			printf("Mean is %ld \n", mul);
			continue;
		}
		if (fother) {
			sprintf(buf,"unknown keyword '%s'", argv[i]);
			error(-1, buf);
		}
		fother = fits_open(argv[i], "r", &other_nrow, &other_ncol);
		strcpy(other_file, argv[i]);
	}
	if (fother) {
		if (other_nrow < (sr + nr))
			error(-1, "divisor file does not have enough rows");
		if (other_ncol < (sc + nc))
			error(-1, "divisor file does not have enough columns");
	}

	fout = fits_open(outfile, new_file?"w":"x", &nrow, &ncol);

	if(new_file)
		fits_copyheader(fh,fout, FITS_CHECK);

	for (i = 0; i < nc; i++)
		other[i] = 1;

	/* use this flag to determine just how much computation
	   we need to use - sometimes compute-bound on PCs */
	op_flag = ((mul != 1L) << 3) + ((constant != 1L) << 2) + 
	          ((dark != 0) << 1) + (fother > 0);

	for (row = sr; row < (sr + nr); row++) {
		fits_get_data(fh, row, sc, data, nc);
		if (fother)
			fits_get_data(fother, row, sc, other, nc);
		for (i = 0; i < nc; i++) {

			/* actually, the following messy construct DOES save 
			   15-25% execution time vs. one do-it-all statement,
			   so keep it in.... */

			switch (op_flag) {
			case 1:
				if (other[i] != 0)
					lresult = (long) data[i] / other[i];
				else
					lresult = MAX_DATA_VAL;
				break;
			case 3:
				if (other[i] != dark)
					lresult = ((long) data[i] - dark)/(other[i] - dark);
				else
					lresult = MAX_DATA_VAL;
				break;
			case 4:
				lresult = ((long) data[i] / constant);
				break;
			case 6:
				lresult = ((long)(data[i] - dark) / constant);
				break;
			case 9:
				if (other[i] != 0)
					lresult = (((long)mul*data[i])/other[i]);
				else
					lresult = MAX_DATA_VAL;
				break;
			case 11:
				if (other[i] != dark)
					lresult = ((mul*(long)(data[i] - dark))/(other[i] - dark));
				else
					lresult = MAX_DATA_VAL;
				break;
			case 12:
				lresult = ((mul*(long)data[i])/constant);
				break;
			case 14:
				if (other[i] != dark)	
					lresult = ((mul*(long)(data[i] - dark))/(other[i] - dark))/constant;
				else
					lresult = MAX_DATA_VAL;
				break;
			default:
				error(-1, "nonsensical options given - think it over and try again\n");
				continue;
			}
			if (lresult >= MAX_DATA_VAL)
				data[i] = MAX_DATA_VAL;
			else if (lresult <= MIN_DATA_VAL)
				data[i] = MIN_DATA_VAL;
			else
				data[i] = (int16) lresult;
		}
		fits_put_data(fout, row, sc, data, nc);
	}
	fits_close(fout);
	exit(0);
}