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


/*	
	subtracts 'fileb' from 'filea', leaves result in 'filea'	

		SUB filea [fileb]  [const=c] [box=] [outfile=] [help]

	sub constant if const is specified

	1/29/91 - added 'floor' to expression for "const", <math.h>   - MWR
	9/9/92  - made arithmetic 32-bit, and added check against over/underflow
	          before converting back to 16-bit for output.  Now probably
	          quite a bit slower, unfortunately.   MWR
	9/22/92 - usage -rrt
	11/17/92 - outfile added -rrt
	1/10/96  - added 'help' option.  MWR
*/

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

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

char *usage = "usage: sub file1 [file2] [const=] [box=] [outfile=] [help]";

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

	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);
			if (getbox(boxno, &sr, &sc, &nr, &nc))
				error(-1, "can't find box");
			check_box(boxno, nrow, ncol);			
			continue;
		}
		if ((gotit = find("outfile", argv[i])) != NULL) {
			strcpy(outfile,gotit);
			new_file=1;
			continue;
		}
		if ((gotit = find("const", argv[i])) != NULL) {
			constant = (int) floor(evaluate(gotit) + 0.5);
			continue;
		}
		fother = fits_open(argv[i], "r", &other_nrow, &other_ncol);
		strcpy(other_file, argv[i]);
		if (nrow != other_nrow)
			error(-1, "number of rows not the same");
		if (ncol != other_ncol)
			error(-1, "number of columns not the same");
	}

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

	if(new_file)
		fits_copyheader(fh,fout, FITS_CHECK);
	for (row = sr; row < (sr + nr); row++) {
		fits_get_data(fh, row, sc, data, nc);
		for (i = 0; i < nc; i++) {
			lresult[i] = data[i];
		}
		if (fother) {
			fits_get_data(fother, row, sc, other, nc);
			for (i = 0; i < nc; i++) 
				lresult[i] -= other[i];
		}
		if (constant) {
			for (i = 0; i < nc; i++) 
				lresult[i] -= constant;
		}
		for (i = 0; i < nc; i++) {
			if (lresult[i] > MAX_DATA_VAL)
					data[i] = MAX_DATA_VAL;
				else if (lresult[i] < MIN_DATA_VAL)
					data[i] = MIN_DATA_VAL;
				else
					data[i] = (int16) lresult[i];
		}
		fits_put_data(fout, row, sc, data, nc);
	}
	fits_close(fout);
	exit(0);
}