# include <stdio.h>

# include <hstio.h>
# include "../stis.h"
# include "calstis7.h"
# include "../stisdq.h"

# define  HST_AREA  45238.93416		/* cm**2, including obscured areas */
# define  H_PLANCK  (6.6260755e-27)	/* Planck's constant, erg * sec */
# define  C_LIGHT   29979245800.	/* cm / sec */
# define  CM_PER_ANGSTROM  (1.e-8)	/* cm / angstrom */

/* This routine corrects the output SCI and ERR data (in-place) for
   telescope throughput and detector sensitivity.  The data before
   correction are expected to be in counts; after correction the BUNIT
   keyword in the science and error extension headers will be set to
   "erg /s /cm**2 /angstrom /arcsec**2".

   Counts are converted to specific intensity by multiplying by:

				h * c * atodgain
	--------------------------------------------------------------
	exptime * R(lambda) * lambda * hst_area * disp * scale * width

   where:

	h = Planck's constant (erg * sec); c = speed of light (cm / sec)
	atodgain = analog to digital gain (electrons / dn)
	exptime = exposure time (sec)
	lambda = wavelength (cm)
	R(lambda) = instrumental response (e.g. output from synphot.calcband)
	hst_area = area of HST (cm**2), _including_ obscured regions
	disp = dispersion (angstrom / pixel)
	scale = image scale along slit (arcsec / pixel)
	width = slit width (arcsec)

   The R(lambda) that we got from the phottab is for an "infinite"
   extraction box height, which is what we need for a diffuse source.

   The DIFF2PT keyword is assigned the value (scale * width * pcorr) / Ta,
   where Ta is the aperture throughput for a point source, scale and width
   are as above, and pcorr is the factor read from the PCTAB, the ratio of
   R(lambda) for an infinite extraction box height divided by R(lambda)
   for the default extraction box.  DIFF2PT will be written to the SCI
   extension header.  Since Ta and pcorr are functions of wavelength and
   we only write one value to the header, we average the throughput over
   the wavelength range of the current extension when computing DIFF2PT.
### (Actually, this range is too large because it covers the full output
   image in the dispersion direction rather than just the pixels mapped
   from points within the input image.)
   For a point source, to convert the specific intensity to flux density
   Flambda, sum the SCI image values along the length of the slit and
   multiply the sum by DIFF2PT.

   Phil Hodge, 1997 Sept 5:
	Include atodgain as a factor.
	plate_scale is now an array of two elements, to allow for unequal
	binning in the two axes.

   Phil Hodge, 1997 Nov 13:
	Assume dispaxis = 1, and remove it from the calling sequence.
	Correct R(lambda) for default extraction box height using pcorr
	in the phot struct, and take account of this in the DIFF2PT factor.

   Phil Hodge, 1998 Jan 26:
	Remove the pcorr factor from the calibration for diffuse source;
	that is, assume R(lambda) is for an infinite extraction box height.
	(We've been assuming all along that R(lambda) has been corrected
	for the aperture throughput.)

   Phil Hodge, 2001 Feb 20:
	Check for zero or negative wavelength (needed for prism).
*/

int AbsFlux (SingleGroup *out, CoordInfo *coord_o,
		PhotInfo *phot, ApInfo *slit,
		double *plate_scale, double atodgain,
		double exptime, double hfactor) {

/* arguments:
SingleGroup *out       io: output data
CoordInfo *coord_o     i: coordinate parameters
PhotInfo *phot         i: photometry info
ApInfo *slit           i: slit info (for throughput)
double plate_scale[2]  i: arcsec / pixel for each axis
double atodgain        i: electrons / adu for CCD
double exptime         i: exposure time in seconds
double hfactor         i: divide to convert back to observed wavelength
*/

	int status;

	double wl;		/* heliocentric wavelength at a pixel */
	double response;	/* instrumental response at a wavelength */
	double pct_factor;	/* correction to infinite extr box height */
	double throughput;	/* factor to correct for slit throughput */
	double sum;		/* for getting average throughput */
	double photfactor;	/* to get inverse sensitivity */
	float correction;	/* combined correction factor */
	float diff_pt;		/* conversion from diffuse source to point */
	float cont_eml;		/* conversion from continuum to emission line */
	int i, j;
	int abs_starti;		/* index to begin search in interp1d */
	int pct_starti;
	int thr_starti;

	double interp1d (double, double *, double *, int, int *);
	int Put_KeyS (Hdr *, char *, char *, char *);
	int Put_KeyF (Hdr *, char *, float, char *);
	int Put_KeyD (Hdr *, char *, double, char *);

	abs_starti = 1;				/* initial values */
	pct_starti = 1;
	thr_starti = 1;
	sum = 0.;

	/* incomplete; will also divide by average throughput */
	diff_pt = coord_o->cdelt[1] * slit->width[0];

	cont_eml = coord_o->cdelt[0] * slit->width[0] /
			plate_scale[0];

	/* The photometry table gives the instrumental response for
	   a point source.  This factor is part of the conversion
	   from count rate to specific intensity.
	*/
	photfactor = H_PLANCK * C_LIGHT * atodgain / (HST_AREA * exptime *
		coord_o->cdelt[0] * coord_o->cdelt[1] * slit->width[0]);

	/* for each pixel in the dispersion direction ... */
	for (i = 0;  i < out->sci.data.nx;  i++) {

	    /* convert back to observed wavelength */
	    wl = (((double)i - coord_o->crpix[0]) * coord_o->cdelt[0] +
			coord_o->crval[0]) / hfactor;	/* Angstroms */
	    if (wl <= 0.) {
		for (j = 0;  j < out->dq.data.ny;  j++) {
		    DQSetPix (out->dq.data, i, j,
				DQPix(out->dq.data,i,j) | CALIBDEFECT);
		}
		continue;
	    }

	    /* interpolate */
	    response = interp1d (wl, phot->wl, phot->thru,
			phot->nelem, &abs_starti);
	    pct_factor = interp1d (wl, phot->wl, phot->pcorr,
			phot->nelem, &pct_starti);
	    throughput = interp1d (wl, slit->wl, slit->thr,
			slit->nelem, &thr_starti);
	    /* We need the average throughput for diff2pt. */
	    sum += (throughput / pct_factor);
	    if (response <= 0.)
		continue;

	    wl *= CM_PER_ANGSTROM;	/* convert from Angstroms to cm */
	    correction = (float) (photfactor / (response * wl));

	    /* for each pixel along the slit ... */
	    if (correction <= 0.) {
		for (j = 0;  j < out->dq.data.ny;  j++) {
		    DQSetPix (out->dq.data, i, j,
				DQPix(out->dq.data,i,j) | CALIBDEFECT);
		}
	    } else {
		for (j = 0;  j < out->sci.data.ny;  j++) {
		    Pix (out->sci.data, i, j) *= correction;
		    Pix (out->err.data, i, j) *= correction;
		}
	    }
	}

	/* average throughput */
	throughput = sum / (double)(out->sci.data.nx);

	/* Update BUNIT in the science and error extension headers. */
	if (status = Put_KeyS (&out->sci.hdr, "BUNIT", phot->bunit,
		"units for flux-calibrated data"))
	    return (status);
	if (status = Put_KeyS (&out->err.hdr, "BUNIT", phot->bunit,
		"units for flux-calibrated data"))
	    return (status);

	if (throughput > 0.)
	    diff_pt /= throughput;
	else
	    diff_pt = -9999.;
	if (status = Put_KeyF (&out->sci.hdr, "DIFF2PT", diff_pt,
		"diffuse to point source conversion factor"))
	    return (status);

	if (status = Put_KeyF (&out->sci.hdr, "CONT2EML", cont_eml,
		"continuum to emission line conversion factor"))
	    return (status);
	if (status = Put_KeyF (&out->sci.hdr, "SCALE_A1", plate_scale[0],
		"arcsec / pixel in dispersion direction"))
	    return (status);
	if (status = Put_KeyD (&out->sci.hdr, "OMEGAPIX",
		plate_scale[0] * plate_scale[1], "pixel area, arcsec^2"))
	    return (status);

	return (0);
}