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


/*	
	Collapses data to generate spectra

	  MASH file SP=i1,i2 [SP=i3,i4] [BK=i6,i7] [NORM] [COLS] [ROWS]

	makes columns of spectrum and background

	9/22/92 - output sent to stdout and not .spc
			- also changed usage, exit(0)
			- the upper limit is NOT included  -rrt
	7/22/93 - add box -rrt
*/

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

#ifdef PROTO
int main(int, char **);
static void mash_rows(int ,int);
static void mash_cols(int, int);
#else
static void mash_rows();
static void mash_cols();
#endif

#define NBIN 5

static int mean_flag;
static int col_flag;
static int nspect, nback;
static int nspect_total, nback_total;
static int slower[NBIN], supper[NBIN], blower[NBIN], bupper[NBIN];
static FITS_HANDLE fh;
static int16 data[NMAX];
static int nrow, ncol;

#define USAGE \
  "usage: mash file sp=i1,i2 [sp=i3,i4] [bk=i6,i7] [norm] [cols] [rows] \n      [box=] [help]"

int
main(argc, argv)
int argc;
char *argv[];
{
	int i, nfield, boxno=-1;
	char *gotit, fname[NBUF], string1[NBUF], string2[NBUF];
	int sr,sc,nr,nc;

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

	for (i = 1; i < argc; i++) {
		if ((gotit = find("sp", argv[i])) != NULL) {
			if (nspect >= NBIN)
				error(-1, "too many spectrum sections to handle");
			nfield = sscanf(gotit,"%[^,],%s", string1, string2);
			if (nfield != 2)
				error(-1, "syntax SP=i1,i2");
			slower[nspect] = (int) (evaluate(string1) + 0.5);
			supper[nspect] = (int) (evaluate(string2) + 0.5);
			nspect_total += (supper[nspect] - slower[nspect]);
			nspect++;
			continue;
		}
		if ((gotit = find("bk", argv[i])) != NULL) {
			if (nback >= NBIN)
				error(-1, "too many background sections to handle");
			nfield = sscanf(gotit,"%[^,],%s", string1, string2);
			if (nfield != 2)
				error(-1, "syntax BK=i1,i2");
			blower[nback] = (int) (evaluate(string1) + 0.5);
			bupper[nback] = (int) (evaluate(string2) + 0.5);
			nback_total += (bupper[nback] - blower[nback] );
			nback++;
			continue;
		}
		if (keyword("norm", argv[i])) {
			mean_flag = 1;
			continue;
		}
		if (keyword("rows", argv[i])) {
			col_flag = 0;
			continue;
		}
		if (keyword("cols", argv[i])) {
			col_flag = 1;
			continue;
		}
		if ((gotit = find("box", argv[i])) != NULL) {
			boxno = (int)evaluate(gotit);
			continue;
		}
		strcpy(fname,argv[i]);		/* generate output spectrum file name*/
	}

	fh = fits_open(fname, "r", &nrow, &ncol);
	if(boxno != -1){
		if (getbox(boxno, &sr, &sc, &nr, &nc))
				error(-1, "can't find box");
		check_box(boxno, nrow, ncol);
	}else{
		sr=0;
		sc=0;
		nr=nrow;
		nc=ncol;
	}
	if (nspect_total == 0)		/* correct for number of rows */
		error(-1,"no rows selected ");

	if ((nback_total != 0) && (nspect_total != nback_total))
		fprintf(stderr, 
			"warning... number of spectrum rows not equal to background rows");

	if (col_flag)
		mash_cols(sr,nr);
	else
		mash_rows(sc,nc);

	exit(0);
}

static void
mash_cols(sr,nr)
int sr,nr;
{
	int i, j,row, ibin;
	double spect_sum;

	for (row = sr , j=0 ; j < nr ; j++, row++) {
		fits_get_data(fh, row, 0, data, ncol);
		spect_sum = 0.0;
		for (ibin = 0; ibin < nspect; ibin++)
			for (i = slower[ibin]; i < supper[ibin]; i++)
				spect_sum += (double) data[i];
		for (ibin = 0; ibin < nback; ibin++)
			for (i = blower[ibin]; i < bupper[ibin]; i++)
				spect_sum -= (double) data[i];
		if (mean_flag)
			spect_sum /= nspect_total;
		printf("%d %8.2f\n", row, spect_sum);
	}
}

static void
mash_rows(sc,nc)
int sc,nc;
{
	int i, row, ibin;
	static double spect_sum[NMAX];

	for (ibin = 0; ibin < nspect; ibin++) {
		for (row = slower[ibin] ; row < supper[ibin]; row++) {
			fits_get_data(fh, row, 0, data, ncol);
			for (i = 0; i < nc; i++)
				spect_sum[i] += (double) data[i+sc];
		}
	}
	for (ibin = 0; ibin < nback; ibin++) {
		for (row = blower[ibin]; row < bupper[ibin]; row++) {
			fits_get_data(fh, row, 0, data, ncol);
			for (i = 0; i < nc; i++)
				spect_sum[i] -= (double) data[i+sc];
		}
	}
	for (i = 0; i < nc; i++) {
		if (mean_flag)
			spect_sum[i] /= nspect_total;
		printf("%d %8.2f\n", i+sc, spect_sum[i]);
	}
}