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


/***************************************************************************
 * Divide an image by a "vector" (an image with one dimension equal to 1)
 *
 *    divvec file vector row|col [flat] [mul=] [outfile=]
 *
 * If 'flat' is given, then multiply result by mean of vector
 * If 'mul' is given, multiply each datum by the given constant value
 * If 'outfile', put result into the given file.
 *
 */

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


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

static char usage[] = 
   "divvec file vector row|col [flat] [mul=] [outfile=] [help]";

int
main(argc, argv)
int argc;
char *argv[];
{
	FITS_HANDLE fh, fvec,fout;
	static int i, nother;
	static unsigned int nr, nc, sr, sc, nrow, ncol, row, vec_nrow, vec_ncol;
	static int16 data[NMAX], vecdata[NMAX], datum;
	static int flat_flag, mul_flag;
	static double temp;
	static long mul = 1;
	static int lresult[NMAX];
	static int do_row, do_col;
	char *gotit, *meanstring, buf[NBUF], vecfile[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]);
	fvec = -1;
	nother = 0;
	flat_flag = 0;
	mul_flag = 0;
	nr = nrow;
	nc = ncol;	
	sr = 0;
	sc = 0;
	do_row = 0;
	do_col = 0;

	for(i = 2; i < argc; i++) {
		if (keyword("help", argv[i])) {
			printf("%s\n", usage);
			exit(0);
		}
		if ((gotit = find("outfile", argv[i])) != NULL) {
			strcpy(outfile,gotit);
			new_file=1;
			continue;
		}
		if (keyword("row", argv[i])){
			do_row = 1;
			continue;
		}
		if (keyword("col", argv[i])){
			do_col = 1;
			continue;
		}
		if (keyword("flat", argv[i])){
			flat_flag = 1;
			sprintf(buf, "mn_%s", vecfile);
			if ((meanstring = getsym(buf)) == NULL) {
				sprintf(buf, "divvec: can't find mean of '%s'", vecfile);
				error(1, buf);
			}
			sscanf(meanstring, "%lf", &temp);
			mul = (long) floor(temp + 0.5);
			printf("Mean is %ld \n", mul);
			continue;
		}
		if ((gotit = find("mul", argv[i])) != NULL) {
			mul_flag = 1;
			mul = (long) floor(evaluate(gotit) + 0.5);
			continue;
		}
		/* arg must be vector's file name if we got here */
		if (nother++ > 1)
			error(-1, "divvec: too many files given on command line");
		strcpy(vecfile, argv[i]);
		fvec = fits_open(vecfile, "r", &vec_nrow, &vec_ncol);
		if ((vec_nrow > 1) && (vec_ncol > 1)) {
			error(0, "divvec: warning: vector has > 1 row and col");
		}
	}
	if (nother != 1) {
		error(1, usage);
	}

	if (mul_flag + flat_flag > 1) {
		error(1, "divvec: can't use both 'mul' and 'flat' together");
	}

	if (do_row + do_col != 1) {
		error(1, "divvec: you must specify exactly one of 'row' or 'col'");
	}
	if (do_row) {
		if (nrow != vec_nrow) {
			error(-1, "divvec: number of rows not the same");
		}
	}
	if (do_col) {
		if (ncol != vec_ncol) {
			error(-1, "divvec: 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);
			
	if (do_col) {
		fits_get_data(fvec, 0, sc, vecdata, nc);
	}

	for (row = sr; row < (sr+nr); row++) {
		fits_get_data(fh, row, sc, data, nc);
		if (mul != 1) {
			for (i = 0; i < nc; i++) {
				lresult[i] = ((long)data[i])*mul;
			}
		}
		else {
			for (i = 0; i < nc; i++) {
				lresult[i] = data[i];
			}
		}
		if (do_row) {
			fits_get_data(fvec, row, sc, &datum, 1);
			for (i = 0; i < nc; i++) {
				lresult[i] /= (long) datum;
			}
		}
		else { /* docol */
			for (i = 0; i < nc; i++) {
				lresult[i] /= (long) vecdata[i];
			}
		}
		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(fvec);
	fits_close(fout);
	exit(0);
}