# include <stdio.h>

# include <xtables.h>
# include "../stis.h"
# include "calstis6.h"
# include "../stisdq.h"

/* Generates the FLUX array from the NET array, and also modifies the ERROR 
   array in-place for telescope throughput and detector sensitivity. The 
   data before correction are expected to be in c/s; after correction its
   units will be "erg /s /cm**2 /angstrom".

   Counts are converted to specific intensity by multiplying by:

                        h * c * Ta
            ------------------------------------
            R(lambda) * lambda * hst_area * disp

   where:

	h = Planck's constant (erg * sec); c = speed of light (cm / 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)
	Ta = aperture throughput for a point source

   Ta is included in R(lambda), and we're multiplying by it in order to
   remove it for the diffuse-source case.

   For now, the extraction box height correction (from PCTAB) assumes that
   the default extraction box (from XTRACTAB) was used. This works in
   pipeline mode, but may fail when calstis6 is run from the command line
   with a non-default extraction box height. 




   Revision history:
   ----------------
   26 Feb 97  -  Implemented, with some code adapted from a similar 
                 routine in calstis7 (I.Busko)
   10 Apr 97  -  Changes after code review (IB):
                 - C_LIGHT in scientific notation
   03 Sep 97  -  Add ATODGAIN correction (IB)
   28 Jan 98  -  Add extraction box height correction (IB)
   09 Apr 99  -  Check for zeroed throughput - OPR 38749 (IB)
   17 Dec 99  -  Flux normalization factor in opt. extr. rejection (IB) 
   17 Aug 00  -  Fix dispersion computation in hires data (IB)
   26 Feb 01  -  Got rid of EvalDisp (IB)
   06 Mar 01  -  Correct wavelengths back to observed values (IB)
*/

void AbsFlux6 (StisInfo6 *sts, DispRelation *disp, RowContents *out, 
               PhotInfo *phot, 
               PhotInfo *photc, ApInfo *slit, double hfactor, 
               double atodgain, int optimal, double *pfactor, double ltm) {

/* arguments:
StisInfo6 *sts         i: calibration switches and info
DispRelation *disp     i:  dispersion relation
RowContents *out       io: output data
PhotInfo *phot         i: photometry info
PhotInfo *photc        i: for extraction box height correction
ApInfo *slit           i: slit info (for throughput)
double hfactor         i: divide to convert back to observed wavelength
double atodgain        i: CCD ATODGAIN value
int optimal            i: is optimal extraction ?
double *pfactor;       i: flux normalization factor from opt. extr. profile
double ltm;	       i: LTM1_1 value 
*/

	double wl;		/* observed wavelength at a pixel */
	double response;	/* instrumental response at a wavelength */
	double pct_factor;	/* correction to infinite extr box height */
	double pct_factor2;	/* correction for current extr box height */
	double throughput;	/* factor to correct for slit throughput */
	double photfactor;	/* to get inverse sensitivity */
	double dispersion;	/* dispersion in A/pix */
	double dummy;
	float correction;	/* combined correction factor */
	float ecorrection;	/* combined correction factor for the error */
	int i;
	int abs_starti;		/* index to begin search in interp1d */
	int pct_starti;
	int thr_starti;
	double hf;

	/* void EvalDisp6 (DispRelation *, int, double, double *, double *); */
	double interp1d (double, double *, double *, int, int *);

	abs_starti = 1;				/* initial values */
	pct_starti = 1;
	thr_starti = 1;

	photfactor = H_PLANCK * C_LIGHT / HST_AREA;

	if (sts->heliocorr == PERFORM)
	    hf = sts->hfactor;
	else
	    hf = 1.0;

	/* Loop over flux array. */
	for (i = 0;  i < out->npts;  i++) {

	    /* Convert back to observed wavelenght. */
	    wl = out->wave[i] / hfactor;

	    /* Interpolate in the reference tables. */
	    response    = interp1d (wl, phot->wl, phot->thru, phot->nelem, 
                                    &abs_starti);
	    pct_factor  = interp1d (wl, phot->wl, phot->pcorr, phot->nelem, 
                                    &pct_starti);
	    pct_factor2 = interp1d (wl, photc->wl, photc->pcorr, photc->nelem, 
                                    &pct_starti);
	    throughput  = interp1d (wl, slit->wl, slit->thr, slit->nelem, 
                                    &thr_starti);

	    /* Check for zeroed throughput (OPR 38749) */
	    if (throughput <= 0.0) {
	        out->flux[i]  = 0.0F;
	        out->error[i] = 0.0F;
	        out->dq[i]   |= CALIBDEFECT;
	        continue;
	    }

	    /* Compute dispersion. */
	    /* EvalDisp6 (disp, out->sporder, wl, &dummy, &dispersion); */
	    if (i > 0 && i < out->npts - 1)
	        dispersion = (out->wave[i+1] - out->wave[i-1]) / 2.0;
            else if (i == 0)
	        dispersion = out->wave[1] - out->wave[0];
	    else
	        dispersion = out->wave[out->npts-1] - out->wave[out->npts-2];
	    if (dispersion <= 0.0) {
	        out->flux[i]  = 0.0F;
	        out->error[i] = 0.0F;
	        out->dq[i]   |= CALIBDEFECT;
	        continue;
	    }
	    dispersion /= hf;

	    /* Compute flux. */
	    if (response <= 0.0) {
	        out->flux[i]  = 0.0F;
	        out->error[i] = 0.0F;
	        out->dq[i]   |= CALIBDEFECT;
	    } else {
	        correction = (float) (photfactor / (response * throughput * 
                                      wl * dispersion * CM_PER_ANGSTROM));
                correction *= atodgain;  /* added 9/3/97 (IB) */
                correction *= pct_factor / pct_factor2;
	        ecorrection = correction;
	        if (optimal)
	            correction *= pfactor[i];
	        if (ecorrection <= 0.) {
	            out->flux[i]  = 0.0F;
	            out->error[i] = 0.0F;
	            out->dq[i]   |= CALIBDEFECT;
	        } else {
	            out->flux[i] = out->net[i] * correction;
	            out->error[i] = out->error[i] * ecorrection;
	        }
	    }
	}
}