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


/*
 * Program to create "trails" in images, to simulate the effects
 * of imperfect charge transfer.  
 *
 * Usage:  trail  fitsfile  [outfile=]
 *
 * where "outfile" is optional name of output image with trails.
 * Default filename for output file is "trail.fts".
 *
 * MWR 1/31/1998
 */



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


/* 
 * We assume that charge is "trailed" uniformly over NTRAIL pixels.
 * ALPHA is the coefficient by which we multiply a pixel's value to 
 * determine how much charge moves into _each_ of the following NTRAIL pixels.
 */

#define NTRAIL      20
#define ALPHA        0.00036



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

#define USAGE \
  "usage: trail imagefile [outfile] [help]"


void
main(argc,argv)
int argc;
char *argv[];
{
	FITS_HANDLE fh, fout;
	static int i;
	static int nr, nc, sr, sc, nrow, ncol, row;
	int j, start_trail, end_trail, num_trail;
	long excess, total_excess;
	static int16 data[NMAX];
	static long lresult[NMAX];
	char *gotit,outfile[NBUF];

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

	fh = fits_open(argv[1], "r", &nrow, &ncol);
	strcpy(outfile, "trail.fts");
	nr = nrow;
	nc = ncol;	
	sr = 0;
	sc = 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);
			continue;
		}
		error(1, "trail: unknown arg");
	}

	fout = fits_open(outfile, "w", &nrow, &ncol);
	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];
		}

		/* 
		 * now move "charge" down the row -- we assume that
		 * pixels with small column values are read out first.
	 	 */
		for (i = 0; i < nc; i++) {
			excess = (long) ((floor) (lresult[i]*ALPHA));
			start_trail = (i + 1);
			end_trail = (i + NTRAIL);
			if (end_trail >= nc) {
				end_trail = (nc - 1);
			}
			num_trail = (end_trail - start_trail) + 1;
			total_excess = num_trail*excess;
			lresult[i] -= total_excess;

			for (j = start_trail; j <= end_trail; j++) {
				lresult[j] += excess;
			}
		}
		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(fout);
	exit(0);
}