/* This file contains:
	doDark
	NUVFactor
	AvgSciVal
*/

# include <stdio.h>
# include <stdlib.h>		/* calloc */
# include <math.h>		/* fabs */

# include <hstio.h>
# include "../stis.h"
# include "calstis1.h"
# include "../stiserr.h"
# include "../stistemperature.h"	/* defines DARKRATE */

static double NUVFactor (double);
static void AvgSciVal (SingleGroup *, short, float *);

/* This routine subtracts the dark image from x (in-place).
   For CCD data, the dark image is multiplied by the exposure time and
   divided by the gain before subtracting.  The dark time is just the
   exposure time; it DOES NOT INCLUDE the idle time since the last
   flushing of the chip or the readout time.

   For MAMA data, the dark image is just multiplied by the exposure time
   before subtracting.

   Phil Hodge, 1997 Oct 1:
	Change avg from 1 to 0, so bin2d will sum within bins instead of
	averaging.

   Phil Hodge, 1998 Feb 6:
	Remove the include for "../stissizes.h".

   Phil Hodge, 1998 Mar 13:
	Change calling sequence of DoppConv.

   Phil Hodge, 1998 Aug 6:
	Add border=0 to calling sequence of DoppConv.

   Phil Hodge, 1998 Oct 7:
	Include NUVFactor, and scale the dark by that factor.

   Phil Hodge, 1998 Nov 13:
	In NUVFactor, scale the expression from the Instrument Handbook
	by a constant, rather than by the mean of the dark image.  The
	calling sequence was changed, also.

   Phil Hodge, 1999 Mar 10:
	DARKRATE (used by NUVFactor) is now defined in ../stistemperature.h
	instead of in this file.
*/

int doDark (StisInfo1 *sts, SingleGroup *x, float *meandark) {

/* arguments:
StisInfo1 *sts     i: calibration switches, etc
SingleGroup *x    io: image to be calibrated; written to in-place
float *meandark    o: mean of dark image values subtracted
*/

	int status;

	SingleGroup y, z;	/* y and z are scratch space */
	float *ds;		/* Doppler smearing array */
	double factor = 1.;	/* scale factor, depending on temperature */
	int nds, d0;		/* size of ds, index of zero point */
	int extver = 1;		/* get this imset from dark image */
	int rx, ry;		/* for binning dark down to size of x */
	int x0, y0;		/* offsets of sci image */
	int same_size;		/* true if no binning of ref image required */
	int high_res;		/* true if high-res pixels in dispersion dir */
	int avg = 0;		/* bin2d should sum values within each bin */

	int FindBin (StisInfo1 *, SingleGroup *, SingleGroup *,
		int *, int *, int *, int *, int *, int *);
	int multk2d (SingleGroup *, float);
	int sub2d (SingleGroup *, SingleGroup *);
	int bin2d (SingleGroup *, int, int, int, int, int, SingleGroup *);
	int MakeDopp (double, double, double, double, double, int,
		float *, int *, int *);
	int DoppConv (SingleGroup *, int, float *, int, int);

	if (sts->darkcorr != PERFORM)
	    return (0);

	initSingleGroup (&y);

	/* Get the dark image data. */
	getSingleGroup (sts->dark.name, extver, &y);
	if (hstio_err())
	    return (OPEN_FAILED);

	/* Compare binning of science image and reference image;
	   get same_size and high_res flags, and get info about
	   binning and offset for use by bin2d.
	*/
	if (status = FindBin (sts, x, &y,
		&same_size, &high_res, &rx, &ry, &x0, &y0))
	    return (status);

	/* Do we need to do Doppler convolution? */
	if (sts->doppcorr == PERFORM) {

	    if (!high_res) {
		printf (
		"ERROR    Doppler convolution (DOPPCORR) was specified, \\\n");
		printf ("ERROR    but %s is binned to low-res pixels.\n",
			sts->dark.name);
		return (SIZE_MISMATCH);
	    }

	    /* Allocate space for the Doppler smearing array, making it
		larger than we will need.  The actual size nds will be
		updated by MakeDopp.
	    */
	    nds = 2 * (sts->doppmag + 1) + 21;
	    ds = (float *) calloc (nds, sizeof (float));

	    if (status = MakeDopp (sts->doppzero, sts->doppmag, sts->orbitper,
			sts->expstart, sts->exptime, sts->dispsign,
			ds, &nds, &d0))
		return (status);

	    /* Convolve y with the Doppler smearing function. */
	    if (status = DoppConv (&y, 0, ds, nds, d0))
		return (status);

	    free (ds);
	}

	/* Multiply the dark image by the exposure time and divide by the
	   atodgain (or just by exposure time for the MAMAs), and
	   subtract it from x.  For NUV MAMA data, there's an additional
	   scaling factor that depends on the temperature.
	*/

	if (same_size) {

	    /* No binning required. */

	    /* For NUV MAMA data, get scaling factor. */
	    factor = NUVFactor (sts->temperature);

	    if (status = multk2d (&y, factor * sts->exptime / sts->atodgain))
		return (status);
	    AvgSciVal (&y, sts->sdqflags, meandark);
	    if (status = sub2d (x, &y)) {
		printf ("ERROR    (darkcorr) size mismatch.\n");
		return (status);
	    }
	    freeSingleGroup (&y);

	} else {

	    /* Bin the dark image down to the actual size of x. */

	    initSingleGroup (&z);
	    allocSingleGroup (&z, x->sci.data.nx, x->sci.data.ny);
	    if (status = bin2d (&y, x0, y0, rx, ry, avg, &z)) {
		printf ("ERROR    (darkcorr) size mismatch.\n");
		return (status);
	    }
	    freeSingleGroup (&y);			/* done with y */
	    factor = NUVFactor (sts->temperature);
	    if (status = multk2d (&z, factor * sts->exptime / sts->atodgain))
		return (status);
	    AvgSciVal (&z, sts->sdqflags, meandark);
	    if (status = sub2d (x, &z))
		return (status);
	    freeSingleGroup (&z);			/* done with z */
	}

	if (sts->verbose && factor != 1.) {
	    printf (
"         Dark reference image was scaled by the factor %.6g, \\\n", factor);
	    printf ("         in addition to the exposure time.\n");
	}

	return (0);
}

/* The dark count rate for the NUV MAMA varies with temperature.  For this
   detector, sts->temperature was gotten earlier from a header keyword.  In
   this routine we will determine the scale factor by which the dark image
   should be scaled.  The factor will be one for other detectors or if
   the temperature could not be obtained from the header; those cases are
   flagged by a value of temperature less than or equal to zero.
   See the "MAMA Darks" section in Chapter 7 of the STIS Instrument Handbook.
*/

static double NUVFactor (double temperature) {

/* argument:
double temperature  i: from header keyword, but converted to Kelvin
the function value is the factor by which the dark image should be multiplied
*/

	double factor;

	if (temperature <= 0.)
	    factor = 1.;
	else
	    factor = (9.012e20 / DARKRATE) * exp (-12710. / temperature);

	return (factor);
}


/* This routine computes the average of the science data values
   that are not flagged as bad by the data quality array.
*/

static void AvgSciVal (SingleGroup *y, short sdqflags,
		float *meandark) {

	double sum;
	int numgood;			/* number of good pixels */
	int i, j;
	short flagval;			/* data quality flag value */

	numgood = 0;
	sum = 0.;
	for (j = 0;  j < y->sci.data.ny;  j++) {
	    for (i = 0;  i < y->sci.data.nx;  i++) {
		flagval = DQPix (y->dq.data, i, j);
		if ( ! (sdqflags & flagval) ) {
		    /* no serious flag bit set */
		    sum += Pix (y->sci.data, i, j);
		    numgood++;
		}
	    }
	}

	if (numgood > 0)
	    *meandark = sum / (double) numgood;
	else
	    *meandark = 0.;
}