Initial commit

1 files changed, 148 insertions, 0 deletions

diff --git a/src/analysis/cf_reflux.c b/src/analysis/cf_reflux.c
new file mode 100644
index 0000000..1a2511f
--- /dev/null
+++ b/src/analysis/cf_reflux.c
@@ -0,0 +1,148 @@
+/*****************************************************************************
+ *		Johns Hopkins University
+ *		Center For Astrophysical Sciences
+ *		FUSE
+ *
+ * Synopsis:	cf_reflux input_file output_file aeff_cal1 [aeff_cal2]
+ *
+ * Description: Program applies a new flux calibration to an extracted
+ *		spectral file.  If two effective-area files are given,
+ *		program interpolates between them based on the MDJ
+ *		of the exposure.
+ *		Program writes a history line to output file header.
+ *
+ * Arguments:	input_file	FUSE extracted spectral file (fcal.fit)
+ *		output_file	Copy of input file with modified flux cal.
+ *		aeff_cal1	Effective-area curve
+ *		aeff_cal2	Optional effective-area curve, needed only
+ *				if you wish to interpolate.
+ *
+ * Returns:     none
+ *
+ * History:     06/02/05 1.1    wvd	Based on programs cf_uninterp and
+ *					cf_extract_spectra.
+ *
+ ****************************************************************************/
+
+#include <math.h>
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "calfuse.h"
+
+static char CF_PRGM_ID[] = "cf_reflux";
+static char CF_VER_NUM[] = "1.1";
+
+int main(int argc, char *argv[])
+{
+    char  aper_act[FLEN_VALUE], datestr[FLEN_VALUE];
+    int   status=0, hdutype;
+    int   aperture, timeref;
+    long  i, frow=1, felement=1, nout;
+    float exptime, wpc;
+    float *wave, *flux, *error, *weights, *bkgd;
+    unsigned char *channel;
+    fitsfile *infits, *outfits;
+
+    /* Enter a timestamp into the log. */
+    cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");
+
+    /* Initialize error checking */
+    cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+
+    if (argc < 4) {
+	printf("Usage: cf_reflux input_file output_file aeff_cal1 [aeff_cal2]\n");
+	return 0;
+    }
+
+    /* Copy input spectral file to output.  Close input file. */
+    FITS_open_file(&infits, argv[1], READONLY, &status);
+    FITS_create_file(&outfits, argv[2], &status);
+    FITS_copy_hdu(infits, outfits, 0, &status);
+    FITS_movrel_hdu(infits, 1, &hdutype, &status);
+    FITS_copy_hdu(infits, outfits, 0, &status);
+    FITS_close_file(infits, &status);
+
+    /* Modify output file keywords. */
+    FITS_movabs_hdu(outfits, 1, &hdutype, &status);
+    FITS_update_key(outfits, TSTRING, "FILENAME", argv[2], NULL, &status);
+    FITS_update_key(outfits, TSTRING, "AEFF1CAL", argv[3], NULL, &status);
+    if (argc == 5)
+	FITS_update_key(outfits, TSTRING, "AEFF2CAL", argv[4], NULL, &status);
+    else
+	FITS_update_key(outfits, TSTRING, "AEFF2CAL", argv[3], NULL, &status);
+
+    /* Add HISTORY lines to the output file */
+    fits_get_system_time(datestr, &timeref, &status);
+    strcat(datestr, "   File recalibrated using cf_reflux");
+    FITS_write_history(outfits, datestr, &status);
+
+    /* Through which aperture was this spectrum obtained? */
+    FITS_read_key(outfits, TSTRING, "APER_ACT", aper_act, NULL, &status);
+         if (!strcmp(aper_act, "HIRS_LIF")) aperture = 1;
+    else if (!strcmp(aper_act, "MDRS_LIF")) aperture = 2;
+    else if (!strcmp(aper_act, "LWRS_LIF")) aperture = 3;
+    else if (!strcmp(aper_act, "HIRS_SIC")) aperture = 5;
+    else if (!strcmp(aper_act, "MDRS_SIC")) aperture = 6;
+    else if (!strcmp(aper_act, "LWRS_SIC")) aperture = 7;
+    else {
+	printf("Unable to interpret header keyword APER_ACT = %s\n", aper_act);
+        FITS_close_file(outfits, &status);
+	return 0;
+    }
+
+    /* Read WAVE, FLUX, ERROR, WEIGHTS, and BKGD arrays */
+    FITS_movabs_hdu(outfits, 2, &hdutype, &status);
+    nout = cf_read_col(outfits, TFLOAT, "WAVE", (void **) &wave);
+    nout = cf_read_col(outfits, TFLOAT, "FLUX", (void **) &flux);
+    nout = cf_read_col(outfits, TFLOAT, "ERROR", (void **) &error);
+    nout = cf_read_col(outfits, TFLOAT, "WEIGHTS", (void **) &weights);
+    nout = cf_read_col(outfits, TFLOAT, "BKGD", (void **) &bkgd);
+    channel = (unsigned char *) cf_malloc(sizeof(unsigned char) * nout);
+
+    /* Compute inputs to flux-calibration routine */
+    for (i = 0; i < nout; i++) {
+	channel[i] = (unsigned char) aperture;
+	if (flux[i] != 0.) error[i] /= flux[i];		/* Relative error */
+	flux[i] = weights[i] - bkgd[i];			/* Units are counts */
+    }
+    
+    /* Apply flux calibration */
+    FITS_movabs_hdu(outfits, 1, &hdutype, &status);
+    FITS_update_key(outfits, TSTRING, "FLUX_COR", "PERFORM", NULL, &status) ;
+    cf_convert_to_ergs(outfits, nout, flux, flux, channel, wave);
+
+    /* Convert flux and error arrays to units of erg/cm2/s/A. */
+    FITS_movabs_hdu(outfits, 1, &hdutype, &status);
+    FITS_read_key(outfits, TFLOAT, "EXPTIME", &exptime, NULL, &status);
+    FITS_read_key(outfits, TFLOAT, "WPC", &wpc, NULL, &status);
+    if (exptime < 1) {
+	printf ("EXPTIME = %g.  Exiting.\n", exptime);
+        FITS_close_file(outfits, &status);
+	return 0;
+    }
+    for (i = 0; i < nout; i++) {
+        flux[i]  /= exptime * wpc;
+        error[i] *= flux[i];
+    }
+ 
+    /* Write new flux and error arrays to the output file. */
+    FITS_movabs_hdu(outfits, 2, &hdutype, &status);
+    FITS_write_col(outfits, TFLOAT, 2, frow, felement, nout, flux, &status);
+    FITS_write_col(outfits, TFLOAT, 3, frow, felement, nout, error, &status);
+ 
+    /* Close output file. */
+    FITS_close_file(outfits, &status);
+
+    /* Free memory. */
+    free(wave);
+    free(flux);
+    free(error);
+    free(weights);
+    free(bkgd);
+    free(channel);
+
+    /* Enter a timestamp into the log. */
+    cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done processing");
+    return 0;
+}