/* N_REFDATA: Contains various routines for reading reference data,
** as well as initializing and freeing the reference data structures.
**
** Revision history:
** H.Bushouse	Sept. 1995	Build 1
** H.Bushouse	Aug.  1996	Rewritten for Build 2 (Version 2.0)
** H.Bushouse	27-Jan-1997	Modified n_getNlinData to check for dummy ref
**				file before reading ref data. Modified
**				n_getDarkImage to check for dummy ref file
**				before reading ref exposure time list. Modified
**				n_getPhotData and n_getPedigree to set the new
**				ref.dummy variable. Modified n_getNlinData to
**				handle new data types for NlinData structure.
**				(Version 2.1)
** H.Bushouse	25-Apr-1997	Modified n_getZoffImage to read SAMPTIME for
**				for zeroth-read. (Version 2.2)
** H.Bushouse	13-Jun-1997	n_getRefData loads NOISFILE for use by CRIDCALC
**				(Version 2.3)
** H.Bushouse	31-Jul-1997	Added group number to n_getDarkImage arguments
**				to support new n_zsigcalc routine (Version 3.0)
** H.Bushouse	18-Aug-1997	Added NlinData.zsci and zerr elements to
**				n_initRefData, n_freeRefData, and n_getNlinData
**				routines (Version 3.0)
** H.Bushouse	02-Oct-1997	n_getRefData no longer loads NOISFILE for use
**				by CRIDCALC (Version 3.0)
** H.Bushouse	14-Jan-1998	Fixed bug in n_getDarkImage that was causing a
**				memory error for MultiAccums with NSAMP=2:
**				instead of keying on "ngroup" for release of
**				memory for "etime", use "nic->group".
**				(Version 3.1.1)
** H.Bushouse	09-Mar-1998	Modified n_getDarkImage to look for exposure
**				time matches only to the nearest 10th of a
**				second (Version 3.2)
** H.Bushouse	05-Oct-1998	Modified n_getRefData to load ref files needed
**				by non-primary steps. Moved part of the dummy
**				ref file handling out of n_getPedigree and into
**				n_getRefData (Version 3.3)
** H.Bushouse	09-Feb-1999	Updated use of getKey routines for HSTIO v2.1
**				(Version 3.2.2)
** H.Bushouse	05-Aug-1999	Modified n_getNlinData to load up to 3 coeff
**				images and up to 6 error images, to accomodate
**				2nd-order correction scheme (Version 3.3)
** H.Bushouse	20-Aug-1999	Changed n_math fn's from type "int" to "void"
**				(Version 3.3)
*/

# include <math.h>
# include <ctype.h>
# include <stdio.h>
# include <stdlib.h>
# include <string.h>

# include <xtables.h>	/* defines TABLE I/O functions */
# include <hstio.h>	/* defines HST I/O functions */
# include "calnic.h"	/* defines NICMOS data structures */
# include "calnica.h"	/* defines CALNICA data structures */

# define  FATAL   1
# define  WARNING 2

int miss_file;		/* missing ref file counter */
int miss_kwd;		/* missing keyword counter */
int mismatch;		/* ref data mismatch counter */

/* N_GETREFDATA: Loads all reference data into memory. Checks for
** missing reference files and reference file data that do not match
** the observing mode of the science data being processed. Returns with
** an error status if any ref files are missing or incompatibile.
** All reference data are loaded when processing the first group of an
** input data file; only a new dark image is loaded when processing 
** subsequent groups of a MULTIACCUM file.
*/

int n_getRefData (NicInfo *nic, SingleNicmosGroup *zoff,
	SingleNicmosGroup *mask, SingleNicmosGroup *nois, NlinData *nlin,
	SingleNicmosGroup *dark, SingleNicmosGroup *flat, PhotData *phot) {

/* Arguments:
**	nic	io: NICMOS info structure
**	zoff	 o: MULTIACCUM zeroth-read image
**	mask	 o: bad pixel mask
**	nois	 o: read-noise image
**	nlin	 o: non-linearity image
**	dark	 o: dark current image
**	flat	 o: flat field image
**	phot	 o: photometry parameters
*/

	/* Function definitions */
	int n_getZoffImage (NicInfo *, SingleNicmosGroup *);
	int n_getMaskImage (NicInfo *, SingleNicmosGroup *);
	int n_getNoisImage (NicInfo *, SingleNicmosGroup *);
	int n_getDarkImage (NicInfo *, SingleNicmosGroup *, int);
	int n_getNlinData  (NicInfo *, NlinData *);
	int n_getFlatImage (NicInfo *, SingleNicmosGroup *);
	int n_getPhotData  (NicInfo *, PhotData *);

	/* Initialize the error counters */
	miss_file = 0;
	miss_kwd  = 0;
	mismatch  = 0;

	/* Load the DARK current image */
	if (nic->DARK.corr == PERFORM || nic->NOIS.corr == PERFORM ||
	   (nic->ZSIG.corr == PERFORM && nic->group == nic->ngroups)) {
	    if (n_getDarkImage (nic, dark, nic->group))
		return (status);
	    if (nic->DARK.ref.dummy) {
		if (nic->DARK.corr == PERFORM) {
		    nic->DARK.corr = SKIP;
		    sprintf (MsgText, "%s will be skipped", nic->DARK.swname);
		    n_warn  (MsgText);
		}
		if (nic->NOIS.corr == PERFORM) {
		    nic->NOIS.corr = SKIP;
		    sprintf (MsgText, "%s will be skipped", nic->NOIS.swname);
		    n_warn  (MsgText);
		}
		if (nic->ZSIG.corr == PERFORM) {
		    nic->ZSIG.corr = SKIP;
		    sprintf (MsgText, "%s will be skipped", nic->ZSIG.swname);
		    n_warn  (MsgText);
		}
	    }
	}

	/* If the static ref data have already been loaded, return.  */
	if (nic->RefDataLoaded)
	    goto ErrCheck;

	/* Load the MULTIACCUM zeroth-read image */
	if (nic->ZOFF.corr == PERFORM) {
	    if (n_getZoffImage (nic, zoff))
		return (status);
	}

	/* Load the bad pixel MASK image */
	if (nic->MASK.corr == PERFORM || nic->ZSIG.corr == PERFORM) {
	    if (n_getMaskImage (nic, mask))
		return (status);
	    if (nic->MASK.ref.dummy) {
		if (nic->MASK.corr == PERFORM) {
		    nic->MASK.corr = SKIP;
		    sprintf (MsgText, "%s will be skipped", nic->MASK.swname);
		    n_warn  (MsgText);
		}
		if (nic->ZSIG.corr == PERFORM) {
		    nic->ZSIG.corr = SKIP;
		    sprintf (MsgText, "%s will be skipped", nic->ZSIG.swname);
		    n_warn  (MsgText);
		}
	    }
	}

	/* Load the readNOISE image */
	if ((nic->NOIS.corr == PERFORM && nic->obsmode != RAMP) ||
	    nic->ZSIG.corr == PERFORM) {
	    if (n_getNoisImage (nic, nois))
		return (status);
	    if (nic->NOIS.ref.dummy) {
		if (nic->NOIS.corr == PERFORM) {
		    nic->NOIS.corr = SKIP;
		    sprintf (MsgText, "%s will be skipped", nic->NOIS.swname);
		    n_warn  (MsgText);
		}
		if (nic->ZSIG.corr == PERFORM) {
		    nic->ZSIG.corr = SKIP;
		    sprintf (MsgText, "%s will be skipped", nic->ZSIG.swname);
		    n_warn  (MsgText);
		}
	    }
	}

	/* Load the non-linearity coefficients */
	if (nic->NLIN.corr == PERFORM || nic->ZSIG.corr == PERFORM) {
	    if (n_getNlinData (nic, nlin))
		return (status);
	    if (nic->NLIN.ref.dummy) {
		if (nic->NLIN.corr == PERFORM) {
		    nic->NLIN.corr = SKIP;
		    sprintf (MsgText, "%s will be skipped", nic->NLIN.swname);
		    n_warn  (MsgText);
		}
		if (nic->ZSIG.corr == PERFORM) {
		    nic->ZSIG.corr = SKIP;
		    sprintf (MsgText, "%s will be skipped", nic->ZSIG.swname);
		    n_warn  (MsgText);
		}
	    }
	}

	/* Load the FLAT field image */
	if (nic->FLAT.corr == PERFORM) {
	    if (n_getFlatImage (nic, flat))
		return (status);
	    if (nic->FLAT.ref.dummy) {
		nic->FLAT.corr = SKIP;
		sprintf (MsgText, "%s will be skipped", nic->FLAT.swname);
		n_warn  (MsgText);
	    }
	}

	/* Load the photometry parameters */
	if (nic->PHOT.corr == PERFORM) {
	    if (n_getPhotData (nic, phot))
		return (status);
	    if (nic->PHOT.ref.dummy) {
		nic->PHOT.corr = SKIP;
		sprintf (MsgText, "%s will be skipped", nic->PHOT.swname);
		n_warn  (MsgText);
	    }
	}

	/* Load the background model parameters */
	if (nic->BACK.corr == PERFORM) {
	    /* no routine yet */
	}

ErrCheck:
	/* Error check */
	if (miss_file > 0) {
	    sprintf (MsgText, "Calibration reference file missing");
	    n_error (MsgText);
	    status = 1;
	} else if (mismatch > 0) {
	    sprintf (MsgText, "Incompatible reference data");
	    n_error (MsgText);
	    status = 1;
	} else if (miss_kwd > 0) {
	    status = 1;
	} else {

	    /* Set flag indicating that the static ref data have been loaded */
	    nic->RefDataLoaded = True;
	    status = 0;
	}

	/* return */
	return (status);
}
 
/* N_GETZOFFIMAGE: Load the data for a MULTIACCUM zeroth-read image.
** This is simply the last image group in the input file being processed.
** It is not necessary to check for PEDIGREE and DESCRIP keywords or to
** check the image size, as with normal reference images.
*/

int n_getZoffImage (NicInfo *nic, SingleNicmosGroup *zoff) {

/* Arguments:
**	nic	i: NICMOS info structure
**	zoff 	o: zeroth-read image data
*/

	/* Initialize the image data structure */
	initSingleNicmosGroup (zoff);

	/* Read the data */
	if (getSingleNicmosGroup (nic->ZOFF.ref.name, nic->nsamp, zoff))
	    status = 1;
	if (hstio_err() || status) {
	    n_filerr (nic->ZOFF.ref.name);
	    freeSingleNicmosGroup (zoff);
	    miss_file += 1;
	    return (status = 0);	/* don't abort yet */
	}

	/* Read the exposure time for the zeroth-read */
	nic->exptime[nic->nsamp-1] = 0;
	if (getKeyD (&(zoff->sci.hdr), "SAMPTIME",
		     &(nic->exptime[nic->nsamp-1]))) {
	    n_kwerr ("SAMPTIME", nic->ZOFF.ref.name);
	    return (status = 1);
	}

	/* Successful return */
	return (status = 0);
}

/* N_GETMASKIMAGE: Load the bad pixel MASK reference file image and
** check to make sure it's the appropriate file for using with the
** science data being processed.
*/

int n_getMaskImage (NicInfo *nic, SingleNicmosGroup *mask) {

/* Arguments:
**	nic	i: NICMOS info structure
**	mask	o: bad pixel mask image
*/

	/* Local variables */
	Bool missing;		/* missing file flag */

	/* Function definitions */
	int n_getRefImage (CalStep *, int, SingleNicmosGroup *, Bool *);
	void checkKeyI (char *, Hdr *, char *, int, int);

	/* Read the reference image */
	if (n_getRefImage (&nic->MASK, 1, mask, &missing))
	    return (status);

	/* Check for a missing file */
	if (missing) {
	    miss_file += 1;
	    return (status = 0);	/* don't abort yet */
	}

	/* Check the MASK file CAMERA number */
	checkKeyI (nic->MASK.ref.name, mask->globalhdr, "CAMERA", nic->camera,
		   FATAL);

	/* Successful return */
	return (status = 0);
}

/* N_GETNOISIMAGE: Load the readnoise reference file image and
** check to make sure it's the appropriate file for using with the
** science data being processed.
*/

int n_getNoisImage (NicInfo *nic, SingleNicmosGroup *nois) {

/* Arguments:
**	nic	i: NICMOS info structure
**	nois	o: read noise image
*/

	/* Local variables */
	Bool missing;		/* missing file flag */

	/* Function definitions */
	int n_getRefImage (CalStep *, int, SingleNicmosGroup *, Bool *);
	void checkKeyI (char *, Hdr *, char *, int, int);
	void checkKeyS (char *, Hdr *, char *, char *, int);

	/* Read the reference image */
	if (n_getRefImage (&nic->NOIS, 1, nois, &missing))
	    return (status);

	/* Check for a missing file */
	if (missing) {
	    miss_file += 1;
	    return (status = 0);	/* don't abort yet */
	}

	/* Check the NOISE file CAMERA number */
	checkKeyI (nic->NOIS.ref.name, nois->globalhdr, "CAMERA", nic->camera,
		   FATAL);

	/* Check the NOISE file READOUT speed */
	checkKeyS (nic->NOIS.ref.name, nois->globalhdr, "READOUT", nic->readout,
		   FATAL);

	/* Check the NOISE file NREAD value */
	checkKeyI (nic->NOIS.ref.name, nois->globalhdr, "NREAD", nic->nread,
		   WARNING);

	/* Successful return */
	return (status = 0);
}

/* N_GETNLINDATA: Read linearity coefficients and their associated
** data quality flags and node values from the NLINFILE reference file.
** The coefficient, DQ, and node arrays are returned in a single NlinData 
** structure. The super zero read SCI and ERR arrays are also read and loaded
** into the NlinData structure. Also check to make sure that this reference
** file is the appropriate one to use with the science data being processed.
**
** Initially the correction used was a 1st-order polynomial, involving 2
** coefficient and 3 error images. This routine now accomodates a 2nd-order
** correction, using 3 coefficient and 6 error images. It is designed to be
** backwards-compatible with the 1st-order reference files by checking to
** see how many coeff arrays are present. Unused coeff and error arrays are
** zero-filled.
*/

int n_getNlinData (NicInfo *nic, NlinData *nlin) {

/* Arguments:
**	nic 	 i: NICMOS info structure
**	nlin	 o: data structure containing linearity coeff's and DQ flags
*/

	/* Local variables */
	int i, j, k;		/* loop indexes */
	IODescPtr im;		/* image pointer */
	int ncoeff = 0;		/* number of coeff arrays */

	/* Function definitions */
	int n_getPedigree (Hdr *, CalStep *);
	void checkKeyI (char *, Hdr *, char *, int, int);

	/* Allocate memory for the nlin data structure */
	nlin->coeff = (FloatHdrData *)calloc(3,sizeof(FloatHdrData));
	nlin->error = (FloatHdrData *)calloc(6,sizeof(FloatHdrData));
	nlin->dqual = (ShortHdrData *)calloc(1,sizeof(ShortHdrData));
	nlin->nodes = (FloatHdrData *)calloc(2,sizeof(FloatHdrData));
	nlin->zsci  = (FloatHdrData *)calloc(1,sizeof(FloatHdrData));
	nlin->zerr  = (FloatHdrData *)calloc(1,sizeof(FloatHdrData));
	if (nlin->coeff==NULL || nlin->error==NULL || nlin->dqual==NULL ||
	    nlin->nodes==NULL || nlin->zsci ==NULL || nlin->zerr ==NULL) {
	    sprintf (MsgText, "Can't allocate memory for NLIN ref data");
	    n_error (MsgText);
	    return (status = 1);
	}

	/* Open the nonlinearity reference file */
	im = openInputImage (nic->NLIN.ref.name, "", 0);
	if (hstio_err()) {
	    n_openerr (nic->NLIN.ref.name);
	    miss_file += 1;
	    return (status = 0);	/* don't abort yet */
	}

	/* Allocate memory for the file's primary header */
	nlin->globalhdr = (Hdr *)calloc(1,sizeof(Hdr));
	if (nlin->globalhdr == NULL) {
	    sprintf (MsgText, "Can't allocate memory for NLIN ref file header");
	    n_error (MsgText);
	    return (status = 1);
	}

	/* Read the primary header */
	getHeader (im, nlin->globalhdr);
	if (hstio_err()) {
	    n_readerr (nic->NLIN.ref.name);
	    closeImage (im);
	    freeHdr (nlin->globalhdr);
	    return (status = 1);
	}

	/* Close the image */
	closeImage (im);

	/* Read the reffile PEDIGREE and DESCRIP keywords */
	if (n_getPedigree (nlin->globalhdr, &nic->NLIN)) {
	    freeHdr (nlin->globalhdr);
	    return (status);
	}

	/* Check the NLIN file CAMERA number */
	checkKeyI (nic->NLIN.ref.name, nlin->globalhdr, "CAMERA", nic->camera,
		   FATAL);

	/* Check for a dummy file */
	if (nic->NLIN.ref.dummy)
	    return (status = 0);

	/* Read the NEXTEND keyword from the NLINFILE */
        if (getKeyI (nlin->globalhdr, "NEXTEND", &ncoeff)) {
            n_kwerr ("NEXTEND", nic->NLIN.ref.name);
            miss_kwd += 1;
            return (status = 0);
        }

	/* Determine the number of coefficient and error arrays present:
	** if nextend <= 10, assume there are 2 coeffs and 3 errors;
	** else assume there are 3 coeffs and 6 errors. */
	if (ncoeff <= 10)
	    ncoeff = 2;
	else
	    ncoeff = 3;

	/* Get the coefficient images from the NLINFILE */
	for (i=0; i < ncoeff; i++) { 
	     initFloatHdrData (&nlin->coeff[i]);
	     if (getFloatHD (nic->NLIN.ref.name, "COEF", i+1, &nlin->coeff[i]))
		 return (status=1);
	}

	/* Initialize undefined coeffs (up to 3) to zero */
	for (i=ncoeff; i < 3; i++) {
	     initFloatHdrData (&nlin->coeff[i]);
	     allocFloatHdrData (&nlin->coeff[i], SZ_NICIMG, SZ_NICIMG);
	     for (k=0; k < SZ_NICIMG; k++) {
		  for (j=0; j < SZ_NICIMG; j++)
		       Pix(nlin->coeff[i].data,j,k) = 0.0;
	     }
	}	

	/* Get the error images from the NLINFILE */
	ncoeff = (ncoeff + ncoeff*ncoeff) / 2;
	for (i=0; i < ncoeff; i++) {
	     initFloatHdrData (&nlin->error[i]);
	     if (getFloatHD (nic->NLIN.ref.name, "ERR", i+1, &nlin->error[i]))
		 return (status=1);
	}

	/* Initialize undefined errors (up to 6) to zero */
	for (i=ncoeff; i < 6; i++) {
	     initFloatHdrData (&nlin->error[i]);
	     allocFloatHdrData (&nlin->error[i], SZ_NICIMG, SZ_NICIMG);
	     for (k=0; k < SZ_NICIMG; k++) {
		  for (j=0; j < SZ_NICIMG; j++)
		       Pix(nlin->error[i].data,j,k) = 0.0;
	     }
	}	

	/* Get the DQ array from the NLINFILE */
	initShortHdrData (&nlin->dqual[0]);
	if (getShortHD (nic->NLIN.ref.name, "DQ", 1, &nlin->dqual[0]))
	    return (status=1);

	/* Get the node arrays from the NLINFILE */
	for (i=0; i<2; i++) {
	     initFloatHdrData (&nlin->nodes[i]);
	     if (getFloatHD (nic->NLIN.ref.name, "NODE", i+1, &nlin->nodes[i]))
		 return (status=1);
	}

	/* Get the super zero-read SCI array from the NLINFILE (Vsn 3.0) */
	initFloatHdrData (&nlin->zsci[0]);
	if (getFloatHD (nic->NLIN.ref.name, "ZSCI", 1, &nlin->zsci[0]))
	    return (status=1);

	/* Get the super zero-read ERR array from the NLINFILE (Vsn 3.0) */
	initFloatHdrData (&nlin->zerr[0]);
	if (getFloatHD (nic->NLIN.ref.name, "ZERR", 1, &nlin->zerr[0]))
	    return (status=1);

	/* Successful return */
	return (status = 0);
}

/* N_GETFLATIMAGE: Load the flat field reference file image and
** check to make sure it's the appropriate file to use with the
** science data being processed.
*/

int n_getFlatImage (NicInfo *nic, SingleNicmosGroup *flat) {

/* Arguments:
**	nic	i: NICMOS info structure
**	flat	o: flat field image
*/

	/* Local variables */
	Bool missing;		/* missing file flag */

	/* Function definitions */
	int n_getRefImage (CalStep *, int, SingleNicmosGroup *, Bool *);
	void checkKeyI (char *, Hdr *, char *, int, int);
	void checkKeyS (char *, Hdr *, char *, char *, int);

	/* Read the reference image */
	if (n_getRefImage (&nic->FLAT, 1, flat, &missing))
	    return (status);

	/* Check for a missing file */
	if (missing) {
	    miss_file += 1;
	    return (status = 0);	/* don't abort yet */
	}

	/* Check the FLAT file CAMERA number */
	checkKeyI (nic->FLAT.ref.name, flat->globalhdr, "CAMERA", nic->camera,
		   FATAL);

	/* Check the FLAT file FILTER name */
	checkKeyS (nic->FLAT.ref.name, flat->globalhdr, "FILTER", nic->filter,
		   FATAL);

	/* Successful return */
	return (status = 0);
}

/* N_GETDARKDATA: Load a dark current image from the DARKFILE reference file.
** A single image is loaded if the exposure time of the data being calibrated
** matches the exposure time of one of the dark images. Otherwise if the
** exposure time lies between two of the available darks, two darks are
** loaded and a single image is computed by interpolating between them.
** If the exposure time lies outside the range of the available darks, an
** image is computed by extrapolating from the closest available dark.
*/

# define ALLOWDIFF 0.1	/* Max allowed exptime difference */

int n_getDarkImage (NicInfo *nic, SingleNicmosGroup *dark1, int ngroup) {

/* Arguments:
**	nic	 i: NICMOS info structure
**	dark1	 o: dark current image
**	ngroup	 i: group number
*/

	/* Local variables */
	int i;					/* loop index */
	char kword[8+1];			/* keyword string */
	static int ndarks;			/* number of dark images */
	double etime_lower, etime_upper;	/* dark exposure times */
	static double *etime;			/* sci image exposure time */
	SingleNicmosGroup dark2;		/* temporary dark image */
	Bool missing;				/* missing ref file flag */

	/* Function definitions */
	int  n_getRefImage (CalStep *, int, SingleNicmosGroup *, Bool *);
	void checkKeyI (char *, Hdr *, char *, int,    int);
	void checkKeyS (char *, Hdr *, char *, char *, int);
	void dark_interp (SingleNicmosGroup *, SingleNicmosGroup *, 
			  double, double, double);
	void dark_extrap (SingleNicmosGroup *, double, double);

	/* Initialize frame interpolation information */
	nic->DarkType   = 0;
	nic->darkframe1 = 0;
	nic->darkframe2 = 0;

	/* When processing group 1, load the list of ref exposure times */
	if (nic->RefDataLoaded == False) {

	    /* Read the first group of the dark reference file */
	    if (n_getRefImage (&nic->DARK, 1, dark1, &missing))
		return (status);

	    /* Check for a missing file */
	    if (missing) {
		miss_file += 1;
		return (status = 0);	/* don't abort yet */
	    }

	    /* Check the DARK file CAMERA number */
	    checkKeyI (nic->DARK.ref.name, dark1->globalhdr, "CAMERA",
		       nic->camera, FATAL);

	    /* Check the DARK file NREAD number */
	    checkKeyI (nic->DARK.ref.name, dark1->globalhdr, "NREAD",
		       nic->nread, WARNING);

	    /* Check the DARK file SAMP_SEQ value */
	    checkKeyS (nic->DARK.ref.name, dark1->globalhdr, "SAMP_SEQ",
		       nic->sampseq, FATAL);

	    /* Check for a dummy file */
	    if (nic->DARK.ref.dummy)
		return (status = 0);

	    /* Find out how many dark images are in this file */
	    ndarks = 0;
	    if (getKeyI (dark1->globalhdr, "NUMEXPOS", &ndarks)) {
		n_kwerr ("NUMEXPOS", nic->DARK.ref.name);
		return (status = 1);
	    }

	    /* Allocate memory for the list of exposure times */
	    if (etime != NULL) free(etime);
	    if ((etime = (double *)calloc(ndarks, sizeof(double))) == NULL) {
		sprintf (MsgText, "Can't allocate memory in getDarkData");
		n_error (MsgText);
		return (status = 1);
	    }

	    /* Read the list of exposure times */
	    for (i = 0; i < ndarks; i++) {
		 sprintf (kword, "EXPOS_%d\0", i+1);
		 etime[i] = 0;
		 if (getKeyD (dark1->globalhdr, kword, &(etime[i]))) {
		     n_kwerr (kword, nic->DARK.ref.name);
		     return (status = 1);
		 }
	    }
	}

	/* If the exposure time is zero, don't bother loading data */
	if (nic->exptime[ngroup-1] <= 0)
	    return (status = 0);

	/* Find the appropriate match between science and dark exp time */
	if (nic->obsmode == MULTIACCUM &&
	    strncmp (nic->sampseq, "NONE", 4) != 0) {

	    /* If the science data are from a MultiAccum standard exposure
	    ** time sequence, then find the ref file group that matches the
	    ** exposure time of the science group. */
	    for (i = 0; i < ndarks; i++) {
		 if (fabs(etime[i] - nic->exptime[ngroup-1]) <= ALLOWDIFF) {
		     nic->DarkType = MATCH;
		     nic->darkframe1 = i+1;
		     break;
		 }
	    }

	    /* Make sure we found a match */
	    if (nic->darkframe1 == 0) {
		sprintf (MsgText,
			 "Can't find matching dark time in %s for group %d\n",
			 nic->DARK.ref.name, ngroup);
		n_error (MsgText);
		return (status = 1);
	    }

	} else {

	    /* If the science data aren't from a MultiAccum standard exposure
	    ** time sequence, then find an exposure time that matches the
	    ** science data, or two times that bracket it */

	    for (i = 0; i < ndarks; i++) {
		 if (fabs(etime[i] - nic->exptime[ngroup-1]) <= ALLOWDIFF) {
		     nic->DarkType = MATCH;
		     nic->darkframe1 = i+1;
		     break;
		 } else if (etime[i] < nic->exptime[ngroup-1]) {
		     nic->DarkType = INTERP;
		     nic->darkframe1 = i+1;
		     etime_lower = etime[i];
		 } else if (etime[i] > nic->exptime[ngroup-1] &&
		     nic->darkframe2 == 0) {
		     nic->DarkType = INTERP;
		     nic->darkframe2 = i+1;
		     etime_upper = etime[i];
		 }
	    }
	}

	/* Reinitialize for loading next dark image */
	freeSingleNicmosGroup (dark1);
	initSingleNicmosGroup (dark1);

	/* If there's a dark image with a matching exposure time, load it. */
	if (nic->DarkType == MATCH) {
	    if (n_getRefImage (&nic->DARK, nic->darkframe1, dark1, &missing))
		return (status);
	    if (missing) {
		miss_file += 1;
		return (status = 0);
	    }

	/* Otherwise, load two images and interpolate between them */
	} else if (nic->darkframe1 != 0  &&  nic->darkframe2 != 0) {
	    if (n_getRefImage (&nic->DARK, nic->darkframe1, dark1, &missing))
		return (status);
	    if (missing) {
		miss_file += 1;
		return (status = 0);
	    }

	    if (n_getRefImage (&nic->DARK, nic->darkframe2, &dark2, &missing))
		return (status);
	    if (missing) {
		miss_file += 1;
		return (status = 0);
	    }

	    dark_interp (dark1, &dark2, etime_lower, etime_upper,
			 nic->exptime[ngroup-1]);

	    freeSingleNicmosGroup (&dark2);

	/* Otherwise, extrapolate outside the available range */
	} else if (nic->darkframe1 !=0) {
	    nic->DarkType = EXTRAP;
	    if (n_getRefImage (&nic->DARK, nic->darkframe1, dark1, &missing))
		return (status);
	    if (missing) {
		miss_file += 1;
		return (status = 0);
	    }

	    dark_extrap (dark1, etime_lower, nic->exptime[ngroup-1]);

	} else if (nic->darkframe2 !=0) {
	    nic->DarkType = EXTRAP;
	    if (n_getRefImage (&nic->DARK, nic->darkframe2, dark1, &missing))
		return (status);
	    if (missing) {
		miss_file += 1;
		return (status = 0);
	    }

	    dark_extrap (dark1, etime_upper, nic->exptime[ngroup-1]);
	}

	/* Free memory for list of exposure times after last group */
	if (nic->group == 1)
	    if (etime != NULL) free(etime);

	/* Successful return */
	return (status = 0);

}

/* DARK_INTERP: Do a linear interpolation of two dark images */

void dark_interp (SingleNicmosGroup *d1, SingleNicmosGroup *d2, 
	     double etime_d1, double etime_d2, double exptime) {

/* Arguments:
**	d1		io: dark image #1; overwritten by interpolated image
**	d2	 	 i: dark image #2
**	etime_d1   	 i: exposure time of dark image #1
**	etime_d2   	 i: exposure time of dark image #2
**	exptime		 i: desired exposure time of interpolated image
*/
	/* Local variables */
	float frac;

	/* Function definitions */
	void n_aadd  (SingleNicmosGroup *, SingleNicmosGroup *);
	void n_asub  (SingleNicmosGroup *, SingleNicmosGroup *);
	void n_amulk (SingleNicmosGroup *, float);

	/* Compute the fractional exposure time */
	frac = (exptime - etime_d1) / (etime_d2 - etime_d1);

	/* Interpolate the images */
	n_asub  (d2, d1);
	n_amulk (d2, frac);
	n_aadd  (d1, d2);
}

/* DARK_EXTRAP: Do a linear extrapolation of a dark image */

void dark_extrap (SingleNicmosGroup *d1, double etime, double exptime) {

/* Arguments:
**	d1		io: dark image; overwritten by extrapolated image
**	etime		 i: exposure time of dark image
**	exptime		 i: desired exposure time of extrapolated image
*/
	/* Local variables */
	float frac;

	/* Function definitions */
	void n_amulk (SingleNicmosGroup *, float);

	/* Compute the fractional exposure time */
	frac = exptime / etime;

	/* Extrapolate the image */
	n_amulk (d1, frac);
}

/* N_GETPHOTDATA: Read the PHOTTAB reference table and load the
** photometry parameters appropriate for the observing mode of the
** image being calibrated. The parameters will be recorded as header
** keywords in the calibrated output file.
*/

# define NCOLS		8	/* number of columns in PHOTTAB */
# define N_FLOAT_COLS	5	/* number of columns with float data type */

int n_getPhotData (NicInfo *nic, PhotData *phot) {

/* Arguments:
**	nic	io: NICMOS info structure
**	phot	 o: photometry data
*/

	/* Local variables */
	int i;				/* loop index */
	Bool found;			/* PHOTMODE found in table? */
	Bool null[N_FLOAT_COLS];	/* null column entry flags */
	int nrows;			/* number of rows in PHOTTAB */
	int col, row;			/* loop indexes */
	IRAFPointer tp;			/* PHOTTAB table pointer */
	IRAFPointer colptr[NCOLS];	/* PHOTTAB column pointers */
	char pmode[SZ_STRKWVAL+1];	/* PHOTTAB photmode entry */
	char pedigree[SZ_STRKWVAL+1];	/* Pedigree keyword value */
	float pvals[N_FLOAT_COLS];	/* photometric values */

					/* PHOTTAB column names */
	char colname[NCOLS][SZ_COLNAME+1] = {
		"PHOTMODE",
		"PHOTFLAM",
		"PHOTFNU",
		"PHOTZPT",
		"PHOTPLAM",
		"PHOTBW",
		"PEDIGREE",
		"DESCRIP"
	};

	/* Initialize values */
	pvals[0]=0; pvals[1]=0; pvals[2]=0; pvals[3]=0; pvals[4]=0;

	/* Initialize flag */
	found = False;

	/* Open the PHOTTAB reference table */
	tp = c_tbtopn (nic->PHOT.ref.name, 1, 0);
	if (c_iraferr()) {
	    n_openerr (nic->PHOT.ref.name);
	    miss_file += 1;
	    return (status = 0);	/* Don't abort yet */
	}

	/* Get pointers to columns in PHOTTAB */
	for (col = 0; col < NCOLS; col++) {
	     c_tbcfnd1 (tp, colname[col], &(colptr[col]));
	     if (c_iraferr() || colptr[col] == 0) {
		 sprintf (MsgText, "Can't find column %s in %s",
			  colname[col], nic->PHOT.ref.name);
		 n_error (MsgText);
		 c_tbtclo (tp);
		 return (status = 1);
	     }
	}

	/* Find out how many rows are in PHOTTAB */
	nrows = 0;
	nrows = c_tbpsta (tp, TBL_NROWS);
	if (nrows <= 0) {
	    sprintf (MsgText, "Invalid number of rows in %s",
		     nic->PHOT.ref.name);
	    n_error (MsgText);
	    c_tbtclo (tp);
	    return (status = 1);
	}

	/* Build the PHOTMODE string from camera and filter/grating */
	sprintf (phot->mode, "NICMOS,%1d,", nic->camera);
	strncat (phot->mode, nic->filter, strlen(nic->filter));
	strncat (phot->mode, ",DN", 3);
	for (i=0; i<strlen(phot->mode); i++)
	     phot->mode[i] = toupper(phot->mode[i]);

	/* Read each row, searching for a matching PHOTMODE */
	for (row=1; row<=nrows; row++) {

	     /* Get the PHOTMODE value in this row */
	     c_tbegtt (tp, colptr[0], row, pmode, 19);
	     if (c_iraferr()) {
		 sprintf (MsgText, "Can't read PHOTMODE in row %d in %s", row,
			  nic->PHOT.ref.name);
		 n_error (MsgText);
		 c_tbtclo (tp);
		 return (status = 1);
	     }

	     /* Does it match the PHOTMODE of the image? */
	     if (strncmp(phot->mode, pmode, strlen(pmode)) == 0) {
	        /* Yes! Found the right row */
		found = True;
		break;
	     }
	}

	if (found == True) {
	    /* Get the photometric parameters from this row */
	    c_tbrgtr (tp, &colptr[1], pvals, null, N_FLOAT_COLS, row);
	    if (c_iraferr()) {
		sprintf (MsgText, "Can't read row %d in %s", row,
			 nic->PHOT.ref.name);
		n_error (MsgText);
		c_tbtclo (tp);
		return (status = 1);
	    }

	    /* Get the PEDIGREE and DESCRIP values from this row */
	    nic->PHOT.ref.pedigree[0] = '\0';
	    c_tbegtt (tp, colptr[6], row, nic->PHOT.ref.pedigree, SZ_STRKWVAL);
	    if (c_iraferr()) {
		sprintf (MsgText, "Can't read PEDIGREE in row %d in %s", row,
			 nic->PHOT.ref.name);
		n_error (MsgText);
		c_tbtclo (tp);
		return (status = 1);
	    }

	    nic->PHOT.ref.descrip[0] = '\0';
	    c_tbegtt (tp, colptr[7], row, nic->PHOT.ref.descrip, SZ_STRKWVAL);
	    if (c_iraferr()) {
		sprintf (MsgText, "Can't read DESCRIP in row %d in %s", row,
			 nic->PHOT.ref.name);
		n_error (MsgText);
		c_tbtclo (tp);
		return (status = 1);
	    }

	    /* Is this a DUMMY reference file entry? */
	    pedigree[0] = '\0';
	    for (i=0; i<strlen(nic->PHOT.ref.pedigree); i++)
		 pedigree[i] = toupper(nic->PHOT.ref.pedigree[i]);
	    if (strncmp (pedigree, "DUMMY", 5) == 0) {
		nic->PHOT.ref.dummy = True;
		sprintf (MsgText, "PEDIGREE=DUMMY in %s", nic->PHOT.ref.name);
		n_warn (MsgText);
		c_tbtclo (tp);
		return (status = 0);
	    } else
		nic->PHOT.ref.dummy = False;

	} else {
	    /* Report failure to find a PHOTMODE match */
	    sprintf (MsgText, "Can't find PHOTMODE \"%s\" in %s", phot->mode,
		     nic->PHOT.ref.name);
	    n_error (MsgText);
	    c_tbtclo (tp);
	    return (status = 1);
	}

	/* Close the PHOTTAB reference file */
	c_tbtclo (tp);

	/* Save the photometry parameters in the NICMOS info structure */
	phot->flam = pvals[0];
	phot->fnu  = pvals[1];
	phot->zpt  = pvals[2];
	phot->plam = pvals[3];
	phot->bw   = pvals[4];

	/* Successful return */
	return (status = 0);
}

/* N_GETREFIMAGE: Load the data from a reference image. */

int n_getRefImage (CalStep *step, int extver, SingleNicmosGroup *image, 
		   Bool *missing) {

/* Arguments:
**	step	io: calibration step info structure
**	extver	 i: extension version number to load
**	image	 o: reference image data
**	missing	 o: missing file flag
*/

	/* Function definitions */
	int n_ckImgSize (SingleNicmosGroup *);
	int n_getPedigree (Hdr *, CalStep *);

	/* Initialize the image data structure */
	initSingleNicmosGroup (image);

	/* Read the data */
	if (getSingleNicmosGroup (step->ref.name, extver, image))
	    status = 1;
	if (hstio_err() || status) {
	    n_filerr (step->ref.name);
	    freeSingleNicmosGroup (image);
	    *missing = True;
	    return (status = 0);	/* Don't abort yet */
	} else {
	    *missing = False;
	}

	/* Check the size of the image */
	if (n_ckImgSize (image)) {
	    freeSingleNicmosGroup (image);
	    return (status);
	}

	/* Read the reffile PEDIGREE and DESCRIP keywords */
	if (n_getPedigree (image->globalhdr, step)) {
	    freeSingleNicmosGroup (image);
	    return (status);
	}

	/* Successful return */
	return (status = 0);
}

/* N_GETPEDIGREE: Read the PEDIGREE and DESCRIP keywords
** from a reference file header.
*/

int n_getPedigree (Hdr *refhdr, CalStep *step) {

/* Arguments:
**	refhdr	 i: reference file header structure
**	step	io: calibration step info structure
*/

	/* Local variables */
	int i;				/* loop index */
	char pedigree[SZ_STRKWVAL+1];	/* pedigree keyword value */

	/* Try to read the PEDIGREE keyword */
	if (getKeyS (refhdr, "PEDIGREE", step->ref.pedigree)) {
	    n_kwerr ("PEDIGREE", step->ref.name);
	    return (status = 1);
	}

	/* Is this a DUMMY reference file? */
	pedigree[0] = '\0';
	for (i = 0; i < strlen(step->ref.pedigree); i++)
	     pedigree[i] = toupper(step->ref.pedigree[i]);
	if (strncmp (pedigree, "DUMMY", 5) == 0) {
	    step->ref.dummy = True;
	    sprintf (MsgText, "PEDIGREE=DUMMY in %s", step->ref.name);
	    n_warn (MsgText);
	} else
	    step->ref.dummy = False;

	/* Try to read the DESCRIP keyword */
	if (getKeyS (refhdr, "DESCRIP",  step->ref.descrip))
	    status = 0;   /* not an error since it's not a required keyword */

	/* Successful return */
	return (status);
}

/* N_INITREFDATA: Initialize all reference data structures. */

void n_initRefData (SingleNicmosGroup *zoff, SingleNicmosGroup *mask,
	SingleNicmosGroup *nois, NlinData *nlin, SingleNicmosGroup *dark,
	SingleNicmosGroup *flat, PhotData *phot) {

/* Arguments:
**	zoff	o: M-ACCUM zero-read image structure
**	mask	o: Bad pixel mask image structure
**	nois	o: Detector read-noise image structure
**	nlin	o: Non-linearity coefficients data structure
**	dark	o: Dark current image structure
**	flat	o: Flat field image structure
**	phot	o: Photometry data structure
*/

	/* Initialize single-group structures */
	initSingleNicmosGroup (zoff);
	initSingleNicmosGroup (mask);
	initSingleNicmosGroup (nois);
	initSingleNicmosGroup (dark);
	initSingleNicmosGroup (flat);

	/* Initialize non-linearity data structure */
	nlin->globalhdr = NULL;
	nlin->coeff     = NULL;
	nlin->error     = NULL;
	nlin->dqual     = NULL;
	nlin->nodes     = NULL;
	nlin->zsci	= NULL;
	nlin->zerr	= NULL;

	/* Initialize photometry data structure */
	phot->mode[0] = '\0';
	phot->flam    = 0;
	phot->fnu     = 0;
	phot->zpt     = 0;
	phot->plam    = 0;
	phot->bw      = 0;

}

/* N_FREEREFDATA: Free all reference data structures and memory. */

void n_freeRefData (SingleNicmosGroup *zoff, SingleNicmosGroup *mask,
	SingleNicmosGroup *nois, NlinData *nlin, SingleNicmosGroup *dark,
	SingleNicmosGroup *flat) {

/* Arguments:
**	zoff	o: M-ACCUM zero-read image structure
**	mask	o: Bad pixel mask image structure
**	nois	o: Detector read-noise image structure
**	nlin	o: Non-linearity coefficients data structure
**	dark	o: Dark current image structure
**	flat	o: Flat field image structure
*/

	/* Local variables */
	int i;				/* loop index */

	/* Free single-group structures */
	freeSingleNicmosGroup (zoff);
	freeSingleNicmosGroup (mask);
	freeSingleNicmosGroup (nois);
	freeSingleNicmosGroup (dark);
	freeSingleNicmosGroup (flat);

	/* Free non-linearity data structure */
	if (nlin->globalhdr != NULL)
	    free (nlin->globalhdr);
	if (nlin->coeff != NULL) {
	    for (i=0; i<3; i++)
		 freeFloatHdrData (&nlin->coeff[i]);
	    nlin->coeff = NULL;
	}
	if (nlin->error != NULL) {
	    for (i=0; i<6; i++)
		 freeFloatHdrData (&nlin->error[i]);
	    nlin->error = NULL;
	}
	if (nlin->dqual != NULL) {
	    freeShortHdrData   (&nlin->dqual[0]);
	    nlin->dqual = NULL;
	}
	if (nlin->nodes != NULL) {
	    freeFloatHdrData (&nlin->nodes[0]);
	    freeFloatHdrData (&nlin->nodes[1]);
	    nlin->nodes = NULL;
	}
	if (nlin->zsci  != NULL) {
	    freeFloatHdrData (&nlin->zsci[0]);
	    nlin->zsci  = NULL;
	}
	if (nlin->zerr  != NULL) {
	    freeFloatHdrData (&nlin->zerr[0]);
	    nlin->zerr  = NULL;
	}

}

void checkKeyI (char *filename, Hdr *header, char *keyword, int scival,
	       int severity) {

	/* Local variables */
	int keyval;

	keyval = 0;

	/* Read the keyword from the ref file */
	if (getKeyI (header, keyword, &keyval)) {
	    n_kwerr (keyword, filename);
	    miss_kwd += 1;
	    status = 0;
	    return;
	}

	/* Does it match the science image value? */
	if (keyval != scival) {
	    sprintf (MsgText, "%s %s=%d does not match science data",
		     filename, keyword, keyval);
	    if (severity == FATAL) {
		mismatch += 1;
		n_error (MsgText);
	    } else
		n_warn  (MsgText);
	}
}

void checkKeyS (char *filename, Hdr *header, char *keyword, char *scival,
	       int severity) {

	/* Local variables */
	char keyval[SZ_STRKWVAL+1];

	keyval[0] = '\0';

	/* Read the keyword from the ref file */
	if (getKeyS (header, keyword, keyval)) {
	    n_kwerr (keyword, filename);
	    miss_kwd += 1;
	    status = 0;
	    return;
	}

	/* Does it match the science image value? */
	if (strncmp (keyval, scival, strlen(scival)) != 0) {
	    sprintf (MsgText, "%s %s=\"%s\" does not match science data",
		     filename, keyword, keyval);
	    if (severity == FATAL) {
		mismatch += 1;
		n_error (MsgText);
	    } else
		n_warn  (MsgText);
	}
}