Initial commit

1 files changed, 336 insertions, 0 deletions

diff --git a/src/fuv/cf_extract_spectra.c b/src/fuv/cf_extract_spectra.c
new file mode 100644
index 0000000..6e06119
--- /dev/null
+++ b/src/fuv/cf_extract_spectra.c
@@ -0,0 +1,336 @@
+/*****************************************************************************
+ *              Johns Hopkins University
+ *              Center For Astrophysical Sciences 
+ *              FUSE
+ *****************************************************************************
+ *
+ * Synopsis:    cf_extract_spectra options intermediate_file  
+ *
+ * Description: This module computes a background model, reads the various
+ *		calibration files, and performs either a standard or optimal
+ *		extraction for the target spectrum in each of the LiF and
+ *		SiC channels.
+ *
+ *		Note: some of the routines called by this module modify values
+ *		in the IDF file header.  We want those changes to appear in
+ *		the final, extracted spectral files, but not in the archived
+ *		IDF, so we make a copy in memory of the IDF header and pass it
+ *		to the various routines.  We delete it at the end.
+ *
+ * Arguments:   input_file              Intermediate data file
+ *
+ * Calibration files:                   
+ *
+ * Returns:     0 on successful completion
+ *
+ * Calls:       cf_apply_filters, cf_scale_bkgd, cf_make_wave_array,
+ *		cf_rebin_and_flux_calibrate, cf_rebin_probability_array,
+ *              cf_standard_or_optimal_extraction, cf_optimal_extraction
+ *              cf_write_extracted_spectrum
+ *
+ * History:     02/13/03   peb    1.1   Begin work
+ *		03/02/03   wvd    1.3	Add cf_standard_or_optimal_extraction()
+ *              03/10/03   peb    1.4   Added verbose_level, unistd.h for
+ *                                      getopt portability, and remove debug
+ *                                      print statements
+ *		04/08/03   wvd    1.6	Copy header of IDF into memory
+ *					so that subsequent routines can
+ *					modify it without changing IDF.
+ *					Delete this virtual header when done.
+ *					Read arrays from timeline table and
+ *					pass to cf_apply_filters.
+ *              05/06/03   rdr    1.7   Read from column ERGCM2 rather than
+ *                                      ERGCM2S
+ *              05/07/03   wvd    1.8   Change ergcm2 to ergcm2 throughout
+ *              05/28/03   rdr    1.9   Read pothole mask
+ *              06/09/03   rdr    1.10  Incorporate the tscreen flag
+ *		06/11/03   wvd    1.11  Pass datatype to cf_read_col
+ *		08/25/03   wvd    1.12 	Change coltype from string to int in
+ *					cf_read_col.
+ *		09/29/03   wvd    1.13 	Move standard_or_optimal so that it
+ *					runs just once.  Don't check return
+ *					values from subroutines.
+ *		10/08/03   wvd    1.14 	Change counts_out array to type long.
+ *		10/16/03   wvd    1.15 	If EXPTIME = 0, generate a pair of
+ *					null-valued spectra and exit.
+ *		10/31/03   wvd    1.16 	Change channel to unsigned char.
+ *		12/12/03   wvd    1.18 	Clean up i/o.
+ *		03/16/04   wvd    1.19 	Don't pass wave_out to optimal
+ *					extraction subroutine.
+ *					Change pycent from int to float.
+ *		03/24/04   wvd    1.20 	Eliminate got_no_data(); fold into
+ *					main routine.  If either EXPTIME = 0
+ *					or NEVENTS = 0, write null spectrum.
+ *		04/26/04   wvd    1.21 	Replace cf_rebin_and_flux_calibrate
+ *					background with cf_rebin_background.
+ *					cf_optimal_extraction now returns
+ *					flux_out and sigma_out in units of
+ *					counts; must flux-calibrate each.
+ *		06/10/04   wvd    1.22 	Don't pass time array from timeline
+ *					table to cf_apply_filters.
+ *              06/11/04   peb    1.23  Add the -r option to override the
+ *                                      default rootname.
+ *		01/28/05   wvd    1.24	If EXP_STAT < 0, generate a pair of
+ *                                      null-valued spectra and exit.
+ *		03/09/05   wvd    1.25	Change cf_ttag_bkgd to cf_scale_bkgd,
+ *					pass weights array to it.
+ *		04/12/05   wvd    1.27	Add -n flag, which forces extraction
+ *					of night-only spectrum.  Its argument
+ *					is the name of a night-only BPM file.
+ *		06/15/05   wvd    1.28  BUG FIX: program always read the
+ *					point-source probability arrays from
+ *					the WGTS_CAL file.  Now uses value
+ *					of extended to determine which HDU
+ *					to read.
+ *		11/23/05   wvd    1.29  Add flag to disable optimal 
+ *					extraction from the command line.
+ *		01/27/06   wvd    1.30  Add flag to force use of optimal 
+ *					extraction from the command line.
+ *		05/19/06   wvd    1.31	Don't discard spectra with non-zero
+ *					values of EXP_STAT.
+ *		06/12/06   wvd    1.32	Change cf_if_warning to cf_verbose.
+ *
+ ****************************************************************************/
+
+#include <unistd.h>
+#include <stdlib.h>
+#include <string.h>
+#include "calfuse.h"
+
+static char CF_PRGM_ID[] = "cf_extract_spectra";
+static char CF_VER_NUM[] = "1.32";
+
+int
+main(int argc, char *argv[])
+{
+    char  *bpm_file=NULL, *outrootname = NULL;
+    unsigned char *channel=NULL, *timeflags=NULL, *loc_flags=NULL;
+    unsigned char *time_status=NULL;
+    int   extended, status=0, optc, j, tscreen=1;
+    int   binx, biny, bnx[2], bny[2], bymin[2], aper[2];
+    int   night_only=FALSE, optimal=FALSE;
+    int   force_optimal=FALSE, override_optimal=FALSE;
+    long  nevents=0, ntimes=0, ngood=0, dtime=0, ntime=0, *good_index=NULL;
+    float *weight=NULL, *x=NULL, *y=NULL, *lambda=NULL;
+    float exptime, *bpmask=NULL ;
+    float *bimage[2]={NULL, NULL};
+
+    fitsfile *memp, *header;
+    
+    char opts[]  = "hn:ofr:sv:";
+    char usage[] = 
+	"Usage:\n"
+	"  cf_extract_spectra [-hsof] [-n bpm_filename] [-r rootname] [-v level]"
+		" idf_file\n";
+    char option[] =
+	"Options:\n"
+	"  -h:  this help message\n"
+	"  -n:  extract night-only spectrum using given bad-pixel map\n"
+	"  -o:  disable optimal extraction\n"
+	"  -f:  force optimal extraction\n"
+        "  -r:  override the default rootname\n"
+        "  -s:  do not perform screening on time flags\n" 
+        "  -v:  verbosity level (=1; 0 is silent)\n" ;
+
+    verbose_level = 1;
+
+    while ((optc = getopt(argc, argv, opts)) != -1) {
+	switch(optc) {
+
+	case 'h':
+	    printf("%s\n%s", usage, option);
+	    return 0;
+        case 'n':
+            bpm_file = optarg;
+	    night_only = TRUE;
+            break;
+        case 'o':
+	    override_optimal = TRUE;
+            break;
+        case 'f':
+	    force_optimal = TRUE;
+            break;
+        case 'r':
+            outrootname = optarg;
+            break;
+	case 'v':
+	    verbose_level = atoi(optarg);
+	    break;
+        case 's':
+	  tscreen = 0 ;
+          break;
+	}
+    }
+    cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+    if (argc <= optind) {
+        printf("%s", usage);
+	return -1;
+    }
+    cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");
+
+    /* Open input data file. */
+    FITS_open_file(&header, argv[optind], READONLY, &status);
+    FITS_read_key(header, TFLOAT, "EXPTIME", &exptime, NULL, &status);
+
+    FITS_movabs_hdu(header, 2, NULL, &status);
+    nevents = cf_read_col(header, TFLOAT, "WEIGHT",  (void **) &weight);
+    nevents = cf_read_col(header, TFLOAT, "X",       (void **) &x);
+    nevents = cf_read_col(header, TFLOAT, "Y",       (void **) &y);
+    nevents = cf_read_col(header, TBYTE,  "CHANNEL", (void **) &channel);
+    nevents = cf_read_col(header, TBYTE,  "TIMEFLGS", (void **) &timeflags);
+    nevents = cf_read_col(header, TBYTE,  "LOC_FLGS", (void **) &loc_flags);
+    nevents = cf_read_col(header, TFLOAT, "LAMBDA",  (void **) &lambda);
+
+    FITS_movabs_hdu(header, 4, NULL, &status);
+    ntimes = cf_read_col(header, TBYTE, "STATUS_FLAGS", (void **) &time_status);
+
+   /* 
+    *  Copy header of IDF into memory, close IDF,
+    *  and pass copy of header to subsequent routines.
+    */
+    FITS_movabs_hdu(header, 1, NULL, &status);
+    FITS_create_file(&memp, "mem://", &status);
+    FITS_copy_hdu(header, memp, 0, &status);
+    FITS_close_file(header, &status);
+    header = memp;
+
+    /* If night-only spectrum is requested, modify header accordingly. */
+    if (night_only) {
+	cf_verbose(3, "Night-only spectrum requested from command line.");
+	FITS_read_key(header,  TFLOAT, "EXPNIGHT", &exptime, NULL, &status);
+	FITS_update_key(header, TFLOAT, "EXPTIME",  &exptime, NULL, &status);
+	FITS_update_key(header, TSTRING, "DAYNIGHT", "NIGHT", NULL, &status);
+	FITS_update_key(header, TSTRING, "BPM_CAL", bpm_file, NULL, &status);
+    }
+
+    extended = cf_source_aper(header, aper);
+
+    /* Determine whether to attempt optimal extraction. */
+    if (override_optimal) {
+	cf_verbose (1, "User set -o flag.  No optimal extraction.");
+	optimal=FALSE;
+    }
+    else if (force_optimal) {
+	cf_verbose (1, "User set -f flag.  Forcing use of optimal extraction.");
+	optimal=TRUE;
+    }
+    else cf_standard_or_optimal_extraction(header, &optimal);
+
+    /* Make list of all photons that satisfy requested screenings. */
+    cf_apply_filters(header, tscreen, nevents, timeflags, loc_flags,
+	ntimes, time_status, &dtime, &ntime, &ngood, &good_index);
+
+    /* Generate model backgrounds for LiF and SiC channels. */
+    cf_scale_bkgd(header, nevents, x, y, weight, channel, timeflags,
+	loc_flags, ngood, good_index, dtime, ntime, &binx, &biny,
+	&bnx[0], &bny[0], &bymin[0], &bimage[0],
+	&bnx[1], &bny[1], &bymin[1], &bimage[1]);
+
+    /* Run once for each of LiF and SiC target channels. */
+    for (j=0; j<2; j++) {
+	int   pny, valid_spectrum=TRUE;
+	float pycent, wpc;
+        long *counts_out=NULL, i, nout=0, nphotons=0;
+	float *wave_out=NULL, *bkgd_out=NULL, *flux_out=NULL, *sigma_out=NULL;
+	float *weights_out=NULL, *barray=NULL, *parray=NULL;
+        short *bpix_out=NULL;
+	unsigned char *channel_temp=NULL;
+
+	/* Generate output wavelength array. */
+	cf_make_wave_array(header, aper[j], &nout, &wave_out);
+
+	/* If no data, generate null-valued output arrays. */
+	for (i = 0; i < ngood; i++) {
+	    if (channel[good_index[i]] == aper[j]) {
+		nphotons++;
+		break;
+	    }
+	}
+	if (exptime < 1. || nphotons < 1) {
+            bkgd_out = (float *) cf_calloc(nout, sizeof(float));
+            bpix_out = (short *) cf_calloc(nout, sizeof(short));
+            counts_out = (long *) cf_calloc(nout, sizeof(long));
+            flux_out = (float *) cf_calloc(nout, sizeof(float));
+            sigma_out = (float *) cf_calloc(nout, sizeof(float));
+            weights_out = (float *) cf_calloc(nout, sizeof(float));
+	    valid_spectrum=FALSE;	
+	    cf_verbose(1, "No data for aperture %d.  Generating "
+		"null-valued output spectrum.", aper[j]);
+	}
+	else {			/* Extract target spectrum. */
+
+	/* Bin background model to match output wavelength array. */
+	    cf_rebin_background(header, aper[j], nout, wave_out,
+		binx, biny, bnx[j], bny[j], bimage[j], &barray);
+
+	/* Use QUAL_CAL files to create a bad-pixel mask for this exposure. */
+            cf_make_mask(header, aper[j], nout, wave_out, bny[j], bymin[j],
+		&bpmask);
+
+	/* Bin the weights/probability array to match output wave array. */
+	    cf_rebin_probability_array(header, extended, aper[j], nout,
+		wave_out, &pny, &pycent, &parray);
+
+	/* Use optimal extraction if possible, standard extraction otherwise. */
+	    cf_optimal_extraction(header, optimal, aper[j],
+		weight, y, channel, lambda, ngood, good_index,
+		barray, bpmask, pny, pycent, parray, nout, wave_out, &flux_out,
+		&sigma_out, &counts_out, &weights_out, &bkgd_out, &bpix_out);
+
+	/* Optimal-extraction routine returns flux_out and sigma_out in units
+	   of counts.  Must flux-calibrate both arrays. */
+	    channel_temp =
+		(unsigned char *) cf_malloc(sizeof(unsigned char) * nout);
+	    memset (channel_temp, aper[j], sizeof(unsigned char) * nout);
+	    FITS_read_key(header, TFLOAT, "WPC", &wpc, NULL, &status);
+	    FITS_update_key(header, TSTRING, "FLUX_COR", "PERFORM", NULL,
+		&status) ;
+	    cf_convert_to_ergs(header, nout, flux_out, flux_out, channel_temp,
+		wave_out);
+	    FITS_update_key(header, TSTRING, "FLUX_COR", "PERFORM", NULL,
+		&status) ;
+	    cf_convert_to_ergs(header, nout, sigma_out, sigma_out, channel_temp,
+		wave_out);
+	    free(channel_temp);
+
+	/* Divide flux and errors by EXPTIME * WPC to get units of erg/cm2/s/A. */
+	    for (i = 0; i < nout; i++) {
+		flux_out[i]  /= exptime * wpc;
+		sigma_out[i] /= exptime * wpc;
+	    }
+	}
+
+	/* Write extracted spectrum to output file. */
+	cf_write_extracted_spectrum(header, aper[j], valid_spectrum,
+		nout, wave_out, flux_out, sigma_out, counts_out,
+		weights_out, bkgd_out, bpix_out, outrootname);
+
+	free(counts_out);
+	free(bkgd_out);
+	free(weights_out);
+	free(sigma_out);
+	free(flux_out);
+	free(wave_out);
+        free(bpix_out) ;
+	free(parray);
+	free(barray);
+    }
+
+    cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+
+    FITS_delete_file(header, &status);
+
+    free(bimage[1]);
+    free(bimage[0]);
+    free(lambda);
+    free(loc_flags);
+    free(timeflags);
+    free(channel);
+    free(y);
+    free(x);
+    free(weight);
+    free(time_status);
+
+    cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done processing");
+    return 0;
+}