# include <stdio.h>
# include <stdlib.h>	/* malloc */
# include <string.h>    /* strchr */

# include <xsynphot.h>		/* for c_phopar */

# include <hstio.h>
# include "../acs.h"
# include "../acsinfo.h"
# include "../acserr.h"

# define NPHOT      4		/* size of phot returned by c_phopar */

static void MultFilter (PhotInfo *);
static float AvgGains (ACSInfo *);

/* This routine gets the absolute flux conversion from PHOTTAB and the
   aperture throughput (filter throughput, really) from APERTAB.  These
   are multiplied together, and then the synphot routine phopar is called
   to determine the inverse sensitivity, reference magnitude (actually a
   constant), pivot wavelength, and RMS bandwidth.  These are written to
   keywords in the primary header.

   Warren Hack, 1998 June 12:
   	Initial ACS version.  
		Changes were made only to GetPhot, except for include file and 
		variable names (sts -> acs2d, StisInfo1 -> ACSInfo).
*/

int doPhot (ACSInfo *acs2d, SingleGroup *x) {

/* arguments:
ACSInfo *acs2d     i: calibration switches, etc
SingleGroup *x    io: image to be calibrated; primary header is modified
*/

	extern int status;

	PhotInfo phot;		/* QE and filter throughput, vs wavelength */
	float atodgain;
	
	float photkey[NPHOT];	/* values of photometry keywords */

	int GetPhot (ACSInfo *, PhotInfo *);
	int GetApThr (ACSInfo *, PhotInfo *);
	int PutKeyFlt (Hdr *, char *, float, char *);

	/* Get values from the photometry table. */
	/* 
		WE ARE GOING TO SELECT ON BOTH FILTERS AT ONCE...
	*/
	if (GetPhot (acs2d, &phot))
	    return (status);
	if (acs2d->photcorr == DUMMY)
	    return (status);

	/* Get values from the aperture throughput table. */
	if (acs2d->filtcorr == PERFORM) {		/* hasn't been reset? */
	    if (GetApThr (acs2d, &phot))
		return (status);
	}

	/* Combine aperture throughput with photometry table info. */
	if (acs2d->filtcorr == PERFORM)		/* still hasn't been reset? */
	    MultFilter (&phot);

	/* Compute photometry values. */
	c_phopar (phot.p_nelem, phot.p_wl, phot.p_thru, photkey);
	if (c_iraferr()) {
	    trlerror ("Error return from c_phopar");
	    return (status = ERROR_RETURN);
	}

	/* If the detector is the CCD, multiply PHOTFLAM by the gain.
	**	The gain values will be averaged to reflect the result of
	**	flat-fielding.
	*/
	atodgain = AvgGains (acs2d);
	if (acs2d->detector != MAMA_DETECTOR)
	    photkey[0] *= atodgain;

	/* Update the photometry keyword values in the SCI header. 
        Revised 10 March 1999 WJH
    */

	if (PutKeyFlt (&x->sci.hdr, "PHOTFLAM", photkey[0],
		"inverse sensitivity"))
	    return (status);

	if (PutKeyFlt (&x->sci.hdr, "PHOTZPT", photkey[1], "zero point"))
	    return (status);

	if (PutKeyFlt (&x->sci.hdr, "PHOTPLAM", photkey[2], "pivot wavelength"))
	    return (status);

	if (PutKeyFlt (&x->sci.hdr, "PHOTBW", photkey[3], "RMS bandwidth"))
	    return (status);

	free (phot.p_wl);
	free (phot.p_thru);
	if (acs2d->filtcorr == PERFORM) {
	    free (phot.f_wl);
	    free (phot.f_thru);
	}

	return (status);
}

/* This routine interpolates the aperture throughput (APERTAB) to the
   same wavelengths as the QE (PHOTTAB), and then multiplies the
   QE throughput in-place by the aperture throughput.
*/

static void MultFilter (PhotInfo *phot) {

	double wl;			/* a wavelength in QE array */
	double filt_throughput;		/* interpolated filter throughput */
	int filt_starti;		/* index to begin search in interp1d */
	int i;
	double interp1d (double, double *, double *, int, int *);

	filt_starti = 1;			/* initial value */

	for (i = 0;  i < phot->p_nelem;  i++) {

	    wl = phot->p_wl[i];

	    filt_throughput = interp1d (wl, phot->f_wl, phot->f_thru,
			phot->f_nelem, &filt_starti);

	    phot->p_thru[i] *= filt_throughput;
	}
}

static float AvgGains (ACSInfo *acs2d) {

	float gain;
    float gn[NAMPS];
    float rn2[NAMPS];   /* We only need it to call get_nsegn */
	int numamps;
	int i;
	char ccdamp[NAMPS+1];
	
	void get_nsegn (int, int, int, int, float *, float *, float *, float *);
    
    /* Initialize internal variables */
	numamps = 0;
	gain = 0.0;
	ccdamp[0] = '\0';
	for (i = 0; i < NAMPS; i++){
        gn[i] = 0.;
        rn2[i] = 0.;
    }
    
    /* Get gain values appropriate for observation */
    get_nsegn (acs2d->detector, acs2d->chip, acs2d->ampx, acs2d->ampy, acs2d->atodgain, acs2d->readnoise, gn, rn2);

	for (i=0; i < NAMPS; i++) {
		/* If one of the amps used is found, ... */
		if (gn[i] > 0.) {
			/* try to use the atodgain value from that amp */
            numamps++;
            gain += gn[i];
		}
	}
	
	if (numamps > 0) 
		return(gain / numamps);
	else 
		return (0.0);
}