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


/***************************************************************************
 * This program is intended to do the usual dark-subtraction and flat-fielding
 * for a drift-scan CCD image, all at once.  We could do the two steps
 * separately, but because XVista writes intermediate results to disk as
 * signed 16-bit ints, we'd risk losing one bit in the intmediate result.
 * So, this program can do it all, using 32-bit ints to keep full precision.
 * 
 * The idea is to make a dark vector and flatfield vector ahead of this,
 * then call like so:
 * 
 *    ccdvec file darkvect flatvect flat_mean row|col [outfile=]
 *
 * The program will go pixel-by-pixel through the image, calculating
 *
 *                       [                      flat_mean   ]
 *             result =  [ (pixel - darkvec) * -----------  ]  - 32768
 *                       [                       flatvec    ]
 *
 * Looking at the above equation from the inside out, and assuming that
 * all data values in the raw image, dark and flat vectors are in the
 * range (-32767 to +32767), we have
 *
 *                 (pixel - darkvec)      should range 0 - 65535
 *         
 *                     flat_mean
 *            mult by  ----------         should range 0 - 65535
 *                      flatvec
 *
 *             subtract 32768             should range -32767 to +32767
 * 
 * so we can write the result out as true 16-bit signed integers.
 */

#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[] = "ccdvec file darkvec flatvec flat_mean row|col [outfile=] help";

int
main(argc, argv)
int argc;
char *argv[];
{
	FITS_HANDLE fh, fdark, fflat, fout;
	static int i;
	static unsigned int nr, nc, sr, sc, nrow, ncol, row;
   static unsigned int dark_nrow, dark_ncol;
   static unsigned int flat_nrow, flat_ncol;
	static int16 data[NMAX], darkvecdata[NMAX], flatvecdata[NMAX];
	static int16 darkdatum, flatdatum;
	static double flat_mean;
	static long lresult[NMAX];
	static int do_row, do_col;
	char *gotit, outfile[NBUF];
	char errmsg[NBUF];
	int new_file=0;

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

	sr = 0;
	sc = 0;
	do_row = 0;
	do_col = 0;

	/* first arg must be file name */
	fh = fits_open(argv[1], "r", &nrow, &ncol);
	strcpy(outfile, argv[1]);
	nr = nrow;
	nc = ncol;	

	/* second arg must be name of dark vector file */
	fdark = fits_open(argv[2], "r", &dark_nrow, &dark_ncol);
	if ((dark_nrow > 1) && (dark_ncol > 1)) {
		error(0, "ccdvec: warning: dark vector has > 1 row and col");
	}

	/* third arg must be name of flat vector file */
	fflat = fits_open(argv[3], "r", &flat_nrow, &flat_ncol);
	if ((flat_nrow > 1) && (flat_ncol > 1)) {
		error(0, "ccdvec: warning: flat vector has > 1 row and col");
	}

	/* fourth arg must be mean of flatfield vector */
	flat_mean = (double) evaluate(argv[4]);
	if (flat_mean == 0) {
		error(1, "ccdvec: flat_mean value is 0");
	}
	if (flat_mean < 0) {
		error(0, "ccdvec: flat_mean value is negative");
	}

	/* fifth arg must be either "row" or "col" */
	if (strcmp(argv[5], "row") == 0) {
		do_row = 1;
	} 
	else if (strcmp(argv[5], "col") == 0) {
		do_col = 1;
	}	
	else {
		error(1, "ccdvec: fifth arg must be 'row' or 'col'");
	}

	for(i = 6; 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 (do_row + do_col != 1) {
		error(1, "ccdvec: you must specify exactly one of 'row' or 'col'");
	}
	if (do_row) {
		if ((nrow != dark_nrow) || (nrow != flat_nrow)) {
			error(-1, "ccdvec: number of rows in dark and/or flat doesn't match image");
		}
	}
	if (do_col) {
		if ((ncol != dark_ncol) || (ncol != flat_ncol)) {
			error(-1, "ccdvec: number of cols in dark and/or flat doesn't match image");
		}
	}

	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(fdark, 0, sc, darkvecdata, nc);
		fits_get_data(fflat, 0, sc, flatvecdata, nc);
	}

	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 (do_row) {
			fits_get_data(fdark, row, sc, &darkdatum, 1);
			fits_get_data(fflat, row, sc, &flatdatum, 1);
			if (flatdatum == 0) {
				sprintf(errmsg, "ccdvec: flatvec value in row %d is 0", row);
				error(0, errmsg);
				for (i = 0; i < nc; i++) {
					if (lresult[i] < 0) {
						lresult[i] = MIN_DATA_VAL;
					} 
					else {
						lresult[i] = MAX_DATA_VAL;
					}
				}
			}
			else {
				for (i = 0; i < nc; i++) {
					lresult[i] = (lresult[i] - darkdatum);
					lresult[i] *= (flat_mean/flatdatum);
					lresult[i] -= (long) 32768;
				}
			}
		}
		else { /* docol */
			for (i = 0; i < nc; i++) {
				if (flatvecdata[i] == 0) {
					sprintf(errmsg, "ccdvec: flatvec value in col %d is 0", i);
					error(0, errmsg);
					if (lresult[i] < 0) {
						lresult[i] = MIN_DATA_VAL;
					} 
					else {
						lresult[i] = MAX_DATA_VAL;
					}
				}
				else {
					lresult[i] = (lresult[i] - darkvecdata[i]);
					lresult[i] *= (flat_mean/flatvecdata[i]);
					lresult[i] -= (long) 32768;
				}
			}
		}
		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(fh);
	fits_close(fdark);
	fits_close(fflat);
	fits_close(fout);
	exit(0);
}