path: root/src/libcf/cf_optimal_extraction.c

/*****************************************************************************
 *              Johns Hopkins University 
 *              Center For Astrophysical Sciences
 *              FUSE
 *****************************************************************************
 *
 * Synopsis:    cf_optimal_extraction(fitsfile *header, int optimal, 
 *              int aperture, float *weight, float *y, unsigned char *channel,
 *              float *lambda, long ngood, long *good_index,
 *		float *barray, float *bpmask, int pny, float pycent,
 *		float *parray, long nout, float *wave_out, float **flux_out,
 *		float **sigma_out, long **counts_out, float **weights_out, 
 *		float **bkgd_out, short **bpix_out)
 *
 * Description: Perform standard or optimal extraction of target spectrum.
 *
 * Note:	We use a modified version of the optimal-extraction
 *		algorithm described by Keith Horne (PASP, 98, 609, 1986).
 *
 * Arguments:   fitsfile  *header       Pointer to IDF FITS file header
 *		int	  optimal	TRUE if optimal extraction requested.
 *		int	  aperture	Target aperture (values 1-8)
 *		float	  *weight       Scale factor for each photon
 *		float	  *y            Final y position of each photon
 *		unsigned char *channel	Channel ID of each photon
 *		float 	  *lambda       Wavelength for each photon
 *		long	  ngood		Length of array good_index 
 *		long	  *good_index	Indices of screened photon events
 *		float     *barray       2-D background array
 *                                       (flux-calibrated and binned)
 *              float     *bpmask       2-D array containing bad pixel locations
 *                                      Same size and location as the background array
 *	        int       pny           Size in Y of weights array
 *	        float     pycent        Y centroid of weights array
 *	        float     parray        2-D weights array (binned)
 *	        long      nout          Length of output arrays
 *	        float     *wave_out	Output wavelength array
 *	        float     **flux_out    Output flux array
 *	        float     **sigma_out   Output error array
 *	        long      **counts_out  Output counts array
 *	        float     **weights_out Output weights array
 *	        float     **bkgd_out 	Output background array
 *              short     **bpix_out    Output bad pixel spectrum
 *
 * Returns:     0 on success
 *
 * History:     01/28/03   1.1   wvd    Initial coding
 *              02/28/03   1.2   peb    Changed flux_out, var_out, counts_out
 *                                      and weight_out argument variables to
 *                                      pointer-to-pointers, so the arrays can
 *                                      be returned to the calling routine.
 *                                      Also tidied the code using lint.
 *		03/11/03   1.3   wvd    Update documentation regarding use of
 *					unsigned char
 *		03/12/03   1.4   wvd    Change weights_out to sum of weights
 *		03/19/03   1.5   wvd    Divide flux_out and var_out by WPC
 *              04/04/03   1.6   rdr    - correct bugs in optimum extraction
 *                                      - adjust criterion for redetermining
 *                                        y centroid
 *		05/07/03   1.7   wvd	Change ergcm2s to ergcm2 throughout. 
 *					Divide by EXPTIME.  Use cf_verbose.
 *		05/22/03   1.8   wvd	If EXPTIME = 0, set flux and error to 0.
 *              05/28/03   1.9   rdr    - For HIST data, spread counts
 *					  along the y dimension of
 *                                        the data and variance arrays
 *                                      - Incorporate the bad pixel mask into
 *                                        optimal extraction algorithm
 *                                      - Return bad-pixel spectrum
 *              09/30/03   1.10  wvd    Negative values of YCENT are 
 *					meaningless.  Use S/N of data to 
 *					determine significance of YCENT.
 *					Exclude airglow lines from S/N.
 *              10/08/03   1.11  wvd    Change counts_out array to type long.
 *					Add 0.5 to it before rounding.
 *					Limit Y coordinates of probability
 *					array before comparing with bkgd.
 *              10/31/03   1.12  wvd    Read centroid quality flags from IDF
 *					header.  Do not attempt to calculate
 *					spectral centroid.
 *              11/03/03   1.13  wvd    Modify calculation of variance estimate.
 *					Modify scheme for smoothing HIST data.
 *					Iterate only to 20% accuracy.
 *              12/05/03   1.14  wvd    Check value of SPECBINY before
 *					smoothing HIST data in the Y dimension.
 *              03/16/04   1.15  wvd    In HIST mode, assume counts are
 *					distributed across one X pixel.
 *					Delete wave_out -- not used.
 *					After 50 iterations, abandon optimal
 *					extraction and use boxcar extraction
 *					instead.
 *					Don't test size of barray,
 *					as it is now the same as parray.
 *					Change pycent from int to float.
 *              04/06/04   1.16  bjg    Include string.h
 *                                      Remove unused variables
 *              04/21/04   1.17  wvd	Set dispapix = fabs(dispapix).
 *					Iterate to 5% accuracy.
 *					In boxcar extraction, use parray to
 *					set extraction limits.
 *					For TTAG data, compute counts_out
 *					directly from photon list.
 *              04/26/04   1.18  wvd	Perform all calculations with
 *					weights array, rather than ergcm2.
 *					Return counts spectrum.
 *					Convert from variance to sigma here,
 *					rather than in cf_write_spectra.
 *					Iterate to 10% accuracy.
 *              05/17/04   1.19  wvd	If weights_out[j] = 0, set variance
 *					= bkgd.  Iterate optimal extraction
 *					to accuracy of 0.01 counts. In boxcar
 *					extraction, reduce parray limit to
 *					1E-4 to better match fluxes from 
 *					optimal extraction for bright stars.
 *              08/13/04   1.20  wvd	Use QUALITY array to set X limits 
 *					of extraction region.  Insure that
 *					sigma_out array is non-zero.
 *              11/22/05   1.21  wvd	Scale variance estimate by bad-pixel
 *					mask to better approximate observed
 *					counts.
 *              06/12/06   1.22  wvd	Call cf_verbose rather than
 *					cf_if_warning.
 *
 ****************************************************************************/

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

#include "calfuse.h"

static char 	CF_PRGM_ID[] = "cf_optimal_extraction";
static char 	CF_VER_NUM[] = "1.22";

int
cf_optimal_extraction(fitsfile *header, int optimal, int aperture, 
	float *weight, float *y, unsigned char *channel,
	float *lambda, long ngood, long *good_index,
	float *barray, float *bpmask, 
	int pny, float pycent, float *parray, long nout, float *wave_out,
	float **flux_out, float **sigma_out, long **counts_out,
	float **weights_out, float **bkgd_out, short **bpix_out)

{
    char	instmode[FLEN_VALUE], yquality_str[FLEN_KEYWORD], 
		ycent_str[FLEN_KEYWORD], yquality[FLEN_VALUE];
    int   	status=0;
    int		nloop;			/* Iterations of optimal
						extraction loop */
    int		ycent_uncertain = FALSE;/* TRUE if centroid quality
						less than HIGH */

    long	i, 			/* index through photon list */
		ii,			/* ii = good_index[i] */
		j, 			/* X index through 2-D arrays */
		jmin, jmax,		/* limits of output spectrum */
		k; 			/* Y index through 2-D arrays */
    
		 
    float	exptime,		/* exposure time */
		tdead,			/* mean dead-time correction */
		w0,			/* minimum value of output wave array */
		wpc,			/* wavelength increment per output bin */
		*flux_old,		/* flux values from previous iteration */
		*flux_box,		/* flux values from boxcar extraction */
		*var_out,		/* 1-D variance array */
		*var_box,		/* variance values from boxcar extraction */

		ycent;			/* Y centroid of target spectrum */

    double      numerator, denominator, /* Used in optimal-extraction calculations */
		var_num, weights_num,
		bkgd_num,
		*data,			/* 2-D data array */
		*variance,		/* 2-D variance array */
		*w;			/* weights array for variance calculation */

    int		l, lmin, lmax,		/* These variables are	*/
    		specbiny, specbiny2;	/* are used to smooth	*/
    float 	psum, pval;		/* HIST data in the	*/
    double	vval;			/* Y dimension.		*/

    long	jj, jjj, m, m_min, n;	/* These variables are used to smooth */
    double	edge[2], frac[2], x; 	/* HIST data in the X dimension.      */
    double	dispapix;		/* Mean dispersion per pixel */
    char	wave_file[FLEN_VALUE];  /* Name of WAVE_CAL file */
    int		hdunum;			/* HDU number of WAVE_CAL file */
    fitsfile	*wavefits;		/* Pointer to WAVE_CAL file */

    cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
    cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");

    cf_verbose(3, "Entering cf_optimal_extraction");

    /* 
     *  Check type of observation (HIST or TTAG) and binning factor.
     */
    FITS_movabs_hdu(header, 1, NULL, &status) ;
    FITS_read_key(header, TSTRING, "INSTMODE", instmode, NULL, &status);
    FITS_read_key(header, TINT, "SPECBINY", &specbiny, NULL, &status);
    specbiny2 = specbiny/2 ;
    if (specbiny2 < 1) specbiny2 = 1 ;
    cf_verbose(3, "Aperture = %d, instrument mode = %s, binning in Y = %d",
	aperture, instmode, specbiny) ;

    /*
     *  Read EXPTIME and TOT_DEAD from file header.
     */
    FITS_read_key(header, TFLOAT, "EXPTIME", &exptime, NULL, &status);
    FITS_read_key(header, TFLOAT, "TOT_DEAD", &tdead, NULL, &status);
    if (tdead < 1.0) tdead = 1.0;
    cf_verbose (3, "EXPTIME = %.0f, TOT_DEAD=%g", exptime, tdead);

    /*
     *  Read wavelength parameters W0 and WPC from file header.
     */
    FITS_read_key(header, TFLOAT, "WPC", &wpc, NULL, &status);
    FITS_read_key(header, TFLOAT, "W0", &w0, NULL, &status);
    cf_verbose (3, "Wavelength info: w0=%g, dw=%g", w0, wpc) ;

    /*
     *  Read spectral centroid YCENT and YQUAL from file header.
     */
    sprintf(ycent_str, "YCENT%1d", aperture);
    FITS_read_key(header, TFLOAT, ycent_str, &ycent, NULL, &status);
    sprintf(yquality_str, "YQUAL%1d", aperture);
    FITS_read_key(header, TSTRING, yquality_str, yquality, NULL, &status);
    if (*yquality != 'H') ycent_uncertain = TRUE;
    cf_verbose (3, "%s = %g, quality = %s", ycent_str, ycent, yquality);

    /*
     * If HIST data, read DISPAPIX from header of WAVE_CAL file.
     */
    if (!strncmp(instmode,"H",1)) {
	hdunum = aperture + 1;
	FITS_read_key(header, TSTRING, "WAVE_CAL", wave_file, NULL, &status);
	FITS_open_file(&wavefits, cf_cal_file(wave_file), READONLY, &status);
	FITS_movabs_hdu(wavefits, hdunum, NULL, &status);
	FITS_read_key(wavefits, TDOUBLE, "DISPAPIX", &dispapix, NULL, &status);
	dispapix = fabs(dispapix);
	FITS_close_file(wavefits, &status);
    }

    /*
     *  Allocate space for all arrays.
     */
     data = (double *) cf_calloc(pny*nout, sizeof(double));
     variance = (double *) cf_calloc(pny*nout, sizeof(double));
     w = (double *) cf_calloc(nout, sizeof(double));

     flux_old = (float *) cf_calloc(nout, sizeof(float));
     flux_box = (float *) cf_calloc(nout, sizeof(float));
     var_box  = (float *) cf_calloc(nout, sizeof(float));

     *flux_out = (float *) cf_calloc(nout, sizeof(float));
     var_out = (float *) cf_calloc(nout, sizeof(float));
     *counts_out = (long *) cf_calloc(nout, sizeof(long));
     *weights_out = (float *) cf_calloc(nout, sizeof(float));
     *bkgd_out = (float *) cf_calloc(nout, sizeof(float));
     *bpix_out = (short *) cf_calloc(nout, sizeof(short)) ;

    /* 
     *  The quality array (bpix_out) is the product of the probability
     *  array and the bad-pixel mask.  We use it to set the limits of
     *  the extraction region.
     */
     cf_verbose(3, "Generating bad-pixel spectrum.") ;
     for (j=0; j<nout; j++) {
         float bpix=0. ;
         for (k=0; k<pny; k++)
	   bpix += (parray[k*nout + j] * bpmask[k*nout + j]);
         (*bpix_out)[j] = (short) (100. * bpix + 0.5) ;
     }
     jmin = 0L; jmax = nout;
     for (j=0; j<nout; j++) {
	if ((*bpix_out)[j] > 0) {
	    jmin = j;
	    break;
	}
     }
     for (j = nout-1; j >= 0; j--) {
	if ((*bpix_out)[j] > 0) {
	    jmax = j+1;
	    break;
	}
     }

    /*
     *  Construct data and variance arrays from photon lists.
     *  Arrays are pny pixels in Y (to match probability array) 
     *  and nout pixels in X (to match output wavelength array).  
     *  We must thus shift individual photons by PYCENT - YCENT
     *  pixels in Y.
     */

     cf_verbose(3, "Filling the data array") ;
     for (i = 0; i < ngood; i++) {
	ii = good_index[i];

	if (channel[ii] != aperture) continue;

	j = (long) ((lambda[ii] - w0)/wpc + 0.5);
	if (j < jmin || j >= jmax) continue;

	k = (long) (y[ii] - ycent + pycent + 0.5);
	if (k < 0 || k >= pny) continue;

        if (!strncmp(instmode,"T",1)) {
	   /* TTAG data */
	    (*counts_out)[j]++;
	    data[k*nout + j] += weight[ii];
	    variance[k*nout + j] += weight[ii] * weight[ii];
	}
	else { 
	   /* If HIST data are already smoothed in Y, it's easy: */
	   if (specbiny == 1) {
	       data[k*nout + j] += weight[ii];
	       variance[k*nout + j] += weight[ii] * tdead;
	   }
	   /* If HIST data are binned, smooth in both X and Y dimensions. */
	    else {
		m_min = 0;
		edge[0] = ((double) j - 0.5) * wpc + w0;
		edge[1] = ((double) j + 0.5) * wpc + w0;
		if ((x = (lambda[ii] - edge[0]) / dispapix) < 0.5) {
		   frac[0] = 0.5 - x;
		   frac[1] = 0.5 + x;
		   jj = j - 1;
		   n = 2;
		   if (jj == -1) m_min = 1;
		}
		else if ((x = (edge[1] - lambda[ii]) / dispapix) < 0.5) {
		   frac[0] = 0.5 + x;
		   frac[1] = 0.5 - x;
		   jj = j;
		   n = 2;
		   if (jj + 1 == nout) n = 1;
		}
		else {
		   frac[0] = 1.0;
		   jj = j;
		   n = 1;
		}

		psum = 0.;
		vval = weight[ii] * tdead;
		lmin = k - specbiny2;
		lmax = lmin + specbiny;
		if (lmin < 0) lmin = 0;
		if (lmax > pny) lmax = pny;
		pval = 1. / (lmax - lmin);
		for (l = lmin; l < lmax; l++)
		    psum += parray[l*nout + j];
		for (l = lmin; l < lmax; l++) {
		    if (psum > 0) pval = parray[l*nout + j] / psum;
		    for (m = m_min; m < n; m++) {
			jjj = jj + m;
			data[l*nout + jjj] += pval * frac[m] * weight[ii];
			variance[l*nout + jjj] += pval * frac[m] * vval;
		    }
		}
	    }
	}
     }

    /*
     *  Boxcar extraction.  Use parray to set extraction limits in Y.
     */
     cf_verbose(3, "Doing the boxcar extraction.") ;
     nloop = 0;
     for (j = jmin; j < jmax; j++) {
	for (k = 0; k < pny; k++) {
	   if (parray[k*nout + j] > 1E-4) {
	         (*bkgd_out)[j] += barray[k*nout + j];
	             var_out[j] += variance[k*nout + j];
	      (*weights_out)[j] += data[k*nout + j];
	   }
	}
	(*flux_out)[j] = (*weights_out)[j] - (*bkgd_out)[j];
     }

    /*
     *  Has the user requested optimal extraction?
     */
     cf_verbose(3,"Optimal keyword = %d ", optimal) ;
     if (optimal) {

       /*
	*  If Y centroid is uncertain, abandon optimal extraction.
	*/
	if (ycent_uncertain)
	   cf_verbose(1, "%s uncertain.  Cannot perform optimal extraction.",
		ycent_str);

	else {		/* Proceed with optimal extraction */

 	    cf_verbose(3, "Attempting optimal extraction");

	/*
	 *  Save boxcar versions of flux and variance arrays.
	 */
	for (j = jmin; j < jmax; j++) {
	    flux_box[j] = (*flux_out)[j];
	    var_box[j] = var_out[j];
	}

	/*
	 *  Need array w[j] to scale the variance estimate.
	 *  w[j] is mean scale factor from counts to weights at each wavelength.
	 *  Must calculate for TTAG data; can use TOT_DEAD for HIST data.
	 */
	if (!strncmp(instmode,"T",1)) for (j = jmin; j < jmax; j++) {
	    if ((*counts_out)[j] > 0) {
		for (k = 0; k < pny; k++) w[j] += data[k*nout + j];
		w[j] /= (*counts_out)[j];
	    }
	    else w[j] = tdead;
	}
	else for (j = jmin; j < jmax; j++) w[j] = tdead;

	/*
	 *  Begin big loop.
	 */
	do {
	    /*
	     *  Revise variance estimate
	     *  We can omit w[j] here, as it appears in both the numerator
	     *  and denominator of the flux estimate (below).
	     *  05/17/02 wvd:  If weights_out[j] = 0, set variance = bkgd.
	     *  11/22/05 wvd:  Scale variance estimate by bpmask.
	     */
	     for (j = jmin; j < jmax; j++) {
		if ((*weights_out)[j] > 0)
		   for (k = 0; k < pny; k++)
			variance[k*nout + j] = bpmask[k*nout + j] *
			fabs(parray[k*nout + j] * (*flux_out)[j] + barray[k*nout + j]);
		else
		   for (k = 0; k < pny; k++)
			variance[k*nout + j] = bpmask[k*nout + j] * barray[k*nout + j];
	     }

	    /*
	     *  Optimal extraction
	     */
	     nloop++;
	     for (j = jmin; j < jmax; j++) {
		numerator = denominator = 0.;
		for (k = 0; k < pny; k++) if (variance[k*nout+j] > 0) {
		    numerator += bpmask[k*nout+j] * parray[k*nout + j] * 
                      (data[k*nout + j] - barray[k*nout + j]) / variance[k*nout + j];
		    denominator += bpmask[k*nout+j] * parray[k*nout + j] * 
                       parray[k*nout + j] / variance[k*nout + j];
		}
		flux_old[j] = (*flux_out)[j];
		if (denominator > 0) (*flux_out)[j] = numerator / denominator;
                else (*flux_out)[j] = flux_old[j] = 0. ;
	     }

	    /*
	     *  Repeat loop if any element of flux_out has changed by > 0.01 counts.
	     */
	     for (j = jmin; j < jmax; j++) {
		if (fabs((*flux_out)[j]-flux_old[j]) > 0.01) {
		    cf_verbose(3, "%d\t%d\t%g", nloop, j, (*flux_out)[j]);
		    break;
		}
	     }

	} while (j < jmax && nloop < 50);	/* End of do loop */

	/*
	 *  If nloop < 50, calculate var_out, weights_out, and bkgd_out.
	 *  Note that we must scale var_out by w[j] here.
	 */
	if (nloop < 50) {

	    memset (*bkgd_out, 0, nout*sizeof(float));
	    memset (*weights_out, 0, nout*sizeof(float));
	    memset (var_out, 0, nout*sizeof(float));

	    for (j = jmin; j < jmax; j++) {
		bkgd_num = weights_num = var_num = denominator = 0.;
		for (k = 0; k < pny; k++) if (variance[k*nout+j] > 0) {
		   bkgd_num += bpmask[k*nout+j] * parray[k*nout + j] *
                      barray[k*nout + j] / variance[k*nout + j];
		   weights_num += bpmask[k*nout+j] * parray[k*nout + j] *
                      data[k*nout + j] / variance[k*nout + j];
		   var_num += bpmask[k*nout+j] * parray[k*nout + j];
		   denominator += bpmask[k*nout+j] * parray[k*nout + j] * 
                        parray[k*nout + j] / variance[k*nout + j];
		}
		if (denominator > 0) {
		    (*bkgd_out)[j] = bkgd_num / denominator;
		    (*weights_out)[j] = weights_num / denominator;
		    var_out[j] = var_num / denominator * w[j];
		}
		else {
		    (*bkgd_out)[j] = (*weights_out)[j] = var_out[j] = 0.;
		}
	    }
	}
	/*
	 *  If nloop = 50, use boxcar extraction instead.
	 */
	else {
	    cf_verbose(1, "Optimal extraction failed to converge for aperture %d."
		"  Using boxcar extraction.", aperture);
	    nloop = 0;
	    for (j = jmin; j < jmax; j++) {
		(*flux_out)[j] = flux_box[j];
		var_out[j] = var_box[j];
            }
	}

	}			/* End of if/else test of YCENT */
     }				/* End of optimal extraction block */

    /*
     *  Write number of iterations to file header.
     */
     FITS_update_key(header, TINT, "OPT_EXTR", &nloop, NULL, &status);
     cf_verbose(2, "Number of iterations in optimal extraction = %d", nloop) ;

    /* Convert error array from variance to sigma. */
    for (j = jmin; j < jmax; j ++) {
	if (var_out[j] > 0.) var_out[j] = sqrt(var_out[j]);
	else if ((*bkgd_out)[j] > 0.) var_out[j] = sqrt((*bkgd_out)[j]);
	else var_out[j] = 0.;
    }
    *sigma_out = var_out;

    /*
     *  For HIST data, counts array = weights "uncorrected" for dead time.
     */
     if (!strncmp(instmode,"H",1)) for (j = jmin; j < jmax; j++)
	(*counts_out)[j] = (long) ((*weights_out)[j] / tdead + 0.5);

    /*
     *  Clean up
     */
    free(data);
    free(flux_old);
    free(flux_box);
    free(var_box);
    free(variance);
    free(w);

    cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done processing");
    return status;
}