From d54fe7c1f704a63824c5bfa0ece65245572e9b27 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 4 Mar 2015 21:21:30 -0500 Subject: Initial commit --- src/fuv/cf_extract_spectra.c | 336 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 336 insertions(+) create mode 100644 src/fuv/cf_extract_spectra.c (limited to 'src/fuv/cf_extract_spectra.c') 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 +#include +#include +#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; +} -- cgit