aboutsummaryrefslogtreecommitdiff
path: root/src/fuv
diff options
context:
space:
mode:
authorJoseph Hunkeler <jhunkeler@gmail.com>2015-03-04 21:21:30 -0500
committerJoseph Hunkeler <jhunkeler@gmail.com>2015-03-04 21:21:30 -0500
commitd54fe7c1f704a63824c5bfa0ece65245572e9b27 (patch)
treeafc52015ffc2c74e0266653eecef1c8ef8ba5d91 /src/fuv
downloadcalfuse-d54fe7c1f704a63824c5bfa0ece65245572e9b27.tar.gz
Initial commit
Diffstat (limited to 'src/fuv')
-rw-r--r--src/fuv/Makefile.Linux.orig75
-rw-r--r--src/fuv/Makefile.Linux64.orig74
-rw-r--r--src/fuv/Makefile.MacOSX.orig76
-rw-r--r--src/fuv/Makefile.Solaris.orig75
-rw-r--r--src/fuv/Makefile.orig.orig75
-rw-r--r--src/fuv/cf_assign_wavelength.c178
-rw-r--r--src/fuv/cf_bad_pixels.c1001
-rw-r--r--src/fuv/cf_convert_to_farf.c268
-rw-r--r--src/fuv/cf_countmap.c541
-rw-r--r--src/fuv/cf_extract_spectra.c336
-rw-r--r--src/fuv/cf_flux_calibrate.c173
-rw-r--r--src/fuv/cf_gainmap.c610
-rw-r--r--src/fuv/cf_hist_init.c653
-rw-r--r--src/fuv/cf_remove_motions.c316
-rw-r--r--src/fuv/cf_screen_photons.c298
-rw-r--r--src/fuv/cf_ttag_init.c594
16 files changed, 5343 insertions, 0 deletions
diff --git a/src/fuv/Makefile.Linux.orig b/src/fuv/Makefile.Linux.orig
new file mode 100644
index 0000000..da60b24
--- /dev/null
+++ b/src/fuv/Makefile.Linux.orig
@@ -0,0 +1,75 @@
+
+CALFUSEDIR= ${PWD}/../..
+SHARED= -shared
+FITSVER= 2.470
+
+# Symbols for include directories
+FUSEINCLDIR= -I${CALFUSEDIR}/include
+
+# Symbols used for compiling
+CC= cc
+OPT= -g -Wall -DCFORTRAN -Dg77Fortran -Df2cFortran
+CFLAGS= ${OPT} ${FUSEINCLDIR}
+
+FUSEBINDIR= ${CALFUSEDIR}/bin
+FUSELIBDIR= -L${CALFUSEDIR}/lib
+FUSELIBS= -lcfitsio-${FITSVER} -lsla -lcf
+LIBS= -lc -lm -lnsl -ldl -lgfortran
+LDFLAGS= -Wl,-R${CALFUSEDIR}/lib,-R${CALFUSEDIR}/src/libcf
+
+# Symbols used for creating shared libraries
+
+BINS= cf_hist_init cf_ttag_init cf_convert_to_farf cf_remove_motions \
+ cf_assign_wavelength cf_screen_photons cf_flux_calibrate \
+ cf_extract_spectra cf_countmap cf_gainmap cf_bad_pixels
+
+all: ${BINS}
+ chmod g+w ${BINS}
+
+install: all
+ /bin/cp -p ${BINS} ${FUSEBINDIR}
+
+clean:
+ - /bin/rm -f ${BINS}
+
+distclean:
+ - /bin/rm -f ${BINS}
+ cd ../../bin; /bin/rm -f ${BINS} calfuse
+
+cf_hist_init:
+ ${CC} ${CFLAGS} -o cf_hist_init cf_hist_init.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_ttag_init:
+ ${CC} ${CFLAGS} -o cf_ttag_init cf_ttag_init.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_convert_to_farf:
+ ${CC} ${CFLAGS} -o cf_convert_to_farf cf_convert_to_farf.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_remove_motions:
+ ${CC} ${CFLAGS} -o cf_remove_motions cf_remove_motions.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_assign_wavelength:
+ ${CC} ${CFLAGS} -o cf_assign_wavelength \
+ cf_assign_wavelength.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_screen_photons:
+ ${CC} ${CFLAGS} -o cf_screen_photons cf_screen_photons.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_flux_calibrate:
+ ${CC} ${CFLAGS} -o cf_flux_calibrate cf_flux_calibrate.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_extract_spectra:
+ ${CC} ${CFLAGS} -o cf_extract_spectra cf_extract_spectra.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_countmap:
+ ${CC} ${CFLAGS} -o cf_countmap cf_countmap.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+
+cf_gainmap:
+ ${CC} ${CFLAGS} -o cf_gainmap cf_gainmap.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+
+cf_bad_pixels:
+ ${CC} ${CFLAGS} -o cf_bad_pixels cf_bad_pixels.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+
diff --git a/src/fuv/Makefile.Linux64.orig b/src/fuv/Makefile.Linux64.orig
new file mode 100644
index 0000000..6fbead0
--- /dev/null
+++ b/src/fuv/Makefile.Linux64.orig
@@ -0,0 +1,74 @@
+
+CALFUSEDIR= ${PWD}/../..
+SHARED= -shared
+
+# Symbols for include directories
+FUSEINCLDIR= -I${CALFUSEDIR}/include
+
+# Symbols used for compiling
+CC= cc
+OPT= -g -Wall -DCFORTRAN -Dg77Fortran -Df2cFortran
+CFLAGS= ${OPT} ${FUSEINCLDIR}
+
+FUSEBINDIR= ${CALFUSEDIR}/bin
+FUSELIBDIR= -L${CALFUSEDIR}/lib
+FUSELIBS= -lsla -lcf
+LIBS= -lc -lm -lnsl -ldl -lgfortran -lcfitsio
+LDFLAGS= -Wl,-R${CALFUSEDIR}/lib,-R${CALFUSEDIR}/src/libcf
+
+# Symbols used for creating shared libraries
+
+BINS= cf_hist_init cf_ttag_init cf_convert_to_farf cf_remove_motions \
+ cf_assign_wavelength cf_screen_photons cf_flux_calibrate \
+ cf_extract_spectra cf_countmap cf_gainmap cf_bad_pixels
+
+all: ${BINS}
+ chmod g+w ${BINS}
+
+install: all
+ /bin/cp -p ${BINS} ${FUSEBINDIR}
+
+clean:
+ - /bin/rm -f ${BINS}
+
+distclean:
+ - /bin/rm -f ${BINS}
+ cd ../../bin; /bin/rm -f ${BINS} calfuse
+
+cf_hist_init:
+ ${CC} ${CFLAGS} -o cf_hist_init cf_hist_init.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_ttag_init:
+ ${CC} ${CFLAGS} -o cf_ttag_init cf_ttag_init.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_convert_to_farf:
+ ${CC} ${CFLAGS} -o cf_convert_to_farf cf_convert_to_farf.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_remove_motions:
+ ${CC} ${CFLAGS} -o cf_remove_motions cf_remove_motions.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_assign_wavelength:
+ ${CC} ${CFLAGS} -o cf_assign_wavelength \
+ cf_assign_wavelength.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_screen_photons:
+ ${CC} ${CFLAGS} -o cf_screen_photons cf_screen_photons.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_flux_calibrate:
+ ${CC} ${CFLAGS} -o cf_flux_calibrate cf_flux_calibrate.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_extract_spectra:
+ ${CC} ${CFLAGS} -o cf_extract_spectra cf_extract_spectra.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_countmap:
+ ${CC} ${CFLAGS} -o cf_countmap cf_countmap.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+
+cf_gainmap:
+ ${CC} ${CFLAGS} -o cf_gainmap cf_gainmap.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+
+cf_bad_pixels:
+ ${CC} ${CFLAGS} -o cf_bad_pixels cf_bad_pixels.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+
diff --git a/src/fuv/Makefile.MacOSX.orig b/src/fuv/Makefile.MacOSX.orig
new file mode 100644
index 0000000..acbfaf7
--- /dev/null
+++ b/src/fuv/Makefile.MacOSX.orig
@@ -0,0 +1,76 @@
+
+CALFUSEDIR= ${PWD}/../..
+FITSVER= 2.470
+MACOSX_DEPLOYMENT_TARGET= 10.2
+
+# Symbols for include directories
+FUSEINCLDIR= -I${CALFUSEDIR}/include
+
+# Symbols used for compiling
+CC= cc
+OPT= -O3 -Wall -DCFORTRAN -Dg77Fortran -Df2cFortran
+CFLAGS= ${OPT} ${FUSEINCLDIR}
+
+FUSEBINDIR= ${CALFUSEDIR}/bin
+FUSELIBDIR= -L${CALFUSEDIR}/lib
+# FUSELIBS= -lcfitsio-${FITSVER} -lsla -lcf
+FUSELIBS= -lcfitsio -lsla -lcf
+LIBS= -lc -lm -ldl -L/sw/lib/ -lgfortran
+
+BINS= cf_hist_init cf_ttag_init cf_convert_to_farf cf_remove_motions \
+ cf_assign_wavelength cf_screen_photons cf_flux_calibrate \
+ cf_extract_spectra cf_countmap cf_gainmap cf_bad_pixels
+
+all: ${BINS}
+ chmod g+w ${BINS}
+
+install:
+ MACOSX_DEPLOYMENT_TARGET="${MACOSX_DEPLOYMENT_TARGET}"; \
+ export MACOSX_DEPLOYMENT_TARGET; \
+ make all
+ /bin/cp -p ${BINS} ${FUSEBINDIR}
+
+clean:
+ - /bin/rm -f ${BINS}
+
+distclean:
+ - /bin/rm -f ${BINS}
+ cd ../../bin; /bin/rm -f ${BINS} calfuse
+
+cf_hist_init:
+ ${CC} ${CFLAGS} -o cf_hist_init cf_hist_init.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS}
+cf_ttag_init:
+ ${CC} ${CFLAGS} -o cf_ttag_init cf_ttag_init.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS}
+cf_convert_to_farf:
+ ${CC} ${CFLAGS} -o cf_convert_to_farf cf_convert_to_farf.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS}
+cf_remove_motions:
+ ${CC} ${CFLAGS} -o cf_remove_motions cf_remove_motions.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS}
+cf_assign_wavelength:
+ ${CC} ${CFLAGS} -o cf_assign_wavelength \
+ cf_assign_wavelength.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS}
+cf_screen_photons:
+ ${CC} ${CFLAGS} -o cf_screen_photons cf_screen_photons.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS}
+cf_flux_calibrate:
+ ${CC} ${CFLAGS} -o cf_flux_calibrate cf_flux_calibrate.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS}
+cf_extract_spectra:
+ ${CC} ${CFLAGS} -o cf_extract_spectra cf_extract_spectra.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS}
+cf_countmap:
+ ${CC} ${CFLAGS} -o cf_countmap cf_countmap.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS}
+
+cf_gainmap:
+ ${CC} ${CFLAGS} -o cf_gainmap cf_gainmap.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS}
+
+cf_bad_pixels:
+ ${CC} ${CFLAGS} -o cf_bad_pixels cf_bad_pixels.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS}
+
diff --git a/src/fuv/Makefile.Solaris.orig b/src/fuv/Makefile.Solaris.orig
new file mode 100644
index 0000000..528e6bc
--- /dev/null
+++ b/src/fuv/Makefile.Solaris.orig
@@ -0,0 +1,75 @@
+
+CALFUSEDIR= ${PWD}/../..
+SHARED= -G
+FITSVER= 2.470
+
+# Symbols for include directories
+FUSEINCLDIR= -I${CALFUSEDIR}/include
+
+# Symbols used for compiling
+CC= cc
+OPT= -O -DCFORTRAN -KPIC
+CFLAGS= ${OPT} ${FUSEINCLDIR}
+
+FUSEBINDIR= ${CALFUSEDIR}/bin
+FUSELIBDIR= -L${CALFUSEDIR}/lib
+FUSELIBS= -lcfitsio-${FITSVER} -lsla -lcf
+LIBS= -lc -lm -lnsl -ldl -lsocket -lsunmath
+LDFLAGS= -Wl,-R${CALFUSEDIR}/lib,-R${CALFUSEDIR}/src/libcf
+
+# Symbols used for creating shared libraries
+
+BINS= cf_hist_init cf_ttag_init cf_convert_to_farf cf_remove_motions \
+ cf_assign_wavelength cf_screen_photons cf_flux_calibrate \
+ cf_extract_spectra cf_countmap cf_gainmap cf_bad_pixels
+
+all: ${BINS}
+ chmod g+w ${BINS}
+
+install: all
+ /bin/cp -p ${BINS} ${FUSEBINDIR}
+
+clean:
+ - /bin/rm -f ${BINS}
+
+distclean:
+ - /bin/rm -f ${BINS}
+ cd ../../bin; /bin/rm -f ${BINS} calfuse
+
+cf_hist_init:
+ ${CC} ${CFLAGS} -o cf_hist_init cf_hist_init.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_ttag_init:
+ ${CC} ${CFLAGS} -o cf_ttag_init cf_ttag_init.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_convert_to_farf:
+ ${CC} ${CFLAGS} -o cf_convert_to_farf cf_convert_to_farf.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_remove_motions:
+ ${CC} ${CFLAGS} -o cf_remove_motions cf_remove_motions.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_assign_wavelength:
+ ${CC} ${CFLAGS} -o cf_assign_wavelength \
+ cf_assign_wavelength.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_screen_photons:
+ ${CC} ${CFLAGS} -o cf_screen_photons cf_screen_photons.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_flux_calibrate:
+ ${CC} ${CFLAGS} -o cf_flux_calibrate cf_flux_calibrate.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_extract_spectra:
+ ${CC} ${CFLAGS} -o cf_extract_spectra cf_extract_spectra.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_countmap:
+ ${CC} ${CFLAGS} -o cf_countmap cf_countmap.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+
+cf_gainmap:
+ ${CC} ${CFLAGS} -o cf_gainmap cf_gainmap.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+
+cf_bad_pixels:
+ ${CC} ${CFLAGS} -o cf_bad_pixels cf_bad_pixels.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+
diff --git a/src/fuv/Makefile.orig.orig b/src/fuv/Makefile.orig.orig
new file mode 100644
index 0000000..528e6bc
--- /dev/null
+++ b/src/fuv/Makefile.orig.orig
@@ -0,0 +1,75 @@
+
+CALFUSEDIR= ${PWD}/../..
+SHARED= -G
+FITSVER= 2.470
+
+# Symbols for include directories
+FUSEINCLDIR= -I${CALFUSEDIR}/include
+
+# Symbols used for compiling
+CC= cc
+OPT= -O -DCFORTRAN -KPIC
+CFLAGS= ${OPT} ${FUSEINCLDIR}
+
+FUSEBINDIR= ${CALFUSEDIR}/bin
+FUSELIBDIR= -L${CALFUSEDIR}/lib
+FUSELIBS= -lcfitsio-${FITSVER} -lsla -lcf
+LIBS= -lc -lm -lnsl -ldl -lsocket -lsunmath
+LDFLAGS= -Wl,-R${CALFUSEDIR}/lib,-R${CALFUSEDIR}/src/libcf
+
+# Symbols used for creating shared libraries
+
+BINS= cf_hist_init cf_ttag_init cf_convert_to_farf cf_remove_motions \
+ cf_assign_wavelength cf_screen_photons cf_flux_calibrate \
+ cf_extract_spectra cf_countmap cf_gainmap cf_bad_pixels
+
+all: ${BINS}
+ chmod g+w ${BINS}
+
+install: all
+ /bin/cp -p ${BINS} ${FUSEBINDIR}
+
+clean:
+ - /bin/rm -f ${BINS}
+
+distclean:
+ - /bin/rm -f ${BINS}
+ cd ../../bin; /bin/rm -f ${BINS} calfuse
+
+cf_hist_init:
+ ${CC} ${CFLAGS} -o cf_hist_init cf_hist_init.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_ttag_init:
+ ${CC} ${CFLAGS} -o cf_ttag_init cf_ttag_init.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_convert_to_farf:
+ ${CC} ${CFLAGS} -o cf_convert_to_farf cf_convert_to_farf.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_remove_motions:
+ ${CC} ${CFLAGS} -o cf_remove_motions cf_remove_motions.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_assign_wavelength:
+ ${CC} ${CFLAGS} -o cf_assign_wavelength \
+ cf_assign_wavelength.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_screen_photons:
+ ${CC} ${CFLAGS} -o cf_screen_photons cf_screen_photons.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_flux_calibrate:
+ ${CC} ${CFLAGS} -o cf_flux_calibrate cf_flux_calibrate.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_extract_spectra:
+ ${CC} ${CFLAGS} -o cf_extract_spectra cf_extract_spectra.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+cf_countmap:
+ ${CC} ${CFLAGS} -o cf_countmap cf_countmap.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+
+cf_gainmap:
+ ${CC} ${CFLAGS} -o cf_gainmap cf_gainmap.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+
+cf_bad_pixels:
+ ${CC} ${CFLAGS} -o cf_bad_pixels cf_bad_pixels.c \
+ ${FUSELIBDIR} ${FUSELIBS} ${LIBS} ${LDFLAGS}
+
diff --git a/src/fuv/cf_assign_wavelength.c b/src/fuv/cf_assign_wavelength.c
new file mode 100644
index 0000000..d84c7a0
--- /dev/null
+++ b/src/fuv/cf_assign_wavelength.c
@@ -0,0 +1,178 @@
+/*****************************************************************************
+ * Johns Hopkins University
+ * Center For Astrophysical Sciences
+ * FUSE
+ *****************************************************************************
+ *
+ * Synopsis: cf_assign_wavelength options intermediate_data_file
+ *
+ * Description: Assign a wavelength to each photon, correcting for
+ * spacecraft motion and instrumental astigmatism.
+ * The final wavelength scale is heliocentric.
+ *
+ * By default all corrections are performed. Command-line
+ * options allow one or more corrections to be omitted.
+ *
+ * Arguments: input_file FARF-corrected intermediate data file
+ *
+ * Calibration files: ASTG_CAL, WAVE_CAL
+ *
+ * Returns: 0 on successful completion
+ *
+ * Calls: cf_astigmatism, cf_dispersion, cf_doppler_and_heliocentric
+ *
+ * History: 11/05/02 1.1 peb Begin work
+ * 12/16/02 1.2 wvd Install subroutines.
+ * 02/12/03 1.3 wvd Change ORBITAL_VELOCITY to ORBITAL_VEL
+ * 03/10/03 1.4 peb Added verbose_level, changed channel to
+ * char *, and added unistd.h for getopt
+ * portability
+ * 03/12/03 1.5 wvd If lambda[k] is undefined,
+ * set channel[k] = 0 in IDF
+ * 03/23/03 1.6 wvd If -w option is set, write astigmatism-
+ * corrected X values to IDF.
+ * 06/11/03 1.7 wvd Pass datatype to cf_read_col and
+ * cf_write_col.
+ * 08/25/03 1.8 wvd Change coltype from string to int in
+ * cf_read_col and cf_write_col.
+ * 09/16/03 1.9 wvd Write comment lines to IDF header if
+ * -w option is requested.
+ * 10/31/03 1.10 wvd Change channel to unsigned char.
+ * 12/03/03 1.11 wvd Initialize astig_return to 1.
+ * 04/06/04 1.12 bjg Function returns EXIT_SUCCESS
+ * Fix option[]
+ * 06/21/04 1.13 wvd If -w flag is set, run the program
+ * even if it's been run before.
+ * 05/20/05 1.14 wvd Clean up i/o.
+ * 05/15/06 1.15 wvd Divide cf_astigmatism_and_dispersion
+ * into two routines. Note that wave-
+ * length calibration is performed even
+ * if astigmatism correction is not.
+ *
+ ****************************************************************************/
+
+#include <unistd.h>
+#include <stdlib.h>
+#include "calfuse.h"
+
+int
+main(int argc, char *argv[])
+{
+ char CF_PRGM_ID[] = "cf_assign_wavelength";
+ char CF_VER_NUM[] = "1.15";
+
+ char comment[FLEN_CARD], datestr[FLEN_CARD];
+ unsigned char *channel=NULL;
+ int astig_corr=1, doppler_corr=1, status=0, optc, timeref;
+ int astig_return=1, wavecal = FALSE;
+ long nevents, nseconds;
+ float *time=NULL, *x=NULL, *y=NULL, *lambda=NULL, *timeline=NULL;
+ float *velocity=NULL;
+ fitsfile *header;
+
+ char opts[] = "hadwv:";
+ char usage[] =
+ "Usage:\n"
+ " cf_assign_wavelength [-hadw] [-v level] idf_file\n";
+ char option[] =
+ "Options:\n"
+ " -h: this help message\n"
+ " -v: verbose mode (default is 1; 0 is silent)\n"
+ " -a: no astigmatism correction\n"
+ " -d: no doppler correction\n"
+ " -w: write astigmatism-corrected X values to IDF\n";
+
+ verbose_level = 1;
+
+ while ((optc = getopt(argc, argv, opts)) != -1) {
+ switch(optc) {
+
+ case 'h':
+ printf("%s\n%s", usage, option);
+ return 0;
+ case 'v':
+ verbose_level = atoi(optarg);
+ break;
+ case 'a':
+ astig_corr = 0;
+ break;
+ case 'd':
+ doppler_corr = 0;
+ break;
+ case 'w':
+ wavecal = TRUE;
+ break;
+ }
+ }
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+ if (argc <= optind) {
+ printf("%s", usage);
+ cf_if_error("Incorrect number of program arguments");
+ }
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");
+
+ FITS_open_file(&header, argv[optind], READWRITE, &status);
+
+ FITS_movabs_hdu(header, 2, NULL, &status);
+ nevents = cf_read_col(header, TFLOAT, "TIME", (void **) &time);
+ 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, TFLOAT, "LAMBDA", (void **) &lambda);
+
+ FITS_movabs_hdu(header, 4, NULL, &status);
+ nseconds = cf_read_col(header, TFLOAT, "TIME", (void **) &timeline);
+ nseconds = cf_read_col(header, TFLOAT, "ORBITAL_VEL", (void **) &velocity);
+ FITS_movabs_hdu(header, 1, NULL, &status);
+
+ if (wavecal) {
+ FITS_update_key(header, TSTRING, "WAVE_COR", "PERFORM", NULL, &status);
+ FITS_update_key(header, TSTRING, "DOPP_COR", "PERFORM", NULL, &status);
+ }
+
+ /* Correct X coordinate for 2-D astigmatism. */
+ if (astig_corr)
+ astig_return = cf_astigmatism(header, nevents, x, y, channel);
+
+ /* Assign wavelength to each photon event according to channel and X */
+ cf_dispersion(header, nevents, x, channel, lambda);
+
+ /* Correct for spacecraft motion and shift to heliocentric wavelength scale. */
+ if (doppler_corr)
+ cf_doppler_and_heliocentric(header, nevents, time, channel,
+ lambda, nseconds, timeline, velocity);
+
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+
+ FITS_movabs_hdu(header, 2, NULL, &status);
+ cf_write_col(header, TBYTE, "CHANNEL", (void *) channel, nevents);
+ cf_write_col(header, TFLOAT, "LAMBDA", (void *) lambda, nevents);
+
+ /* If -w option is requested and astigmatism correction was applied,
+ write astigmatism-corrected X array to IDF. */
+ if (wavecal && !astig_return) {
+ cf_write_col(header, TFLOAT, "X", (void *) x, nevents);
+ FITS_movabs_hdu(header, 1, NULL, &status);
+ sprintf(comment, "Writing astigmatism-corrected X array to IDF.");
+ cf_if_warning(comment);
+ FITS_write_comment(header, " ", &status);
+ FITS_write_comment(header, comment, &status);
+ fits_get_system_time(datestr, &timeref, &status);
+ sprintf(comment, "CalFUSE v%s %.10s", CALFUSE_VERSION, datestr);
+ FITS_write_comment(header, comment, &status);
+ FITS_write_comment(header, " ", &status);
+ }
+
+ FITS_close_file(header, &status);
+
+ free(velocity);
+ free(timeline);
+ free(lambda);
+ free(channel);
+ free(y);
+ free(x);
+ free(time);
+
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done processing");
+ return EXIT_SUCCESS;
+}
diff --git a/src/fuv/cf_bad_pixels.c b/src/fuv/cf_bad_pixels.c
new file mode 100644
index 0000000..a10ed20
--- /dev/null
+++ b/src/fuv/cf_bad_pixels.c
@@ -0,0 +1,1001 @@
+/*****************************************************************************
+ * Johns Hopkins University
+ * Center For Astrophysical Sciences
+ * FUSE
+ *****************************************************************************
+ *
+ * Synopsis: cf_bad_pixels options input_file
+ *
+ * Description: Determines a map of the bad pixels on the detector, after
+ * correction for image motions during the observation. The
+ * procedure works by generating a series of pseudo photons
+ * corresponding to bad pixels within the area of the detector
+ * seen by the target aperture. This list of pseudo photons is
+ * repeated for each second of the observation where the data
+ * is expected to be acceptable (i.e., it excludes times of
+ * bursts, SAA passages, etc). The list is then run through
+ * various routines in the pipeline which correct for photon
+ * motions. A map is then generated from the final photon
+ * list.
+ *
+ * Arguments: input_file FARF-corrected intermediate data
+ * file
+ * Calibration files: QUAL_CAL
+ *
+ * Returns: file containing the bad pixel map
+ *
+ *
+ * History: 04/08/03 1.0 rdr Begin work
+ * 04/30/03 1.1 wvd Install
+ * 05/22/03 1.2 wvd Modify call to cf_set_photon_flags
+ * 05/28/03 1.3 rdr Made adjustments for HIST data and
+ * correct some bugs
+ * 05/04/03 1.4 wvd Reset JITR_COR and FLAG_COR to PERFORM.
+ * Pass weights to cf_set_photon_flags.
+ * 09/06/03 1.5 rdr Incorporate tscreen flag.
+ * 06/11/03 1.6 wvd Pass datatype to cf_read_col and
+ * cf_write_col.
+ * 08/22/03 1.7 wvd Change channel array to unsigned char.
+ * Delete GTI from call to
+ * cf_satellite_jitter.
+ * 08/25/03 1.8 wvd Change coltype from string to int in
+ * cf_read_col.
+ * 10/17/03 1.9 wvd If EXPTIME = 0, exit without
+ * generating a bad-pixel file.
+ * 11/05/03 1.10 wvd Change chan_pix to unsigned char
+ * 11/05/03 1.11 wvd Use fits_create_tbl alone to make
+ * binary table extension.
+ * 12/01/03 1.12 bjg Update FILENAME, FILETYPE and IDF_FILE
+ * keywords.
+ * 12/01/03 1.13 bjg Minor changes.
+ * 12/21/03 1.14 wvd Remove underscore from idf and bpm
+ * filenames.
+ * 02/24/04 1.15 rdr change maximum size of the output
+ * arrays in cf_combine_pothole_data
+ * 03/03/04 1.16 rdr make sure that lif_cnt and sic_cnt arrays
+ * are non-zero
+ * 04/06/04 1.19 bjg Include ctype.h and strings.h
+ * Change formats to match arg types in
+ * printf
+ * Change ( a<b & a>c ) to ((a<b)&&(a>c))
+ * in cf_combine_pothole_data
+ * Create potholes from HIST FILL_DATA
+ * 06/07/04 1.20 wvd Add exp_jitr to cf_satellite_jitter(),
+ * delete timeline from cf_apply_filters.
+ * 07/21/04 1.21 wvd Delete unused variable nmax in
+ * cf_generate_pseudo_photons
+ * 10/29/04 1.22 bjg Added check for selected potholes that
+ * fall out the extraction window
+ * 12/02/04 1.23 wvd If a processing step was skipped for
+ * the data file, skip it also for the
+ * bpm file.
+ * 01/27/05 1.24 wvd Ignore LIMB, SAA, and BRST when
+ * determining good times for HIST data.
+ * 02/01/05 1.25 wvd In cf_generate_pseudo_photons, pad
+ * extraction windows by ten pixels
+ * when making pothole list.
+ * In cf_combine_pothole_data, toss any
+ * potholes that drift out of extraction
+ * window. Call cf_get_extraction_limits
+ * only when the aperture changes.
+ * 02/17/05 1.26 wvd Increase estimated size of output
+ * pixel arrays (nmax) by 20%.
+ * Add additional diagnostic output.
+ * 03/15/05 1.27 wvd Replace cf_get_extraction_limits with
+ * cf_extraction_limits
+ * 03/16/05 1.28 wvd Pass srctype to cf_extraction_limits.
+ * 03/22/05 1.29 wvd Change TIME_SUNRISE and TIME_SUNSET
+ * from floats to shorts.
+ * 03/22/05 1.30 wvd In cf_combine_pothole_data, use fabs()
+ * when looking for changes in the time
+ * array.
+ * 04/12/05 1.31 wvd Add -n flag, which forces creation
+ * of night-only BPM file. Its argument
+ * is the name of the output BPM file.
+ * This file name is NOT written to the
+ * IDF file header.
+ * 06/03/05 1.32 wvd Include math.h
+ * 11/25/05 1.33 wvd Use cf_x2lambda and cf_get_potholes
+ * rather than writing new routines here.
+ * Don't add random numbers to output
+ * pixel coordinates.
+ * 05/15/06 1.34 wvd Divide cf_astigmatism_and_dispersion
+ * into two routines. Note that wave-
+ * length calibration is performed even
+ * if astigmatism correction is not.
+ * Add ASTG_COR to list of updated
+ * keywords.
+ * 06/13/06 1.35 wvd In cf_combine_pothole_data,
+ * initialize nfill to 0 before each
+ * iteration and ignore bad pixels that
+ * lie outside of the extraction window.
+ * 03/20/07 1.36 wvd In cf_combine_pothole_data, if the
+ * entire dead spot falls outside of the
+ * aperture, stop and move on to the
+ * next one.
+ * 03/21/07 1.37 wvd Replace "break" with "continue" in
+ * cf_combine_pothole_data.
+ * 04/03/07 1.38 wvd Scale nmax by 1.5. Call cf_error_init
+ * after each call to cf_extraction_limits.
+ * Make CF_PRGM_ID and CF_VER_NUM static
+ * variables.
+ * 04/03/07 1.39 bot In cf_generate_pseduo_photons, l.198
+ * and l.219 added a test on LOCATION_SHLD
+ * for FILL DATA.
+ * 07/18/08 1.40 wvd If EXP_STAT = 2, target is bright
+ * earth/airglow. Mask out limb-angle flag.
+ *
+ ****************************************************************************/
+
+#include <strings.h>
+#include <string.h>
+#include <stdio.h>
+#include <unistd.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <math.h>
+#include "calfuse.h"
+
+static char CF_PRGM_ID[] = "cf_bad_pixels";
+static char CF_VER_NUM[] = "1.40";
+
+
+/****************************************************************************
+ *
+ * CF_GENERATE_PSEUDO_PHOTONS
+ *
+ * Procedure to generate an array of pseudo photons to represent the
+ * locations of the detector potholes
+ *
+ ****************************************************************************/
+
+int cf_generate_pseudo_photons(fitsfile *header, long ngood, long *good_index,
+ float *timeline, unsigned char *statflag, float *lif_cnt,
+ float *sic_cnt, long *nevents, float **time, float **weight,
+ float **x, float **y, unsigned char **channel, float **rx,
+ float **ry, unsigned char **timeflag, unsigned char **loc_flag) {
+
+ char instmode[FLEN_CARD] ;
+ long nevts, npot, nfill=0, npot_sel;
+ int *bpndx, npot2, status=0, hdutype ;
+ int srctype, ap[2], bymax, bymin, pad=10 ;
+ short *ylow, *yhigh;
+ long npts, i, j, k, xndx, num, num_tot, ndx, tndx ;
+ float *xpot, *ypot, *rxpot, *rypot ;
+ float *xpot2, *ypot2, *rxpot2, *rypot2 ;
+ float *xval, *yval, *rxval, *ryval, *xout, *yout, *rxout, *ryout ;
+ float *time_out, *weight_out, tval, ctot_lif, ctot_sic ;
+ float wval_lif, wval_sic ;
+ short binx, biny, xmin, xmax;
+ unsigned char *chan, *chan_out;
+ unsigned char *tflag_out, tflag_val ;
+ char * hdu2_loc_flgs;
+ float * hdu2_xfarf, * hdu2_yfarf;
+
+
+ cf_verbose(3, "Generating pseudo-photons") ;
+
+ /* read the instrument mode for the data */
+ FITS_movabs_hdu(header, 1, &hdutype, &status);
+ FITS_read_key(header, TSTRING, "INSTMODE", instmode, NULL, &status) ;
+ FITS_read_key(header, TSHORT, "SPECBINX", &binx, NULL, &status) ;
+ FITS_read_key(header, TSHORT, "SPECBINY", &biny, NULL, &status) ;
+
+ /* read the header and determine the target apertures for the exposure */
+ srctype = cf_source_aper(header, ap) ;
+ cf_verbose(3,"Target apertures = %d and %d ",ap[0],ap[1]) ;
+
+ /* read in the pothole positions */
+ npot2 = cf_get_potholes(header, &xpot2, &ypot2, &rxpot2, &rypot2) ;
+ cf_verbose(3, "number of potholes found = %d ",npot2) ;
+
+ if (!(strncasecmp(instmode,"HIST",4))) {
+ FITS_movabs_hdu(header, 2, &hdutype, &status);
+ nevts=cf_read_col(header,TBYTE,"LOC_FLGS",(void **) &hdu2_loc_flgs);
+ nevts=cf_read_col(header,TFLOAT,"XFARF",(void **) &hdu2_xfarf);
+ nevts=cf_read_col(header,TFLOAT,"YFARF",(void **) &hdu2_yfarf);
+ for (j=0;j<nevts;++j) {
+ if (!((hdu2_loc_flgs[j] & LOCATION_FILL) == 0) &&
+ ((hdu2_loc_flgs[j] & LOCATION_SHLD) == 0)) ++nfill;
+ }
+ cf_verbose(3, "number of FILL DATA potholes found = %d ",nfill) ;
+ npot=npot2+nfill;
+
+ xpot=(float *) cf_calloc(npot, sizeof(float)) ;
+ ypot=(float *) cf_calloc(npot, sizeof(float)) ;
+ rxpot=(float *) cf_calloc(npot, sizeof(float)) ;
+ rypot=(float *) cf_calloc(npot, sizeof(float)) ;
+
+ for (j=0;j<npot2;++j){
+ xpot[j]=xpot2[j];
+ ypot[j]=ypot2[j];
+ rxpot[j]=rxpot2[j];
+ rypot[j]=rypot2[j];
+ }
+
+ k=npot2;
+
+ for (j=0;j<nevts;++j){
+ if (!((hdu2_loc_flgs[j] & LOCATION_FILL) == 0) &&
+ ((hdu2_loc_flgs[j] & LOCATION_SHLD) == 0)) {
+
+ xpot[k]=hdu2_xfarf[j];
+ ypot[k]=hdu2_yfarf[j];
+ rxpot[k]=binx;
+ rypot[k]=biny;
+ k++;
+
+ }
+ }
+ free(hdu2_loc_flgs) ;
+ free(hdu2_xfarf) ;
+ free(hdu2_yfarf) ;
+ }
+ else {
+ npot=npot2;
+ xpot=xpot2;
+ ypot=ypot2;
+ rxpot=rxpot2;
+ rypot=rypot2;
+ }
+
+ /* allocate space to hold the pothole list (for 1 second of data) */
+ xval = (float *) cf_calloc(npot, sizeof(float)) ;
+ yval = (float *) cf_calloc(npot, sizeof(float)) ;
+ rxval = (float *) cf_calloc(npot, sizeof(float)) ;
+ ryval = (float *) cf_calloc(npot, sizeof(float)) ;
+ chan = (unsigned char *) cf_calloc(npot, sizeof(char)) ;
+ num=-1 ;
+
+
+ cf_verbose(3, "For initial pothole selection, "
+ "we pad extraction windows by %d pixels in Y.", pad);
+
+ /* do each aperture separately */
+
+ for (j=0 ; j<2; j++) {
+
+ /* get the extraction limits for the relevant apertures */
+ npts = cf_extraction_limits(header, ap[j], srctype,
+ &ylow, &yhigh, &xmin, &xmax) ;
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+
+ /* select potholes located within the extraction window */
+ bpndx = (int *) cf_calloc(npot, sizeof(int)) ;
+ npot_sel = -1 ;
+ for (i=0; i<npot; i++) {
+ xndx = (int) (xpot[i] + 0.5) ;
+ bymax = (int) (ypot[i]+rypot[i]+0.5) ;
+ bymin = (int) (ypot[i]-rypot[i]) ;
+ if ( bymax >= ylow[xndx] - pad && bymin <= yhigh[xndx] + pad ) {
+ npot_sel += 1 ;
+ bpndx[npot_sel] = i ;
+ cf_verbose(3, "pothole %d selected: x=%d, y=%d, rx=%d, ry=%d ", i,
+ cf_nint(xpot[i]),cf_nint(ypot[i]),cf_nint(rxpot[i]),cf_nint(rypot[i]));
+ }
+ }
+
+ /* convert npot_sel from an array index to a number of potholes selected */
+ npot_sel++ ;
+ cf_verbose(3, "number of potholes selected for ap %d = %d ",
+ ap[j], npot_sel) ;
+
+ /* put the locations into the output arrays */
+ for (i=0; i<npot_sel; i++) {
+ num+=1 ;
+ ndx = bpndx[i] ;
+ xval[num] = xpot[ndx] ;
+ yval[num] = ypot[ndx] ;
+ chan[num] = ap[j] ;
+ rxval[num] = rxpot[ndx] ;
+ ryval[num] = rypot[ndx] ;
+ }
+ }
+
+ /* convert num from an index to the number of potholes found */
+ num+=1 ;
+ cf_verbose(3, "total number of potholes found = %d ",num) ;
+
+ /* generate a timeline with only the good times included */
+ num_tot = num * ngood ;
+
+ /* copy the non-zero values into the output arrays */
+ time_out = (float *) cf_calloc(num_tot, sizeof(float)) ;
+ weight_out = (float *) cf_calloc(num_tot, sizeof(float)) ;
+ xout = (float *) cf_calloc(num_tot, sizeof(float)) ;
+ yout = (float *) cf_calloc(num_tot, sizeof(float)) ;
+ rxout = (float *) cf_calloc(num_tot, sizeof(float)) ;
+ ryout = (float *) cf_calloc(num_tot, sizeof(float)) ;
+ chan_out = (unsigned char *) cf_calloc(num_tot, sizeof(char)) ;
+ tflag_out = (unsigned char *) cf_calloc(num_tot, sizeof(char)) ;
+ *loc_flag = (unsigned char *) cf_calloc(num_tot, sizeof(char)) ;
+
+ tndx=-1 ;
+ ctot_lif = 0. ;
+ ctot_sic = 0. ;
+ for (j=0; j<ngood; j++) {
+ ndx=good_index[j] ;
+ tval = timeline[ndx] ;
+ ctot_lif += lif_cnt[ndx] ;
+ ctot_sic += sic_cnt[ndx] ;
+ wval_lif = lif_cnt[ndx] ;
+ wval_sic = sic_cnt[ndx] ;
+ tflag_val = statflag[ndx] ;
+ for (i=0; i<num; i++) {
+ tndx += 1;
+ if (tndx > num_tot -1) {
+ cf_verbose(2, "overflow tndx while filling pothole array ") ;
+ tndx= num_tot-1 ; }
+ time_out[tndx] = tval ;
+ if (chan[i] < 4) weight_out[tndx] = wval_lif ;
+ else weight_out[tndx] = wval_sic ;
+ xout[tndx] = xval[i] ;
+ yout[tndx] = yval[i] ;
+ rxout[tndx] = rxval[i] ;
+ ryout[tndx] = ryval[i] ;
+ chan_out[tndx] = chan[i] ;
+ tflag_out[tndx] = tflag_val ;
+ }
+ }
+
+ /* normalize the weights for TTAG data and set to 1 for HIST data*/
+ if (!strncmp(instmode,"T",1))
+ for (i=0; i<num_tot; i++) {
+ if (chan_out[i] < 4) weight_out[i] /= ctot_lif ;
+ else weight_out[i] /= ctot_sic ;
+ }
+ else for (i=0; i<num_tot; i++) weight_out[i] = 1. ;
+
+ /* assign pointers to output variables */
+ *x = xout ;
+ *y = yout ;
+ *rx = rxout ;
+ *ry = ryout ;
+ *channel = chan_out ;
+ *timeflag = tflag_out ;
+ *nevents = num_tot ;
+ *weight = weight_out ;
+ *time = time_out ;
+
+ free(xval) ;
+ free(yval) ;
+ free(rxval) ;
+ free(ryval) ;
+ free(chan) ;
+ free(bpndx) ;
+
+ return 0;
+
+}
+
+
+/*******************************************************************************
+ *
+ * CF_CORRECT_DOPPLER_MOTIONS
+ *
+ * procedure to correct the pothole centroid positions for doppler effects
+ *
+ ******************************************************************************/
+
+ int cf_correct_doppler_motions(fitsfile *header, long nevents, float *time,
+ float *x, unsigned char *channel, long nseconds, float *timeline,
+ float *velocity)
+{
+ char wave_file[FLEN_FILENAME];
+ float *wavelength, v_helio, vel, lam, dlam, dx;
+ fitsfile *wavefits;
+ int status=0, anynull;
+ int fcol, targ_ap[2], ap, ii;
+ long i, k, ndx ;
+
+ /* get the information needed in the analysis */
+
+ /* first determine the target apertures */
+ cf_source_aper(header, targ_ap) ;
+ cf_verbose(3, "target apertures = %d and %d ",targ_ap[0], targ_ap[1]) ;
+
+ /* read in the wavelength calibration */
+ FITS_read_key(header, TSTRING, "WAVE_CAL", wave_file, NULL, &status);
+ FITS_open_file(&wavefits, cf_cal_file(wave_file), READONLY, &status);
+ wavelength = (float *) cf_calloc(NXMAX, sizeof(float));
+
+ /* read the heliocentric velocity from the header */
+ FITS_read_key(header, TFLOAT, "V_HELIO", &v_helio, NULL, &status);
+ cf_verbose(3, "heliocentric velocity = %f ",v_helio) ;
+
+ /* Go through each target aperture */
+ for (i = 0; i < 2; i++) {
+ ap = targ_ap[i] ;
+
+ FITS_movabs_hdu(wavefits, ap+1, NULL, &status);
+ FITS_get_colnum(wavefits, TRUE, "WAVELENGTH", &fcol, &status);
+ fits_read_col(wavefits, TFLOAT, fcol, 1, 1, NXMAX, NULL, wavelength,
+ &anynull, &status);
+
+ ndx = 0 ;
+ for (k = 0; k < nevents; k++)
+ if (channel[k] == ap) {
+ while (time[k] >= timeline[ndx] ) ndx++ ;
+ if (ndx > nseconds-1 ) ndx=nseconds-1 ;
+ ii = (int) (x[k] + 0.5);
+ if (ii >= 1 && ii < NXMAX) {
+ lam = wavelength[ii];
+ vel = velocity[ndx] + v_helio ;
+ dlam = vel * lam / C ;
+ dx = dlam / ( wavelength[ii] - wavelength[ii-1]) ;
+ x[k] += dx ;
+ if (vel > 100 || vel < -100 )
+ cf_if_warning ("at k=%ld, time=%f, lam=%f, vel=%f, "
+ "dlam=%f, dx=%f \n",
+ k, time[k], lam, vel, dlam, dx) ;
+ }
+ }
+
+ }
+ free(wavelength);
+ FITS_close_file(wavefits, &status);
+
+ return status;
+}
+
+
+/*******************************************************************************
+ *
+ * CF_COMBINE_POTHOLE_DATA
+ *
+ * procedure to combine the centroid data as a function of time with the
+ * pothole size information to generate a list of affected pixels on the
+ * detector
+ *
+ ******************************************************************************/
+
+ long cf_combine_pothole_data(fitsfile *header, long nevents, float *time,
+ float *weight, float *x, float *rx, float *y, float *ry,
+ unsigned char *channel, unsigned char *timeflag, float **xpix,
+ float **ypix, unsigned char **chan_pix, float **wt_pix) {
+
+ int srctype, expstat, aperture[2];
+ int npot, xndx, yndx, ap, nfill ;
+ short sxmin, sxmax, *ylow, *yhigh ;
+ int xdim1, ydim1, xmin1, xmax1, ymin1, ymax1 ;
+ int xminp, xmaxp, yminp, ymaxp ;
+ long i, j, k, l, m, num, numt ;
+ long nout, nmax, npts1, npts2, xdim, ydim, ndx, ndxp ;
+ float *xpixt, *ypixt, *wt_pixt, xmin, xmax, ymin, ymax ;
+ float *array1, *array2, wtot, wval, wmax ;
+ float xv2, yv2, rx2, ry2, pos, xval, yval ;
+ float rxp, ryp, xv, yv, time0 ;
+ unsigned char *chan_pixt;
+
+ int status=0;
+ char instmode[FLEN_VALUE];
+ unsigned char TEMPORAL_MASK;
+
+ cf_verbose(3, "Combining pothole data") ;
+
+ /*
+ * Read INSTMODE keyword. If in HIST mode, mask out
+ * LIMB, SAA, and BRST flags. TEMPORAL_DAY is the default.
+ */
+ TEMPORAL_MASK = TEMPORAL_DAY;
+ FITS_read_key(header, TSTRING, "INSTMODE", instmode, NULL, &status);
+ if (!strncmp(instmode, "HIST", 4)) {
+ TEMPORAL_MASK |= TEMPORAL_LIMB;
+ TEMPORAL_MASK |= TEMPORAL_SAA;
+ TEMPORAL_MASK |= TEMPORAL_BRST;
+ }
+ /*
+ * If EXP_STAT = 2, target is bright earth/airglow. Mast out limb-angle flag.
+ */
+ FITS_read_key(header, TINT, "EXP_STAT", &expstat, NULL, &status);
+ if (expstat == 2) TEMPORAL_MASK |= TEMPORAL_LIMB;
+
+ /* Read source type from header. */
+ srctype = cf_source_aper(header, aperture) ;
+
+ /* determine the number of potholes by looking at the number
+ of entries for time=time[0] */
+
+ time0 = time[0] ;
+ i = 1 ;
+ while (fabs((time[i] - time0)) < 0.5) i++ ;
+ cf_verbose(3, "Number of potholes found = %d ", npot = i) ;
+
+ nmax = 0;
+ /* estimate size of the array needed to store the data and open
+ arrays to contain it */
+ for (i=0; i<npot; i++) nmax += (60.+2.*rx[i]) * (30.+2*ry[i]) ;
+ nmax *= 1.5;
+ cf_verbose(3, "Estimated storage space needed for output array = %d ",
+ nmax) ;
+ xpixt = (float *) cf_calloc(nmax, sizeof(float)) ;
+ ypixt = (float *) cf_calloc(nmax, sizeof(float)) ;
+ wt_pixt = (float *) cf_calloc(nmax, sizeof(float)) ;
+ chan_pixt = (unsigned char *) cf_calloc(nmax, sizeof(char)) ;
+
+ /* initialize nout = total number of pixels affected by potholes */
+ nout = -1 ;
+
+ /* process each pothole individually */
+
+ ap = -1;
+ for (i=0; i<npot; i++) {
+
+ nfill = 0;
+
+ /* If channel has changed, get the new extraction limits */
+ if (ap != (int) channel[i]) {
+ ap = (int) channel[i] ;
+ (void) cf_extraction_limits(header, ap, srctype,
+ &ylow, &yhigh, &sxmin, &sxmax);
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+ }
+
+ /* determine the range in X and Y for the pothole centroid positions */
+ xmin = 17000. ;
+ xmax = 0. ;
+ ymin = 1024. ;
+ ymax = 0. ;
+ numt = 0 ;
+ for (j=i; j<nevents ; j+= npot) {
+ numt ++ ;
+ if ((x[j] > xmax) && (x[j] < 16383)) xmax = x[j] ;
+ if ((x[j] < xmin) && (x[j] > 0)) xmin = x[j] ;
+ if ((y[j] > ymax) && (y[j] < 1023)) ymax = y[j] ;
+ if ((y[j] < ymin) && (y[j] > 0)) ymin = y[j] ;
+ }
+ cf_verbose(3,"For pothole %.1d, xmin=%.1d, "
+ "xmax=%.1d, ymin=%.1d, ymax=%.1d",
+ i, cf_nint(xmin), cf_nint(xmax), cf_nint(ymin), cf_nint(ymax)) ;
+
+ /* set up a 2D array which contains all of the movement in the pothole */
+ xdim = (long) (xmax - xmin + 1.5) ;
+ ydim = (long) (ymax - ymin + 1.5) ;
+ npts1 = xdim * ydim ;
+ cf_verbose(3, "Setting up array1: xdim=%d, ydim=%d, npts=%d ",
+ xdim, ydim, npts1) ;
+ array1 = (float *) cf_calloc(npts1, sizeof(float)) ;
+
+ /* Determine the 2D distribution of positions for the pothole */
+ num = 0 ;
+ for (j=i ; j<nevents ; j+= npot)
+ if (!(timeflag[j] & ~TEMPORAL_MASK)) {
+ num++ ;
+ xndx = cf_nint(x[j] - xmin) ;
+ if (xndx < 0 ) xndx = 0 ;
+ if (xndx > xdim-1) xndx = xdim-1 ;
+ yndx = cf_nint(y[j] - ymin) ;
+ if (yndx < 0 ) yndx = 0 ;
+ if (yndx > ydim-1) yndx = ydim-1 ;
+ ndx = xndx + yndx*xdim ;
+ if (ndx > npts1-1) ndx=npts1-1 ;
+ array1[ndx] += weight[j] ;
+ }
+
+ if (i == 0 ) cf_verbose(3, "%d good seconds out of a total of %d ",
+ num, numt) ;
+
+ wtot=0. ;
+ for(j=0; j<ydim; j++)
+ for (k=0; k<xdim; k++) wtot += array1[j*xdim + k] ;
+ cf_verbose(3, "wtot for array 1= %d", cf_nint(wtot)) ;
+
+ /* specify the size of the pothole */
+ rxp = rx[i] ;
+ rx2 = rxp * rxp ;
+ ryp = ry[i] ;
+ ry2 = ryp * ryp ;
+ cf_verbose(3,"pothole size: rx=%d, ry=%d ",cf_nint(rxp),cf_nint(ryp));
+
+ /* set up a target array to contain the image of the full pothole, after
+ reconstruction */
+ xmax1 = xmax + rxp ;
+ xmin1 = xmin - rxp ;
+ xdim1 = (int) (xmax1 - xmin1 + 1.5) ;
+ xndx = (long) x[i] ;
+ ymax1 = ymax + ryp ;
+ ymin1 = ymin - ryp ;
+ cf_verbose(3,"area of the detector covered by pothole :") ;
+ cf_verbose(3, " xmin=%d, xmax=%d, ymin=%d, ymax=%d ",
+ xmin1,xmax1,ymin1,ymax1) ;
+ cf_verbose(3,"extraction window extends from %d to %d ",ylow[xndx], yhigh[xndx]) ;
+
+ /* Include only pixels that lie within extraction window. */
+ if (ymin1 < ylow[xndx]) ymin1 = ylow[xndx];
+ if (ymax1 > yhigh[xndx]) ymax1 = yhigh[xndx];
+ ydim1 = (int) (ymax1 - ymin1 + 1.5) ;
+
+ /* If no pixels overlap the extraction window, skip to the next one. */
+ if (ydim1 <= 0) {
+ cf_verbose(3, "This pothole does not overlap the extraction window. We'll skip it.");
+ continue;
+ }
+
+ npts2 = xdim1 * ydim1 ;
+ cf_verbose(3, "setting up array2: xdim = %d, ydim = %d, npts2 = %d ",
+ xdim1, ydim1, npts2) ;
+ array2 = (float *) cf_calloc(npts2, sizeof(float)) ;
+
+ /* fill in array2
+ - scan all of the positions specified in array1 */
+ for (j=0 ; j<ydim ; j++)
+ for (k=0 ; k<xdim ; k++){
+ /* specify the properties of the pothole at this location */
+ xval = k + xmin ;
+ yval = j + ymin ;
+ ndx = j*xdim + k ;
+ if (ndx > npts1-1) ndx=npts1-1 ;
+ wval = array1[ndx] ;
+
+ if (wval > 0) {
+ /* set limits to the region of array2 affected by the pothole */
+ xminp = (int) (xval - rxp - xmin1 + 0.5) ;
+ xmaxp = (int) (xval + rxp - xmin1 + 0.5) ;
+ yminp = (int) (yval - ryp - ymin1 + 0.5) ;
+ ymaxp = (int) (yval + ryp - ymin1 + 0.5) ;
+ if (yminp < 0) yminp = 0 ;
+ if (ymaxp > ydim1-1) ymaxp = ydim1-1 ;
+ if (xminp < 0) xminp = 0 ;
+ if (xmaxp > xdim1-1) xmaxp = xdim1-1 ;
+
+ /* fill in the array for this pothole location */
+ for (l=yminp; l<=ymaxp; l++)
+ for (m=xminp; m <= xmaxp ; m++) {
+ ndxp = l * xdim1 + m ;
+ if (ndxp > npts2-1) ndxp = npts2 -1 ;
+ xv = m + xmin1 - xval ;
+ yv = l + ymin1 - yval ;
+ xv2= xv * xv ;
+ yv2 = yv * yv ;
+ pos = (xv2/rx2) + (yv2/ry2) ;
+ if ( pos <= 1) array2[ndxp] += wval ;
+ nfill++;
+ }
+ }
+ }
+
+ cf_verbose(3,"number of pixels affected by this pothole = %d ",nfill);
+
+ /* determine the maximum weights across the pothole map */
+ wmax=0. ;
+ for(j=0; j<ydim1; j++)
+ for (k=0; k<xdim1; k++) {
+ ndx = j * xdim1 + k ;
+ if (array2[ndx] > wmax) wmax = array2[ndx] ;
+ }
+ cf_verbose(3, "wmax for array2 = %d ", cf_nint(wmax)) ;
+
+ /* normalize the values of the affected pixels to the maximum and put
+ their locations into the output array */
+ for(j=0; j<ydim1; j++)
+ for (k=0; k<xdim1; k++) {
+ ndx = j * xdim1 + k ;
+ if (array2[ndx] > 0) {
+ if (nout < nmax-1) {
+ nout ++ ;
+ xpixt[nout] = k + xmin1;
+ ypixt[nout] = j + ymin1;
+ wt_pixt[nout] = array2[ndx]/wmax ;
+ chan_pixt[nout] = channel[i] ;
+ }
+ else {
+ cf_if_warning("Array overflow. Truncating pseudo-photon list.");
+ break;
+ }
+ }
+ }
+
+
+ cf_verbose(3, "%d points tabulated after processing pothole %d\n",nout,i);
+
+ free(array1) ;
+ free(array2) ;
+ }
+
+ /* assign output pointers */
+ *xpix = xpixt ;
+ *ypix = ypixt ;
+ *wt_pix = wt_pixt ;
+ *chan_pix = chan_pixt ;
+
+ cf_verbose(3, "Total number of pixels tabulated = %d ", nout) ;
+
+ return nout ;
+}
+
+
+/****************************************************************************/
+
+int
+main(int argc, char *argv[])
+{
+ unsigned char *timeflag=NULL, *statflag=NULL, *locflag=NULL;
+ unsigned char *loc_flag=NULL, *channel=NULL ;
+ int grating_motion=1, mirror_motion=1, fpa_pos=1;
+ int jitter=1, doppler_motion=1, astig=1, tscreen=1, night_only=FALSE;
+ int status=0, hdutype, optc;
+ long nevents, nseconds, i, nrows=1, npixels ;
+ long ngood, *good_index, dtime, ntime ;
+ float exptime, *time=NULL, *x=NULL, *y=NULL, *velocity=NULL;
+ float *rx=NULL, *ry=NULL, *xpix, *ypix, *wt_pix ;
+ float *lif_cnt=NULL, *sic_cnt=NULL, *lambda=NULL ;
+ float *timeline=NULL, *weight=NULL;
+ short *tsunrise=NULL, *tsunset=NULL ;
+ fitsfile *header, *bpmfits, *memp;
+ char rootname[FLEN_VALUE]={'\0'}, bpm_file[FLEN_FILENAME]={'\0'};
+ char idf_file[FLEN_FILENAME]={'\0'}, det[FLEN_VALUE]={'\0'};
+ char *keywords[] = {"GRAT_COR","FPA__COR","MIRR_COR","ASTG_COR","WAVE_COR"};
+ char keyval[FLEN_VALUE];
+ unsigned char *chan_pix ;
+
+ char extname[]="POTHOLE_DATA", fmt_float[FLEN_VALUE], fmt_byte[FLEN_VALUE] ;
+ char instmode[FLEN_VALUE] ;
+ int tfields=5 ;
+ char *ttype[]={ "X", "Y", "CHANNEL", "WEIGHT", "LAMBDA"} ;
+ char *tform[5] ; /* we will define the tform once the number of
+ array elements is known */
+ char *tunit[]={ "pixels", "pixels", "unitless", "unitless", "Angstroms"} ;
+
+ char opts[] = "hgmfjdasn:v:";
+ char usage[] =
+ "Usage:\n"
+ " cf_bad_pixels [-hvgmfjdas] [-n bpm_filename] [-v level] idffile\n";
+ char option[] =
+ "Options:\n"
+ " -h: this help message\n"
+ " -v: verbosity level (=1; 0 is silent)\n"
+ " -g: no grating motion correction \n"
+ " -m: no mirror motion correction \n"
+ " -n: Create BPM file for nighttime fraction of exposure.\n"
+ " Argument is the name of the output BPM file.\n"
+ " -f: no fpa motion correction \n"
+ " -j: no jitter correction \n"
+ " -d: no doppler correction \n"
+ " -a: no astigmatism correction \n"
+ " -s: no screening on time flags (if 0)" ;
+
+ verbose_level = 1;
+
+ while ((optc = getopt(argc, argv, opts)) != -1) {
+ switch(optc) {
+
+ case 'h':
+ printf("%s\n%s", usage, option);
+ return 0;
+ case 'v':
+ verbose_level = atoi(optarg);
+ break;
+ case 'g':
+ grating_motion = 0 ;
+ break;
+ case 'm':
+ mirror_motion = 0 ;
+ break ;
+ case 'n':
+ sprintf(bpm_file,"%s", optarg) ;
+ night_only = TRUE;
+ break;
+ case 'f':
+ fpa_pos = 0 ;
+ break ;
+ case 'j':
+ jitter = 0 ;
+ break ;
+ case 'a':
+ astig = 0 ;
+ break;
+ case 'd':
+ doppler_motion = 0 ;
+ 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 the input IDF file. */
+ FITS_open_file(&header, argv[optind], READWRITE, &status);
+ FITS_read_key(header, TSTRING, "INSTMODE", instmode, NULL, &status) ;
+
+ /* If EXPTIME < 1, exit without generating a bad-pixel file. */
+ FITS_read_key(header, TFLOAT, "EXPTIME", &exptime, NULL, &status);
+ if (night_only)
+ FITS_read_key(header, TFLOAT, "EXPNIGHT", &exptime, NULL, &status);
+ if (exptime < 1.) {
+ cf_verbose(1, "EXPTIME = %g. "
+ "Will not attempt to generate bad-pixel map.", exptime);
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done processing");
+ return 0;
+ }
+
+ /* Read timeline data from the input file */
+ FITS_movabs_hdu(header, 4, &hdutype, &status);
+ nseconds = cf_read_col(header, TFLOAT, "TIME", (void **) &timeline);
+ nseconds = cf_read_col(header, TBYTE, "STATUS_FLAGS", (void **) &statflag);
+ nseconds = cf_read_col(header, TSHORT, "TIME_SUNRISE", (void **) &tsunrise);
+ nseconds = cf_read_col(header, TSHORT, "TIME_SUNSET", (void **) &tsunset) ;
+ nseconds = cf_read_col(header, TFLOAT, "ORBITAL_VEL", (void **) &velocity);
+ nseconds = cf_read_col(header, TFLOAT, "LIF_CNT_RATE", (void **) &lif_cnt) ;
+ nseconds = cf_read_col(header, TFLOAT, "SIC_CNT_RATE", (void **) &sic_cnt) ;
+ FITS_movabs_hdu(header, 1, NULL, &status);
+
+ /*
+ * If night-only spectrum is requested from the command line, copy
+ * IDF header into memory, close the IDF, modify copy, and pass it
+ * to subsequent routines.
+ */
+ if (night_only) {
+ cf_verbose(3, "Night-only spectrum requested from command line.");
+ FITS_create_file(&memp, "mem://", &status);
+ FITS_copy_hdu(header, memp, 0, &status);
+ if (!strncmp(instmode, "H",1)) {
+ FITS_movabs_hdu(header, 2, NULL, &status);
+ FITS_copy_hdu(header, memp, 0, &status);
+ }
+ FITS_close_file(header, &status);
+ header = memp;
+
+ FITS_movabs_hdu(header, 1, NULL, &status);
+ FITS_update_key(header, TFLOAT, "EXPTIME", &exptime, NULL, &status);
+ FITS_update_key(header, TSTRING, "DAYNIGHT", "NIGHT", NULL, &status);
+ }
+ else { /* Generate name of output file */
+ FITS_read_key(header, TSTRING, "ROOTNAME", rootname, NULL, &status) ;
+ FITS_read_key(header, TSTRING, "DETECTOR", det, NULL, &status) ;
+ det[1]=tolower(det[1]) ;
+ if (!strncmp(instmode, "T",1))
+ sprintf(bpm_file,"%11s%2sttagfbpm.fit",rootname,det) ;
+ else
+ sprintf(bpm_file,"%11s%2shistfbpm.fit",rootname,det) ;
+ }
+
+ /* Write name of BPM file to header of IDF. */
+ cf_verbose(3,"Output bad pixel map to file %s", bpm_file) ;
+ FITS_update_key(header, TSTRING, "BPM_CAL", bpm_file, NULL, &status);
+
+ /* Create the output file and copy header from the input */
+ FITS_create_file(&bpmfits, bpm_file, &status) ;
+ FITS_copy_hdu(header, bpmfits, 0, &status) ;
+
+ FITS_update_key(bpmfits, TSTRING, "FILENAME", bpm_file,
+ NULL, &status) ;
+ FITS_update_key(bpmfits, TSTRING, "FILETYPE", "BAD PIXEL MAP",
+ NULL, &status) ;
+
+ FITS_read_key(header, TSTRING, "FILENAME", idf_file, NULL, &status) ;
+ FITS_update_key(bpmfits, TSTRING, "IDF_FILE", idf_file,
+ NULL, &status) ;
+
+ /* Reset analysis keywords in the top level of the output file */
+ FITS_movabs_hdu(bpmfits, 1, &hdutype, &status);
+ for (i = 0; i < 5; i++) {
+ FITS_read_key(bpmfits, TSTRING, keywords[i], keyval, NULL, &status);
+ if (!strncmp(keyval, "COMPLETE", 8))
+ FITS_update_key(bpmfits, TSTRING, keywords[i], "PERFORM", NULL,
+ &status) ;
+ else
+ FITS_update_key(bpmfits, TSTRING, keywords[i], "OMIT", NULL,
+ &status) ;
+ }
+
+ /* Make sure that the count rates are non-zero */
+ for (i=0; i<nseconds; i++) {
+ if (lif_cnt[i] < 0.01) lif_cnt[i] = 0.01 ;
+ if (sic_cnt[i] < 0.01) sic_cnt[i] = 0.01 ;
+ }
+
+ /* generate a null array for locflag */
+ locflag = (unsigned char *) cf_calloc(nseconds, sizeof(unsigned char)) ;
+
+ /* select only the good times */
+ FITS_movabs_hdu(header, 1, &hdutype, &status);
+ cf_apply_filters(header, tscreen, nseconds, statflag, locflag, nseconds,
+ statflag, &dtime, &ntime, &ngood, &good_index) ;
+
+ free(locflag) ;
+
+ cf_verbose(3, "nseconds=%d, dtime=%d, ntime=%d, ngood=%d ",
+ nseconds, dtime, ntime, ngood) ;
+
+ /* Generate the pseudo photons */
+ cf_generate_pseudo_photons(header, ngood, good_index, timeline, statflag,
+ lif_cnt, sic_cnt, &nevents, &time, &weight, &x,
+ &y, &channel, &rx, &ry, &timeflag, &loc_flag) ;
+
+ if (night_only)
+ FITS_delete_file(header, &status);
+ else
+ FITS_close_file(header, &status);
+
+ cf_verbose(3, "number of pseudo photon events = %d ",nevents) ;
+
+ /* Call routines to remove motions */
+ if (grating_motion) {
+ cf_verbose(4,"Correcting for grating motions") ;
+ cf_grating_motion(bpmfits, nevents, time, x, y, channel,
+ nseconds, timeline, tsunrise) ; }
+ if (fpa_pos) {
+ cf_verbose(4,"Correcting for fpa motions ") ;
+ cf_fpa_position(bpmfits, nevents, x, channel) ;}
+ if (mirror_motion) {
+ cf_verbose(4,"Correcting for mirror motions") ;
+ cf_mirror_motion(bpmfits, nevents, time, x, y, channel,
+ nseconds, timeline, tsunset) ; }
+ if (jitter) {
+ cf_verbose(4,"Correcting for satellite jitter ") ;
+ cf_satellite_jitter(bpmfits, nevents, time, x, y, channel, nseconds,
+ timeline, statflag) ;
+ }
+
+ /* correct for the Doppler motions */
+ if (doppler_motion)
+ cf_correct_doppler_motions(bpmfits, nevents, time, x, channel,
+ nseconds, timeline, velocity) ;
+
+ /* combine the pothole centroid information and generate a pothole map by
+ including the pothole sizes and shapes */
+ npixels = cf_combine_pothole_data(bpmfits, nevents, time, weight,
+ x, rx, y, ry, channel, timeflag, &xpix, &ypix, &chan_pix, &wt_pix) ;
+
+ /* generate a blank wavelength array */
+ lambda = (float *) cf_calloc( npixels, sizeof(float) ) ;
+
+ /* correct for astigmatism and assign wavelengths */
+ if (astig) cf_astigmatism(bpmfits, npixels, xpix, ypix, chan_pix) ;
+ cf_dispersion(bpmfits, npixels, xpix, chan_pix, lambda);
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+
+ /* create a table and output the results to the file */
+ /* first set tform values */
+ cf_verbose(3, "Writing pseudo-pixel information to output file.");
+ sprintf(fmt_float, "%1ldE", npixels) ;
+ sprintf(fmt_byte, "%1ldB", npixels) ;
+ tform[0]=fmt_float ;
+ tform[1]=fmt_float ;
+ tform[2]=fmt_byte ;
+ tform[3]=fmt_float ;
+ tform[4]=fmt_float ;
+
+ /* put the data into the table */
+ FITS_create_tbl(bpmfits, BINARY_TBL, nrows, tfields, ttype, tform,
+ tunit, extname, &status) ;
+ cf_verbose(3, " Photon list created.");
+ FITS_write_col(bpmfits, TFLOAT, 1, 1L, 1L, npixels, xpix, &status) ;
+ cf_verbose(3, " X array has been written.");
+ FITS_write_col(bpmfits, TFLOAT, 2, 1L, 1L, npixels, ypix, &status) ;
+ cf_verbose(3, " Y array has been written.");
+ FITS_write_col(bpmfits, TBYTE, 3, 1L, 1L, npixels, chan_pix, &status) ;
+ cf_verbose(3, " Channel array has been written.");
+ FITS_write_col(bpmfits, TFLOAT, 4, 1L, 1L, npixels, wt_pix, &status) ;
+ cf_verbose(3, " Weights array has been written.");
+ FITS_write_col(bpmfits, TFLOAT, 5, 1L, 1L, npixels, lambda, &status) ;
+ cf_verbose(3, " Lambda array has been written.");
+
+ FITS_close_file(bpmfits, &status);
+
+ free(statflag);
+ free(timeline);
+ free(timeflag);
+ free(rx) ;
+ free(ry) ;
+ free(y);
+ free(x);
+ free(time);
+ free(xpix) ;
+ free(ypix) ;
+ free(chan_pix) ;
+ free(wt_pix) ;
+ free(lambda) ;
+
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done processing");
+ return 0;
+}
diff --git a/src/fuv/cf_convert_to_farf.c b/src/fuv/cf_convert_to_farf.c
new file mode 100644
index 0000000..a50dd2e
--- /dev/null
+++ b/src/fuv/cf_convert_to_farf.c
@@ -0,0 +1,268 @@
+/*****************************************************************************
+ * Johns Hopkins University
+ * Center For Astrophysical Sciences
+ * FUSE
+ *****************************************************************************
+ *
+ * Synopsis: cf_convert_to_farf options intermediate_file
+ *
+ * Description: Corrects the RAW data for IDS and electronic deadtime
+ * and for thermal, count rate, geometric, and PHA distortions.
+ * The resulting data are in the Flight Aligned Reference
+ * Frame (FARF), the location of an event on an ideal
+ * detector. The input file is modified in place.
+ *
+ * By default all corrections are performed. Command line
+ * options allow one or more corrections to be omitted.
+ *
+ * Arguments: input_file Intermediate data file
+ *
+ * Calibration files:
+ *
+ * Returns: 0 on successful completion
+ *
+ * Calls: cf_ids_dead_time, cf_electronics_dead_time,
+ * cf_thermal_distort, cf_count_rate_y_distort,
+ * cf_geometric_distort, cf_pha_x_distort
+ *
+ * History: 08/08/02 peb 1.1 Begin work
+ * 10/27/02 peb 1.3 Added IDS and electronics dead time
+ * corrections and time-dependent
+ * correction for count_rate_y_distort.
+ * Changed the help option to -h.
+ * Corrected function description.
+ * 11/12/02 peb 1.4 Added argument to cf_thermal_distort,
+ * cf_count_rate_y_distort, and
+ * cf_geometric_distort calls, so
+ * stim-pulse flag can be set.
+ * Added function to flag events in
+ * inactive regions.
+ * 12/09/02 wvd 1.5 Set keyword TOT_DEAD to mean value
+ * for exposure.
+ * 01/28/03 wvd 1.6 Added cf_check_digitizer
+ * 03/10/03 peb 1.7 Added verbose_level option and unistd.h
+ * so getopt.h will be portable.
+ * 03/10/03 peb 1.8 Changed locflgs to unsigned char *
+ * 03/10/03 peb 1.9 Changed pha to unsigned char *
+ * 06/11/03 wvd 1.10 Pass datatype to cf_read_col and
+ * cf_write_col.
+ * 06/16/03 rdr 1.11 Correct datatype in reading LOC_FLGS
+ * 08/01/03 wvd 1.12 Add cf_apply_dead_time to properly
+ * correct both TTAG and HIST data.
+ * 08/04/03 wvd 1.13 Convert count-rate arrays to shorts.
+ * 08/25/03 wvd 1.14 Change coltype from string to int in
+ * cf_read_col and cf_write_col.
+ * 11/26/03 wvd 1.15 Read aic_rate and fec_rate as floats.
+ * 02/09/04 wvd 1.16 Don't pass aic_rate to
+ * cf_apply_dead_time().
+ * 02/27/04 rdr 1.17 Added weight to cf_thermal_distortion
+ * 03/02/05 wvd 1.18 Must assign channel numbers and pulse-
+ * height values to HIST data before walk
+ * correction.
+ * Read XFARF and YFARF from XRAW and YRAW.
+ * 05/20/05 wvd 1.19 Clean up i/o.
+ * 11/02/06 wvd 1.20 Add call to cf_time_xy_distort.
+ * 12/29/06 wvd 1.21 Move cf_screen_fifo_overflow from
+ * cf_screen_photons. Rename it
+ * cf_fifo_dead_time. Set all dead-time
+ * keywords to 1.0 by default.
+ * 02/08/08 wvd 1.22 For HIST data, write modified PHA
+ * values to the IDF.
+ *
+ ****************************************************************************/
+
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+#include "calfuse.h"
+
+int
+main(int argc, char *argv[])
+{
+ char CF_PRGM_ID[] = "cf_convert_to_farf";
+ char CF_VER_NUM[] = "1.22";
+
+ unsigned char *pha=NULL, *locflags=NULL;
+ char instmode[FLEN_VALUE], aperture[FLEN_VALUE];
+ int check_digitizer=1, ids_deadtime=1, elec_deadtime=1, fifo_deadtime=1;
+ int thermal_corr=1;
+ int countrate_corr=1, geometric_corr=1, pha_corr=1, active_flag=1;
+ int temporal_corr=1, status=0, hdutype, last_call, optc, pad;
+ long i, nevents, nseconds;
+ float *elec_dtc=NULL, *ids_dtc=NULL, unity = 1.0;
+ float *time=NULL, *xfarf=NULL, *yfarf=NULL, *weight=NULL;
+ float *timeline=NULL;
+ float *aic_rate=NULL, *fec_rate=NULL;
+ fitsfile *header;
+
+ char opts[] = "hceiftrsgpav:";
+ char usage[] =
+ "Usage:\n"
+ " cf_convert_to_farf [-hceiftrsgpa] [-v level] idf_file\n";
+ char option[] =
+ "Options:\n"
+ " -h: this help message\n"
+ " -v: verbosity level (default is 1; 0 is silent)\n"
+ " -c: no check of digitizer values\n"
+ " -e: no electronics deadtime correction\n"
+ " -i: no IDS deadtime correction\n"
+ " -f: no FIFO deadtime correction\n"
+ " -t: no thermal distortion correction\n"
+ " -r: no count-rate distortion correction\n"
+ " -s: no time-dependent distortion correction\n"
+ " -g: no geometric distortion correction\n"
+ " -p: no PHA distortion correction\n"
+ " -a: do not flag events beyond active region\n";
+
+ verbose_level=1;
+
+ while ((optc = getopt(argc, argv, opts)) != -1) {
+ switch(optc) {
+
+ case 'h':
+ printf("%s\n%s", usage, option);
+ return 0;
+ case 'v':
+ verbose_level = atoi(optarg);
+ break;
+ case 'c':
+ check_digitizer = 0;
+ break;
+ case 'i':
+ ids_deadtime = 0;
+ break;
+ case 'e':
+ elec_deadtime = 0;
+ break;
+ case 'f':
+ fifo_deadtime = 0;
+ break;
+ case 't':
+ thermal_corr = 0;
+ break;
+ case 'r':
+ countrate_corr = 0;
+ break;
+ case 's':
+ temporal_corr = 0;
+ break;
+ case 'g':
+ geometric_corr = 0;
+ break;
+ case 'p':
+ pha_corr = 0;
+ break;
+ case 'a':
+ active_flag = 0;
+ break;
+ }
+ }
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+ if (argc <= optind) {
+ printf("%s", usage);
+ cf_if_error("Incorrect number of program arguments");
+ }
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");
+
+ FITS_open_file(&header, argv[optind], READWRITE, &status);
+
+ /* Read a couple of keywords from the file header. */
+ FITS_read_key(header, TSTRING, "INSTMODE", instmode, NULL, &status);
+ FITS_read_key(header, TSTRING, "APERTURE", aperture, NULL, &status);
+
+ FITS_movabs_hdu(header, 2, &hdutype, &status);
+ nevents = cf_read_col(header, TFLOAT, "TIME", (void **) &time);
+ nevents = cf_read_col(header, TBYTE, "PHA", (void **) &pha);
+ nevents = cf_read_col(header, TFLOAT, "XRAW", (void **) &xfarf);
+ nevents = cf_read_col(header, TFLOAT, "YRAW", (void **) &yfarf);
+ nevents = cf_read_col(header, TFLOAT, "WEIGHT", (void **) &weight);
+ nevents = cf_read_col(header, TBYTE, "LOC_FLGS", (void **) &locflags);
+
+ FITS_movabs_hdu(header, 4, &hdutype, &status);
+ nseconds = cf_read_col(header, TFLOAT, "TIME", (void **) &timeline);
+ nseconds = cf_read_col(header, TFLOAT, "FEC_CNT_RATE", (void **) &fec_rate);
+ nseconds = cf_read_col(header, TFLOAT, "AIC_CNT_RATE", (void **) &aic_rate);
+ FITS_movabs_hdu(header, 1, &hdutype, &status);
+
+ /* Initialize dead-time correction arrays and header keywords. */
+ FITS_update_key (header, TFLOAT, "DET_DEAD", &unity, NULL, &status);
+ FITS_update_key (header, TFLOAT, "IDS_DEAD", &unity, NULL, &status);
+ FITS_update_key (header, TFLOAT, "TOT_DEAD", &unity, NULL, &status);
+ elec_dtc = (float *) cf_malloc(sizeof(float) * nseconds);
+ ids_dtc = (float *) cf_malloc(sizeof(float) * nseconds);
+ for (i = 0; i < nseconds; i++) elec_dtc[i] = ids_dtc[i] = 1.0;
+
+ if (check_digitizer)
+ cf_check_digitizer(header);
+
+ if (elec_deadtime)
+ cf_electronics_dead_time(header, nseconds, fec_rate, elec_dtc);
+
+ if (ids_deadtime)
+ cf_ids_dead_time(header, nseconds, aic_rate, ids_dtc);
+
+ if (fifo_deadtime)
+ cf_fifo_dead_time(header, nevents, time,
+ nseconds, timeline, aic_rate, ids_dtc);
+
+ if (ids_deadtime || elec_deadtime || fifo_deadtime)
+ cf_apply_dead_time(header, nevents, time, weight,
+ nseconds, timeline, ids_dtc, elec_dtc);
+
+ if (thermal_corr)
+ cf_thermal_distort(header, nevents, xfarf, yfarf, weight, locflags);
+
+ if (countrate_corr)
+ cf_count_rate_y_distort(header, nevents, time, yfarf, locflags,
+ nseconds, timeline, fec_rate);
+ if (temporal_corr)
+ cf_time_xy_distort(header, nevents, xfarf, yfarf, locflags);
+
+ if (geometric_corr)
+ cf_geometric_distort(header, nevents, xfarf, yfarf, locflags);
+
+ /* Must assign PHA values to HIST data before applying walk correction. */
+ if (pha_corr && !strncmp(instmode, "HIST", 4)) {
+ unsigned char *channel;
+ channel = cf_calloc(nevents, sizeof(unsigned char));
+ if (!strncmp(aperture, "RFPT", 4)) pad = 10;
+ else pad = 100;
+ cf_identify_channel(header, nevents, xfarf, yfarf, channel, locflags,
+ pad, last_call=FALSE);
+ cf_modify_hist_pha(header, nevents, pha, channel);
+ free(channel);
+ }
+
+ if (pha_corr)
+ cf_pha_x_distort(header, nevents, pha, xfarf, locflags);
+
+ if (active_flag)
+ cf_active_region(header, nevents, xfarf, yfarf, locflags);
+
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+
+ FITS_movabs_hdu(header, 2, &hdutype, &status);
+ if (pha_corr && !strncmp(instmode, "HIST", 4))
+ cf_write_col(header, TBYTE, "PHA", (void *) pha, nevents);
+ cf_write_col(header, TFLOAT, "WEIGHT", (void *) weight, nevents);
+ cf_write_col(header, TFLOAT, "XFARF", (void *) xfarf, nevents);
+ cf_write_col(header, TFLOAT, "YFARF", (void *) yfarf, nevents);
+ cf_write_col(header, TBYTE, "LOC_FLGS", (void *) locflags, nevents);
+
+ FITS_close_file(header, &status);
+
+ free(fec_rate);
+ free(aic_rate);
+ free(timeline);
+ free(locflags);
+ free(weight);
+ free(yfarf);
+ free(xfarf);
+ free(pha);
+ free(time);
+ free(elec_dtc);
+ free(ids_dtc);
+
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done processing");
+ return 0;
+}
diff --git a/src/fuv/cf_countmap.c b/src/fuv/cf_countmap.c
new file mode 100644
index 0000000..2d62b14
--- /dev/null
+++ b/src/fuv/cf_countmap.c
@@ -0,0 +1,541 @@
+/*****************************************************************************
+ * Johns Hopkins University
+ * Center For Astrophysical Sciences
+ * FUSE
+ *****************************************************************************
+ *
+ * Synopsis: cf_countmap input_files
+ *
+ * Description: Adds the counts from a list of input idf files into
+ * the count map. The accumulated count map cube is
+ * determined by reading the COUNT_CAL keyword in the
+ * input file's header. To prevent multiple copies
+ * of reprocessed input files from being added to
+ * the total exposure map, if the input file has been
+ * previously added, it is skipped over and a message is
+ * printed out.
+ *
+ *
+ * Arguments: input_files Input initialized ttag FITS file names
+ *
+ * Calibration files required: None
+ *
+ * Returns: int 0 = successful completion
+ *
+ * Note: This routine requires significantly more memory than your
+ * average routine. Its design has not been optimized to
+ * reduce memory usage.
+ *
+ * History: 03/25/03 1.1 peb Begin and finish work. Based on
+ * cf_ttag_to_image
+ * 03/31/03 1.2 peb Removed filtering of photons and fixed
+ * exposure time bug when subtracting data
+ * 06/11/03 1.3 wvd Pass datatype to cf_read_col
+ * 08/25/03 1.4 wvd Change coltype from string to int in
+ * cf_read_col.
+ * 08/25/03 1.5 wvd Change CHIDCOR1 to CHID_COR
+ * 02/17/04 1.6 bjg Added scaling option
+ * Routine name changed from
+ * cf_ttag_countmap to cf_countmap.
+ * 02/26/04 1.7 bjg Added option to select coordinates
+ * to use: RAW, FARF (default), final
+ * 04/08/04 1.8 bjg Remove unused function get_xy_columns
+ * 04/27/04 1.9 bjg Change author name of countfile
+ * 08/25/04 1.10 bjg Added lower and upper limits option
+ * for PHA. In force mode, issue a warning
+ * if file already included.
+ * 03/02/05 1.11 wvd Don't write AUTHOR keyword to header.
+ * 08/24/07 1.12 bot Changed nx and ny to long and called
+ * as TLONG ; changed binx and biny to int
+ * and called as TINT.
+ *
+ ****************************************************************************/
+
+#include <unistd.h>
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "calfuse.h"
+
+#define ROOT_LENGTH 13
+#define MAXCHARS 120
+#define MASTER_CAL_FILE "master_calib_file.dat"
+#define IMAGE_EXT 1
+#define INDEX_EXT 2
+
+static char CF_PRGM_ID[] = "cf_countmap";
+static char CF_VER_NUM[] = "1.12";
+
+
+static int
+new_imagfile(char *newfile, char *detector, float effmjd, int binx, int biny)
+{
+ char date[FLEN_CARD]={'\0'};
+ char *ttype[] = {"ROOTNAME", "OBSDATE", "COUNTS", "EXPTIME", "ADDDATE", "SCALE", "PHAMIN", "PHAMAX"};
+ char *tform[] = {"13A" , "19A" , "1J" , "1E" , "19A" , "1E" , "1I" , "1I" };
+ char *tunit[] = {" " , " " , "counts", "seconds", " " , " " , " " , " " };
+ int status=0, version=1, tref = 0, tfields = 8;
+ long axis[2];
+ float exptime=0.;
+ double nevents=0.;
+ fitsfile *imagfits;
+
+ FITS_create_file(&imagfits, newfile, &status);
+
+ axis[0] = NXMAX/binx;
+ axis[1] = NYMAX/biny;
+
+ FITS_create_img(imagfits, FLOAT_IMG, 2, axis, &status);
+ FITS_write_key(imagfits, TINT, "SPECBINX", &binx,
+ "Binsize in detector X coordinate", &status);
+ FITS_write_key(imagfits, TINT, "SPECBINY", &biny,
+ "Binsize in detector Y coordinate", &status);
+ FITS_write_key(imagfits, TSTRING, "CALFTYPE", "COUNT",
+ "Calibration file type", &status);
+ FITS_write_key(imagfits, TINT, "CALFVERS", &version,
+ "Calibration file version", &status);
+ FITS_write_key(imagfits, TSTRING, "DETECTOR", detector,
+ "detector (1A, 1B, 2A, 2B)", &status);
+ FITS_write_key(imagfits, TFLOAT, "EFFMJD", &effmjd,
+ "Date on which file should be applied (MJD)", &status);
+ fits_get_system_time(date, &tref, &status);
+ FITS_write_key(imagfits, TSTRING, "DATE", date,
+ "file creation date (YYYY-MM-DDThh:mm:ss UTC)", &status);
+ FITS_write_key(imagfits, TDOUBLE, "NEVENTS", &nevents,
+ "Accumulated events", &status);
+ FITS_write_key(imagfits, TFLOAT, "EXPTIME", &exptime,
+ "Accumulated exposure time", &status);
+
+ FITS_create_tbl(imagfits, BINARY_TBL, 0, tfields, ttype, tform, tunit,
+ "PROCESSED FILES", &status);
+
+ FITS_close_file(imagfits, &status);
+ return status;
+}
+
+static int
+_check_index(fitsfile *imagfits, char *rootname, float *scale, int *phamin, int *phamax)
+{
+ int status=0, anynull=0, hdutype;
+ float *scal;
+ int *phami, *phama;
+ long j, nrows, frow=1, felem=1;
+ char **filename, strnull[] = " ";
+
+ FITS_movabs_hdu(imagfits, INDEX_EXT, &hdutype, &status);
+ FITS_read_key(imagfits, TLONG, "NAXIS2", &nrows, NULL, &status);
+
+ filename = (char **) cf_malloc(nrows * sizeof(char *));
+ filename[0] = (char *) cf_malloc(nrows*(ROOT_LENGTH+1)*sizeof(char));
+ scal = (float *) cf_malloc(nrows * sizeof(float));
+ phami = (int *) cf_malloc(nrows * sizeof(int));
+ phama = (int *) cf_malloc(nrows * sizeof(int));
+ for (j=1; j<nrows; j++)
+ filename[j] = filename[j-1]+(ROOT_LENGTH+1);
+
+ FITS_read_col(imagfits, TSTRING, 1, frow, felem, nrows, strnull,
+ filename, &anynull, &status);
+ FITS_read_col(imagfits, TFLOAT, 6, frow, felem, nrows, strnull,
+ scal, &anynull, &status);
+ FITS_read_col(imagfits, TINT, 7, frow, felem, nrows, strnull,
+ phami, &anynull, &status);
+ FITS_read_col(imagfits, TINT, 8, frow, felem, nrows, strnull,
+ phama, &anynull, &status);
+
+ for (j=0; j<nrows; j++) {
+ if (strncmp(filename[j], rootname, strlen(rootname)) == 0) {
+ *scale = scal[j];
+ *phamin = phami[j];
+ *phamax = phama[j];
+ free(filename[0]);
+ free(filename);
+ free(scal);
+ return j+1;
+ }
+ }
+ free(filename[0]);
+ free(filename);
+ free(scal);
+ free(phami);
+ free(phama);
+ return 0;
+}
+
+static int
+_update_index(fitsfile *imagfits, char *rootname, char *obsdate,
+ double nevents, float exptime, float scale, int phamin, int phamax)
+{
+ char *strptr[1], adddate[FLEN_CARD]={'\0'};
+ int status=0, hdutype, tref=0;
+ long nrows;
+
+ FITS_movabs_hdu(imagfits, INDEX_EXT, &hdutype, &status);
+ FITS_read_key(imagfits, TLONG, "NAXIS2", &nrows, NULL, &status);
+
+ strptr[0] = rootname;
+ FITS_write_col(imagfits, TSTRING, 1, nrows+1, 1, 1, strptr, &status);
+ strptr[0] = obsdate;
+ FITS_write_col(imagfits, TSTRING, 2, nrows+1, 1, 1, strptr, &status);
+ FITS_write_col(imagfits, TDOUBLE, 3, nrows+1, 1, 1, &nevents, &status);
+ FITS_write_col(imagfits, TFLOAT, 4, nrows+1, 1, 1, &exptime, &status);
+ fits_get_system_time(adddate, &tref, &status);
+ strptr[0] = adddate;
+ fits_write_col(imagfits, TSTRING, 5, nrows+1, 1, 1, strptr, &status);
+ FITS_write_col(imagfits, TFLOAT, 6, nrows+1, 1, 1, &scale, &status);
+ FITS_write_col(imagfits, TINT, 7, nrows+1, 1, 1, &phamin, &status);
+ FITS_write_col(imagfits, TINT, 8, nrows+1, 1, 1, &phamax, &status);
+ return 0;
+}
+
+
+
+static double
+add_events(fitsfile *infits, char *xcolumn, char *ycolumn, float scale,
+ int phamin, int phamax, fitsfile *imagfits)
+{
+ int binx, biny, status=0, hdutype, anynull=0;
+ long nx, ny, nevents, j;
+ float *xpos, *ypos, *weight, *image;
+ char *pha;
+ double ncounts=0.;
+
+ FITS_movabs_hdu(infits, 2, &hdutype, &status);
+ nevents = cf_read_col(infits, TFLOAT, xcolumn, (void **) &xpos);
+ nevents = cf_read_col(infits, TFLOAT, ycolumn, (void **) &ypos);
+ nevents = cf_read_col(infits, TFLOAT, "WEIGHT", (void **) &weight);
+ nevents = cf_read_col(infits, TBYTE, "PHA", (void **) &pha);
+
+ FITS_movabs_hdu(imagfits, IMAGE_EXT, &hdutype, &status);
+ FITS_read_key(imagfits, TLONG, "NAXIS1", &nx, NULL, &status);
+ FITS_read_key(imagfits, TLONG, "NAXIS2", &ny, NULL, &status);
+ FITS_read_key(imagfits, TINT, "SPECBINX", &binx, NULL, &status);
+ FITS_read_key(imagfits, TINT, "SPECBINY", &biny, NULL, &status);
+
+ image = (float *) cf_malloc(sizeof(float)*nx*ny);
+ FITS_read_img(imagfits, TFLOAT, 1, nx*ny, NULL, image, &anynull, &status);
+
+ for (j=0; j<nevents; j++) {
+ if ((pha[j]>=phamin) && (pha[j]<=phamax)){
+ int x, y;
+ x = xpos[j];
+ y = ypos[j];
+ if ((x >= 0) && (x < nx*binx) && (y >= 0) && (y < ny*biny)) {
+ image[(y/biny)*nx + (x/binx)] += weight[j]*scale;
+ ncounts += weight[j]*scale;
+ }
+ }
+ }
+ FITS_write_img(imagfits, TFLOAT, 1, nx*ny, image, &status);
+
+ free(image);
+ free(weight);
+ free(ypos);
+ free(xpos);
+ return ncounts;
+}
+
+static double
+subtract_events(fitsfile *infits, char *xcolumn, char *ycolumn,
+ float scale, int phamin, int phamax, fitsfile *imagfits)
+{
+ int binx, biny, status=0, hdutype, anynull=0;
+ long nx, ny, nevents, j;
+ float *xpos, *ypos, *weight, *image;
+ double ncounts=0.;
+ char *pha;
+
+ FITS_movabs_hdu(infits, 2, &hdutype, &status);
+ nevents = cf_read_col(infits, TFLOAT, xcolumn, (void **) &xpos);
+ nevents = cf_read_col(infits, TFLOAT, ycolumn, (void **) &ypos);
+ nevents = cf_read_col(infits, TFLOAT, "WEIGHT", (void **) &weight);
+ nevents = cf_read_col(infits, TBYTE, "PHA", (void **) &pha);
+
+ FITS_movabs_hdu(imagfits, IMAGE_EXT, &hdutype, &status);
+ FITS_read_key(imagfits, TLONG, "NAXIS1", &nx, NULL, &status);
+ FITS_read_key(imagfits, TLONG, "NAXIS2", &ny, NULL, &status);
+ FITS_read_key(imagfits, TINT, "SPECBINX", &binx, NULL, &status);
+ FITS_read_key(imagfits, TINT, "SPECBINY", &biny, NULL, &status);
+
+ image = (float *) cf_malloc(sizeof(float)*nx*ny);
+ FITS_read_img(imagfits, TFLOAT, 1, nx*ny, NULL, image, &anynull, &status);
+
+ for (j=0; j<nevents; j++) {
+ if ((pha[j]>=phamin) && (pha[j]<=phamax)){
+ int x, y;
+ x = xpos[j];
+ y = ypos[j];
+ if ((x >= 0) && (x < nx*binx) && (y >= 0) && (y < ny*biny)) {
+ image[(y/biny)*nx + (x/binx)] -= weight[j]*scale;
+ ncounts += weight[j]*scale;
+ }
+ }
+ }
+ FITS_write_img(imagfits, TFLOAT, 1, nx*ny, image, &status);
+
+ free(image);
+ free(weight);
+ free(ypos);
+ free(xpos);
+ return ncounts;
+}
+
+
+int main(int argc, char *argv[])
+{
+ char *newfile = NULL, *detector = NULL, buffer[FLEN_CARD]={'\0'};
+ int subtract=0, force=0, newbinx=1, newbiny=1, status=0, optc;
+ float effmjd = 0.,scale=1.,old_scale;
+ int phamin=0, phamax=31;
+ int old_phamin, old_phamax;
+
+ char xcolumn[FLEN_CARD]={'\0'}, ycolumn[FLEN_CARD]={'\0'};
+
+ char opts[] = "hv:fsrcn:d:j:w:x:y:l:u:";
+ char usage[] =
+ "Usage:\n"
+ " cf_countmap [-hfsrc] [-v level] [-n filename] [-d segment]\n"
+ " [-w scale] [-j effmjd] [-x xbin] [-y ybin] \n"
+ " [-l phamin] [-u phamax] idffiles\n";
+ char option[] =
+ "Options:\n"
+ " -h: this help message\n"
+ " -v: verbosity level (=1; 0 is silent)\n"
+ " -f: force processing\n"
+ " -s: subtract files\n"
+ " -n: image map file\n"
+ " -d: image map segment (no default for new file)\n"
+ " -j: effective MJD (=0.0)\n"
+ " -w: weighting factor for image (=1.0)\n"
+ " -x: X-bin size (=1)\n"
+ " -y: Y-bin size (=1)\n"
+ " -r: use raw XY coordinates instead of farf\n"
+ " -c: use final XY coordinates instead of farf\n"
+ " -l: lower limit for pha (=0) \n"
+ " -u: upper limit for pha (=31) \n";
+
+ verbose_level = 1;
+
+ strcpy(xcolumn, "XFARF");
+ strcpy(ycolumn, "YFARF");
+
+ /* Check number of options and arguments */
+ while ((optc = getopt(argc, argv, opts)) != -1) {
+ switch(optc) {
+
+ case 'h':
+ printf("%s\n%s", usage, option);
+ return 0;
+ case 'v':
+ verbose_level = atoi(optarg);
+ break;
+ case 'f':
+ force = 1;
+ break;
+ case 's':
+ subtract = 1;
+ break;
+ case 'n':
+ newfile = optarg;
+ break;
+ case 'd':
+ detector = optarg;
+ break;
+ case 'j':
+ effmjd = atof(optarg);
+ break;
+ case 'w':
+ scale = atof(optarg);
+ break;
+ case 'x':
+ newbinx = atoi(optarg);
+ break;
+ case 'y':
+ newbiny = atoi(optarg);
+ break;
+ case 'r':
+ strcpy(xcolumn, "XRAW");
+ strcpy(ycolumn, "YRAW");
+ break;
+ case 'c':
+ strcpy(xcolumn, "X");
+ strcpy(ycolumn, "Y");
+ break;
+ case 'l':
+ phamin = atoi(optarg);
+ break;
+ case 'u':
+ phamax = atoi(optarg);
+ break;
+ }
+ }
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+ if (argc <= optind) {
+ cf_if_error("%s\nIncorrect number of program arguments", usage);
+ }
+
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin processing");
+
+ while (optind < argc) {
+
+ char imagfile[FLEN_CARD], rootname[FLEN_CARD]={'\0'};
+ char segment[FLEN_CARD] = {'\0'};
+ int hdutype;
+ FILE *fptr;
+ fitsfile *infits, *imagfits;
+
+ FITS_open_file(&infits, argv[optind], READONLY, &status);
+ FITS_movabs_hdu(infits, 1, &hdutype, &status);
+
+ if (newfile) {
+ strcpy(imagfile, newfile);
+
+ if ((fptr = fopen(imagfile, "r")))
+ fclose(fptr);
+ else {
+ if (!detector)
+ cf_if_error("No detector segment specified: "
+ "1A, 1B, 2A, 2B");
+ status = new_imagfile(imagfile, detector, effmjd,
+ newbinx, newbiny);
+ }
+ FITS_open_file(&imagfits, imagfile, READWRITE, &status);
+ FITS_movabs_hdu(imagfits, IMAGE_EXT, &hdutype, &status);
+ detector = buffer;
+ FITS_read_key(imagfits, TSTRING, "DETECTOR", detector, NULL,
+ &status);
+ }
+ else {
+ FITS_read_key(infits, TSTRING, "TCNT_MAP", imagfile, NULL,
+ &status);
+ detector = buffer;
+ FITS_read_key(infits, TSTRING, "DETECTOR", detector, NULL,
+ &status);
+
+ if ((fptr = fopen(cf_hist_file(imagfile), "r"))) {
+ fclose(fptr);
+ } else {
+ char keyword[FLEN_CARD]={'\0'}, segment[FLEN_CARD]={'\n'};
+ char linin[MAXCHARS]={'\0'}, filename[FLEN_CARD]={'\0'};
+ int interp;
+ float aftermjd;
+ FILE *master;
+ /*
+ * Get the effective MJD from the master calibration file
+ */
+ master = fopen(cf_parm_file(MASTER_CAL_FILE), "r");
+ if (master == NULL)
+ cf_if_error("Master calibration database file not found");
+ while (fgets(linin, MAXCHARS, master) != NULL) {
+ /*
+ * Check for comment lines
+ */
+ if ((linin[0] != '#') && (linin[0] != '\n')) {
+ sscanf(linin, "%4c%*2c%2c%*2c%13c%*2c%9f%*2c%1d",
+ keyword, segment, filename, &aftermjd, &interp);
+ if (strcmp(filename, imagfile) == 0) {
+ effmjd = aftermjd;
+ break;
+ }
+ }
+ }
+ fclose(master);
+
+ status = new_imagfile(cf_hist_file(imagfile), detector,
+ effmjd, newbinx, newbiny);
+ }
+ FITS_open_file(&imagfits, cf_hist_file(imagfile), READWRITE,
+ &status);
+ }
+
+ strncpy(rootname, argv[optind], ROOT_LENGTH);
+ FITS_read_key(infits, TSTRING, "DETECTOR", segment, NULL,
+ &status);
+ if (strcmp(detector, segment) != 0) {
+ cf_if_warning("%s segment is not the same as %s.",
+ argv[optind], imagfile);
+ }
+ else if (subtract) {
+ int index;
+ if ((index = _check_index(imagfits, rootname, &scale, &phamin, &phamax))) {
+
+ int hdutype;
+ float exptime=0., acctime=0.;
+ double nevents, accevnt=0.;
+
+ FITS_read_key(infits, TFLOAT, "EXPTIME", &exptime, NULL,
+ &status);
+
+
+ nevents = subtract_events(infits, xcolumn, ycolumn, scale, phamin, phamax, imagfits);
+
+ FITS_movabs_hdu(imagfits, IMAGE_EXT, &hdutype, &status);
+ FITS_read_key(imagfits, TDOUBLE, "NEVENTS", &accevnt,
+ NULL, &status);
+ FITS_read_key(imagfits, TFLOAT, "EXPTIME", &acctime,
+ NULL, &status);
+ accevnt -= nevents;
+ acctime -= exptime;
+ FITS_update_key(imagfits, TDOUBLE, "NEVENTS", &accevnt,
+ NULL, &status);
+ FITS_update_key(imagfits, TFLOAT, "EXPTIME", &acctime,
+ NULL, &status);
+
+ FITS_movabs_hdu(imagfits, INDEX_EXT, &hdutype, &status);
+ FITS_delete_rows(imagfits, index, 1, &status);
+ cf_verbose(1, "%s is subtracted from %s",
+ argv[optind], imagfile);
+ } else {
+ cf_if_warning("%s is not in %s.", argv[optind], imagfile);
+ }
+ } else {
+ if (!force && _check_index(imagfits, rootname, &old_scale, &old_phamin, &old_phamax)) {
+ cf_if_warning("%s is already in %s with scale = %f, phamin = %d, phamax=%d.",
+ argv[optind], imagfile, old_scale, old_phamin, old_phamax);
+ } else {
+ char obsdate[FLEN_CARD]={'\0'}, dateobs[FLEN_CARD]={'\0'};
+ char timeobs[FLEN_CARD]={'\0'};
+
+ float exptime=0., acctime=0.;
+ double nevents, accevnt=0;
+
+ FITS_movabs_hdu(imagfits, IMAGE_EXT, &hdutype, &status);
+ FITS_read_key(infits, TSTRING, "DATEOBS", dateobs, NULL,
+ &status);
+ FITS_read_key(infits, TSTRING, "TIMEOBS", timeobs, NULL,
+ &status);
+ obsdate[0] = '\0';
+ sprintf(obsdate, "%s%s%s", dateobs, "T", timeobs);
+ FITS_read_key(infits, TFLOAT, "EXPTIME", &exptime, NULL,
+ &status);
+
+ if ((force) && _check_index(imagfits, rootname, &old_scale, &old_phamin, &old_phamax))
+ cf_if_warning("FORCE MODE - %s is already in %s with scale = %f, phamin = %d, phamax=%d.",
+ argv[optind], imagfile, old_scale, old_phamin, old_phamax);
+
+ nevents = add_events(infits, xcolumn, ycolumn, scale, phamin, phamax, imagfits);
+
+ FITS_movabs_hdu(imagfits, IMAGE_EXT, &hdutype, &status);
+ FITS_read_key(imagfits, TDOUBLE, "NEVENTS", &accevnt,
+ NULL, &status);
+ FITS_read_key(imagfits, TFLOAT, "EXPTIME", &acctime,
+ NULL, &status);
+ accevnt += nevents;
+ acctime += exptime;
+ FITS_update_key(imagfits, TDOUBLE, "NEVENTS", &accevnt,
+ NULL, &status);
+ FITS_update_key(imagfits, TFLOAT, "EXPTIME", &acctime,
+ NULL, &status);
+
+ status = _update_index(imagfits, rootname, obsdate, nevents,
+ exptime, scale, phamin, phamax);
+ cf_verbose(1, "%s is added to %s", argv[optind], imagfile);
+ }
+ }
+ FITS_close_file(imagfits, &status);
+ FITS_close_file(infits, &status);
+ optind ++;
+ }
+
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Finished processing");
+ return 0;
+}
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;
+}
diff --git a/src/fuv/cf_flux_calibrate.c b/src/fuv/cf_flux_calibrate.c
new file mode 100644
index 0000000..2ad92e6
--- /dev/null
+++ b/src/fuv/cf_flux_calibrate.c
@@ -0,0 +1,173 @@
+/*****************************************************************************
+ * Johns Hopkins University
+ * Center For Astrophysical Sciences
+ * FUSE
+ *****************************************************************************
+ *
+ * Synopsis: cf_flux_calibrate options intermediate_file
+ *
+ * Description: Flux calibrate each photon assigned to a channel.
+ *
+ * ** NEXT STEPS COMMENTED OUT FOR NOW **
+ * For each photon in the target aperture,
+ * apply astigmatism correction to XFARF array and use the
+ * resulting coordinate system to apply flat-field and
+ * worm corrections to photon weights.
+ *
+ * By default all corrections are performed. Command line
+ * options allow one or more corrections to be omitted.
+ *
+ * Arguments: input_file Intermediate data file
+ *
+ * Calibration files: AEFF_CAL
+ *
+ * Returns: 0 on successful completion
+ *
+ * Calls: cf_effective_area, cf_astig_farf, cf_flat_field,
+ * cf_worm_correction
+ *
+ * History: 12/06/02 wvd 1.1 Begin work
+ * 02/12/03 wvd 1.2 Read ERGCM2S from IDF
+ * Move cf_convert_to_ergs to end of
+ * module
+ * 03/11/03 wvd 1.3 Change channel to type char
+ * 03/12/03 wvd 1.4 If ERGCM2S[i] is undefined,
+ * set channel[i] = 0
+ * 03/19/03 peb 1.5 Added verbose_level option.
+ * 05/06/03 rdr 1.6 Read from column ERGCM2 rather
+ * than ERGCM2S
+ * 05/07/03 wvd 1.7 Change ergcm2s to ergcm2
+ * 06/11/03 wvd 1.8 Pass datatype to cf_read_col and
+ * cf_write_col.
+ * Comment out call to cf_astig_farf().
+ * 08/25/03 wvd 1.9 Change coltype from string to int in
+ * cf_read_col and cf_write_col.
+ * 09/17/03 wvd 1.10 Add some documentation.
+ * 09/17/03 wvd 1.11 Change channel to unsigned char.
+ * 05/20/05 wvd 1.12 Clean up i/o.
+ *
+ ****************************************************************************/
+
+#include <unistd.h>
+#include <stdlib.h>
+#include "calfuse.h"
+
+static char CF_PRGM_ID[] = "cf_flux_calibrate";
+static char CF_VER_NUM[] = "1.12";
+
+int
+main(int argc, char *argv[])
+{
+ unsigned char *channel=NULL;
+ int effective_area=1, astig_farf=1, flat_field=1, worm_correction=1;
+ int status=0, hdutype, optc;
+ long nevents, nseconds;
+ float *time=NULL, *ttime=NULL, *xfarf=NULL, *yfarf=NULL, *weight=NULL;
+ float *lambda=NULL, *ycentl=NULL, *ycents=NULL, *ergcm2=NULL;
+ fitsfile *header;
+
+ char opts[] = "heafwv:";
+ char usage[] =
+ "Usage:\n"
+ " cf_flux_calibrate [-heafw] [-v level] idf_file\n";
+ char option[] =
+ "Options:\n"
+ " -h: this help message\n"
+ " -v: verbosity level (default is 1; 0 is silent)\n"
+ " -e: no effective-area calibration\n"
+ " -a: no astigmatism correction of XFARF array\n"
+ " -f: no flat-field correction\n"
+ " -w: no worm correction\n";
+
+ verbose_level = 1;
+
+ while ((optc = getopt(argc, argv, opts)) != -1) {
+ switch(optc) {
+
+ case 'h':
+ printf("%s\n%s", usage, option);
+ return 0;
+ case 'v':
+ verbose_level = atoi(optarg);
+ break;
+ case 'e':
+ effective_area = 0;
+ break;
+ case 'a':
+ astig_farf = 0;
+ break;
+ case 'f':
+ flat_field = 0;
+ break;
+ case 'w':
+ worm_correction = 0;
+ break;
+ }
+ }
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+ if (argc != optind+1) {
+ printf("%s", usage);
+ cf_if_error("Incorrect number of command-line arguments");
+ }
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");
+
+ FITS_open_file(&header, argv[optind], READWRITE, &status);
+
+ FITS_movabs_hdu(header, 2, &hdutype, &status);
+ nevents = cf_read_col(header, TFLOAT, "TIME", (void **) &time);
+ nevents = cf_read_col(header, TFLOAT, "WEIGHT", (void **) &weight);
+ nevents = cf_read_col(header, TFLOAT, "XFARF", (void **) &xfarf);
+ nevents = cf_read_col(header, TFLOAT, "YFARF", (void **) &yfarf);
+ nevents = cf_read_col(header, TBYTE, "CHANNEL", (void **) &channel);
+ nevents = cf_read_col(header, TFLOAT, "LAMBDA", (void **) &lambda);
+ nevents = cf_read_col(header, TFLOAT, "ERGCM2", (void **) &ergcm2);
+
+ FITS_movabs_hdu(header, 4, &hdutype, &status) ;
+ nseconds = cf_read_col(header, TFLOAT, "TIME", (void **) &ttime) ;
+ nseconds = cf_read_col(header, TFLOAT, "YCENT_LIF", (void **) &ycentl) ;
+ nseconds = cf_read_col(header, TFLOAT, "YCENT_SIC", (void **) &ycents) ;
+
+ FITS_movabs_hdu(header, 1, &hdutype, &status);
+
+ /***************************************************************************
+ * Since we currently have no flat-field or worm correction,
+ * these subroutines are commented out.
+ *
+ if (astig_farf)
+ cf_astig_farf(header, nevents, xfarf, yfarf, channel, time,
+ nseconds, ttime, ycentl, ycents);
+ if (flat_field)
+ cf_flat_field(header, nevents, time, weight, xfarf, yfarf, channel,
+ nseconds, ttime, ycentl, ycents);
+
+ if (worm_correction)
+ cf_worm_correction(header, nevents, time, weight, xfarf, yfarf, channel,
+ nseconds, ttime, ycentl, ycents);
+ ***************************************************************************/
+
+ if (effective_area)
+ cf_convert_to_ergs(header, nevents, weight, ergcm2, channel, lambda);
+
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+
+ FITS_movabs_hdu(header, 2, &hdutype, &status);
+ cf_write_col(header, TFLOAT, "WEIGHT", (void *) weight, nevents);
+ cf_write_col(header, TBYTE, "CHANNEL", (void *) channel, nevents);
+ cf_write_col(header, TFLOAT, "ERGCM2", (void *) ergcm2, nevents);
+
+ FITS_close_file(header, &status);
+
+ free(time);
+ free(weight);
+ free(xfarf);
+ free(yfarf);
+ free(channel);
+ free(lambda);
+ free(ergcm2);
+ free(ttime);
+ free(ycentl);
+ free(ycents);
+
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done processing");
+ return 0;
+}
diff --git a/src/fuv/cf_gainmap.c b/src/fuv/cf_gainmap.c
new file mode 100644
index 0000000..63ed3f7
--- /dev/null
+++ b/src/fuv/cf_gainmap.c
@@ -0,0 +1,610 @@
+/*****************************************************************************
+ * Johns Hopkins University
+ * Center For Astrophysical Sciences
+ * FUSE
+ *****************************************************************************
+ *
+ * Synopsis: cf_gainmap input_files
+ *
+ * Description: Adds the counts from a list of input ttag idf files into
+ * the gain map. The accumulated gain map cube is
+ * determined by reading the GAIN_CAL keyword in the
+ * input file's header. To prevent multiple copies
+ * of reprocessed input files from being added to
+ * the total exposure map, the gain map contains an
+ * index list of the files which have been added.
+ * The new file is skipped if it is in the list.
+ *
+ *
+ * Arguments: input_files Input ttag IDF FITS file name
+ *
+ * Calibration files required: None
+ *
+ * Returns: int 0 = successful completion
+ *
+ * Note: This routine requires significantly more memory than your
+ * average routine. Its design has not been optimized to
+ * reduce memory usage.
+ *
+ * History: 03/31/03 1.1 peb Finished work.
+ * Based on cf_ttag_countmap.
+ * 06/11/03 1.2 wvd Pass datatype to cf_read_col
+ * 08/25/03 1.3 wvd Change coltype from string to int in
+ * cf_read_col.
+ * 01/09/04 1.4 bjg Change CHIDCOR1 to CHID_COR
+ * 02/17/04 1.5 bjg Added scaling option
+ * histogram files given in the argument
+ * list are ignored
+ * 02/26/04 1.6 bjg Bug Fix. Routine name changed from
+ * cf_ttag_gainmap to cf_gainmap.
+ * Added option to select coordinates
+ * to use: RAW, FARF (default), final
+ * Print ignored HIST files
+ * 04/08/04 1.7 bjg Remove unused function get_xy_columns
+ * Use cf_verbose instead of printf
+ * 04/27/04 1.8 bjg Added lower and upper limits option
+ * for PHA
+ * Change author name of gainfile
+ * 08/26/04 1.9 bjg In force mode, issue a warning
+ * if file already included.
+ * 03/02/05 1.10 wvd Use FLOAT_IMG to insure that output
+ * images are written as floats.
+ * 03/02/05 1.11 wvd Don't write AUTHOR keyword to header.
+ * 04/11/05 1.12 wvd Force detector string to be upper case.
+ * Add -t option, which prevents creation
+ * of detector image for PHA=31. It's
+ * required if you don't bin the data.
+ * 06/03/05 1.13 wvd Include ctype.h
+ * 09/09/05 1.14 wvd Update list of allowed options.
+ * 08/24/07 1.15 bot Changed nx and ny to long and called
+ * as TLONG ; changed binx and biny to int
+ * and called as TINT.
+ *
+ ****************************************************************************/
+
+#include <ctype.h>
+#include <string.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include "calfuse.h"
+
+#define ROOT_LENGTH 13
+#define MAXCHARS 120
+#define MASTER_CAL_FILE "master_calib_file.dat"
+#define IMAGE_EXT 1
+#define INDEX_EXT 34
+
+static char CF_PRGM_ID[] = "cf_gainmap";
+static char CF_VER_NUM[] = "1.15";
+
+
+static int
+new_gainfile(char *newfile, char *detector, float effmjd, int binx, int biny, int truncate)
+{
+ char date[FLEN_CARD]={'\0'};
+ char *ttype[] = {"ROOTNAME", "OBSDATE", "COUNTS", "EXPTIME", "ADDDATE", "SCALE", "PHAMIN", "PHAMAX"};
+ char *tform[] = {"13A", "19A", "1D", "1E", "19A", "1E", "1I", "1I"};
+ char *tunit[] = {" ", " ", "counts", "seconds", " ", " ", " ", " "};
+ int naxis=2, status=0, version=1, tref=0, tfields=8, phaext;
+ long axis[2], naxis2=0;
+ float exptime=0.;
+ double ncounts=0.;
+ fitsfile *gainfits;
+
+ cf_verbose(3, "Entering subroutine new_gainfile");
+ cf_verbose(3, "Creating file %s", newfile);
+
+ FITS_create_file(&gainfits, newfile, &status);
+
+ FITS_create_img(gainfits, 8, 0, NULL, &status);
+ FITS_write_key(gainfits, TSTRING, "CALFTYPE", "GAIN",
+ "Calibration file type", &status);
+ FITS_write_key(gainfits, TINT, "CALFVERS", &version,
+ "Calibration file version", &status);
+ FITS_write_key(gainfits, TSTRING, "DETECTOR", detector,
+ "detector (1A, 1B, 2A, 2B)", &status);
+ FITS_write_key(gainfits, TFLOAT, "EFFMJD", &effmjd,
+ "Date on which file should be applied (MJD)", &status);
+ fits_get_system_time(date, &tref, &status);
+ FITS_write_key(gainfits, TSTRING, "DATE", date,
+ "file creation date (YYYY-MM-DDThh:mm:ss UTC)", &status);
+ FITS_write_key(gainfits, TDOUBLE, "NEVENTS", &ncounts,
+ "Accumulated weighted events", &status);
+ FITS_write_key(gainfits, TFLOAT, "EXPTIME", &exptime,
+ "Accumulated exposure time", &status);
+
+ axis[0] = NXMAX/binx;
+ axis[1] = NYMAX/biny;
+ for (phaext=0; phaext<32-truncate; phaext++) {
+
+ cf_verbose(3, "Creating extension %d", phaext);
+ FITS_create_img(gainfits, FLOAT_IMG, naxis, axis, &status);
+ FITS_write_key(gainfits, TINT, "PHA", &phaext,
+ "PHA value of image", &status);
+ FITS_write_key(gainfits, TINT, "SPECBINX", &binx,
+ "Binsize in detector X coordinate", &status);
+ FITS_write_key(gainfits, TINT, "SPECBINY", &biny,
+ "Binsize in detector Y coordinate", &status);
+
+ }
+
+ cf_verbose(3, "Creating binary table extension");
+ FITS_create_tbl(gainfits, BINARY_TBL, naxis2, tfields, ttype, tform, tunit,
+ "PROCESSED FILES", &status);
+
+ cf_verbose(3, "Closing file %s", newfile);
+ FITS_close_file(gainfits, &status);
+ cf_verbose(3, "Exiting new_gainfile");
+ return status;
+}
+
+static int
+_check_index(fitsfile *gainfits, char *rootname, float *scale, int *phamin,
+ int *phamax, int truncate)
+{
+ int status=0, anynull=0, hdutype;
+ float *scal;
+ int *phami, *phama;
+ long j, nrows, frow=1, felem=1;
+ char **filename, strnull[] = " ";
+
+ FITS_movabs_hdu(gainfits, INDEX_EXT-truncate, &hdutype, &status);
+ FITS_read_key(gainfits, TLONG, "NAXIS2", &nrows, NULL, &status);
+
+ filename = (char **) cf_malloc(nrows * sizeof(char *));
+ filename[0] = (char *) cf_malloc(nrows*(ROOT_LENGTH+1)*sizeof(char));
+ scal = (float *) cf_malloc(nrows * sizeof(float));
+ phami = (int *) cf_malloc(nrows * sizeof(int));
+ phama = (int *) cf_malloc(nrows * sizeof(int));
+ for (j=1; j<nrows; j++)
+ filename[j] = filename[j-1]+(ROOT_LENGTH+1);
+
+ FITS_read_col(gainfits, TSTRING, 1, frow, felem, nrows, strnull,
+ filename, &anynull, &status);
+ FITS_read_col(gainfits, TFLOAT, 6, frow, felem, nrows, strnull,
+ scal, &anynull, &status);
+ FITS_read_col(gainfits, TINT, 7, frow, felem, nrows, strnull,
+ phami, &anynull, &status);
+ FITS_read_col(gainfits, TINT, 8, frow, felem, nrows, strnull,
+ phama, &anynull, &status);
+
+ for (j=0; j<nrows; j++) {
+ if (strncmp(filename[j], rootname, strlen(rootname)) == 0) {
+ *scale = scal[j];
+ *phamin = phami[j];
+ *phamax = phama[j];
+ free(filename[0]);
+ free(filename);
+ free(scal);
+ free(phami);
+ free(phama);
+ return j+1;
+ }
+ }
+ free(filename[0]);
+ free(filename);
+ free(scal);
+ free(phami);
+ free(phama);
+ return 0;
+}
+
+static int
+_update_index(fitsfile *gainfits, char *rootname, char *obsdate,
+ double ncounts, float exptime, float scale, int phamin, int phamax,
+ int truncate)
+{
+ char *strptr[1], adddate[FLEN_CARD]={'\0'};
+ int status=0, hdutype, tref=0;
+ long nrows;
+
+ FITS_movabs_hdu(gainfits, INDEX_EXT-truncate, &hdutype, &status);
+ FITS_read_key(gainfits, TLONG, "NAXIS2", &nrows, NULL, &status);
+
+ strptr[0] = rootname;
+ FITS_write_col(gainfits, TSTRING, 1, nrows+1, 1, 1, strptr, &status);
+ strptr[0] = obsdate;
+ FITS_write_col(gainfits, TSTRING, 2, nrows+1, 1, 1, strptr, &status);
+ FITS_write_col(gainfits, TDOUBLE, 3, nrows+1, 1, 1, &ncounts, &status);
+ FITS_write_col(gainfits, TFLOAT, 4, nrows+1, 1, 1, &exptime, &status);
+ fits_get_system_time(adddate, &tref, &status);
+ strptr[0] = adddate;
+ fits_write_col(gainfits, TSTRING, 5, nrows+1, 1, 1, strptr, &status);
+ FITS_write_col(gainfits, TFLOAT, 6, nrows+1, 1, 1, &scale, &status);
+ fits_write_col(gainfits, TINT, 7, nrows+1, 1, 1, &phamin, &status);
+ FITS_write_col(gainfits, TINT, 8, nrows+1, 1, 1, &phamax, &status);
+ return 0;
+}
+
+
+static double
+add_events(fitsfile *infits, char *xcolumn, char *ycolumn, float scale,
+ int phamin, int phamax, fitsfile *gainfits)
+{
+ char *pha;
+ int binx, biny, phaext, status=0, hdutype, anynull=0;
+ long nevents, j;
+ long nx, ny;
+ float *xpos, *ypos, *weight, *image;
+ double ncounts=0.;
+
+ cf_verbose(3,"Entering add_events");
+
+ FITS_movabs_hdu(infits, 2, &hdutype, &status);
+ nevents = cf_read_col(infits, TFLOAT, xcolumn, (void **) &xpos);
+ nevents = cf_read_col(infits, TFLOAT, ycolumn, (void **) &ypos);
+ nevents = cf_read_col(infits, TBYTE, "PHA", (void **) &pha);
+ nevents = cf_read_col(infits, TFLOAT, "WEIGHT", (void **) &weight);
+
+ for (phaext=phamin; phaext<=phamax; phaext++) {
+ FITS_movabs_hdu(gainfits, phaext+2, &hdutype, &status);
+ if (hdutype == BINARY_TBL) {
+ cf_if_warning("Can't find image extension for PHA = %d", phaext);
+ break;
+ }
+ FITS_read_key(gainfits, TLONG, "NAXIS1", &nx, NULL, &status);
+ FITS_read_key(gainfits, TLONG, "NAXIS2", &ny, NULL, &status);
+ FITS_read_key(gainfits, TINT, "SPECBINX", &binx, NULL, &status);
+ FITS_read_key(gainfits, TINT, "SPECBINY", &biny, NULL, &status);
+
+ image = (float *) cf_malloc(sizeof(float)*nx*ny);
+ FITS_read_img(gainfits, TFLOAT, 1, nx*ny, NULL, image,
+ &anynull, &status);
+
+ for (j=0; j<nevents; j++) {
+ if (pha[j] == phaext) {
+ int x, y;
+ x = xpos[j];
+ y = ypos[j];
+ if ((x >= 0) && (x < nx*binx) && (y >= 0) && (y < ny*biny)) {
+ image[(y/biny)*nx + (x/binx)] += weight[j]*scale;
+ ncounts += weight[j]*scale;
+ }
+ }
+ }
+ FITS_write_img(gainfits, TFLOAT, 1, nx*ny, image, &status);
+ free(image);
+ }
+
+ free(weight);
+ free(pha);
+ free(ypos);
+ free(xpos);
+ return ncounts;
+}
+
+static double
+subtract_events(fitsfile *infits, char *xcolumn, char *ycolumn,
+ float scale, int phamin, int phamax, fitsfile *gainfits)
+{
+ char *pha;
+ long nx, ny;
+ int phaext, status=0, hdutype, anynull=0, binx, biny;
+ long nevents, j;
+ float *xpos, *ypos, *weight, *image;
+ double ncounts=0.;
+
+ FITS_movabs_hdu(infits, 2, &hdutype, &status);
+ nevents = cf_read_col(infits, TFLOAT, xcolumn, (void **) &xpos);
+ nevents = cf_read_col(infits, TFLOAT, ycolumn, (void **) &ypos);
+ nevents = cf_read_col(infits, TBYTE, "PHA", (void **) &pha);
+ nevents = cf_read_col(infits, TFLOAT, "WEIGHT", (void **) &weight);
+
+ for (phaext=phamin; phaext<=phamax; phaext++) {
+ FITS_movabs_hdu(gainfits, phaext+2, &hdutype, &status);
+ FITS_read_key(gainfits, TLONG, "NAXIS1", &nx, NULL, &status);
+ FITS_read_key(gainfits, TLONG, "NAXIS2", &ny, NULL, &status);
+ FITS_read_key(gainfits, TINT, "SPECBINX", &binx, NULL, &status);
+ FITS_read_key(gainfits, TINT, "SPECBINY", &biny, NULL, &status);
+
+ image = (float *) cf_malloc(sizeof(float) * nx*ny);
+ FITS_read_img(gainfits, TFLOAT, 1, nx*ny, NULL, image,
+ &anynull, &status);
+
+ for (j=0; j<nevents; j++) {
+ if (pha[j] == phaext) {
+ int x, y;
+ x = xpos[j];
+ y = ypos[j];
+ if ((x >=0) && (x < nx*binx) && (y >= 0) && (y < ny*biny)) {
+ image[(y/biny)*nx + (x/binx)] -= weight[j]*scale;
+ ncounts += weight[j]*scale;
+ }
+ }
+ }
+ FITS_write_img(gainfits, TFLOAT, 1, nx*ny, image, &status);
+ free(image);
+ }
+
+ free(weight);
+ free(pha);
+ free(ypos);
+ free(xpos);
+ return ncounts;
+}
+
+
+int main(int argc, char *argv[])
+{
+ char *newfile = NULL, *detector = NULL, buffer[FLEN_CARD]={'\0'}, instmode[FLEN_CARD]={'\0'};
+ int subtract=0, force=0, newbinx=8, newbiny=4, status=0, optc;
+ int truncate=0;
+ float effmjd = 0.,scale=1., old_scale;
+ int phamin=0, phamax=31;
+ int old_phamin, old_phamax;
+
+ char xcolumn[FLEN_CARD]={'\0'}, ycolumn[FLEN_CARD]={'\0'};
+
+ char opts[] = "hv:fsrtcn:d:j:w:x:y:l:u:";
+ char usage[] =
+ "Usage:\n"
+ " cf_gainmap [-hfsrtc] [-v level] [-n filename] [-d segment]\n"
+ " [-w scale] [-j effmjd] [-x xbin] [-y ybin] \n"
+ " [-l phamin] [-u phamax] idffiles\n";
+ char option[] =
+ "Options:\n"
+ " -h: this help message\n"
+ " -v: verbosity level (=1; 0 is silent)\n"
+ " -f: force the files\n"
+ " -s: subtract files\n"
+ " -n: gain map file\n"
+ " -d: gain map segment (no default for new file)\n"
+ " -j: effective MJD (=0.0)\n"
+ " -w: weighting factor for image (=1.0)\n"
+ " -x: X-bin size (=8)\n"
+ " -y: Y-bin size (=4)\n"
+ " -r: use raw XY coordinates instead of farf\n"
+ " -t: truncate map at PHA = 30\n"
+ " -c: use final XY coordinates instead of farf\n"
+ " -l: lower limit for pha (=0) \n"
+ " -u: upper limit for pha (=31) \n";
+
+ verbose_level = 1;
+
+ strcpy(xcolumn, "XFARF");
+ strcpy(ycolumn, "YFARF");
+
+
+
+ /* Check number of options and arguments */
+ while ((optc = getopt(argc, argv, opts)) != -1) {
+ switch(optc) {
+
+ case 'h':
+ printf("%s\n%s", usage, option);
+ return 0;
+ case 'v':
+ verbose_level = atoi(optarg);
+ break;
+ case 'f':
+ force = 1;
+ break;
+ case 's':
+ subtract = 1;
+ break;
+ case 'n':
+ newfile = optarg;
+ break;
+ case 'd':
+ detector = optarg;
+ break;
+ case 'j':
+ effmjd = atof(optarg);
+ break;
+ case 'w':
+ scale = atof(optarg);
+ break;
+ case 'x':
+ newbinx = atoi(optarg);
+ break;
+ case 'y':
+ newbiny = atoi(optarg);
+ break;
+ case 'r':
+ strcpy(xcolumn, "XRAW");
+ strcpy(ycolumn, "YRAW");
+ break;
+ case 't':
+ truncate = 1;
+ break;
+ case 'c':
+ strcpy(xcolumn, "X");
+ strcpy(ycolumn, "Y");
+ break;
+ case 'l':
+ phamin = atoi(optarg);
+ break;
+ case 'u':
+ phamax = atoi(optarg);
+ break;
+ }
+ }
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+ if (argc <= optind)
+ cf_if_error("%s\nIncorrect number of program arguments", usage);
+
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin processing");
+
+ while (optind < argc) {
+
+ char gainfile[FLEN_CARD], rootname[FLEN_CARD]={'\0'};
+ char segment[FLEN_CARD] = {'\0'};
+ int hdutype;
+ FILE *fptr;
+ fitsfile *infits, *gainfits;
+
+ FITS_open_file(&infits, argv[optind], READONLY, &status);
+ FITS_movabs_hdu(infits, 1, &hdutype, &status);
+
+ FITS_read_key(infits, TSTRING, "INSTMODE", instmode, NULL,
+ &status);
+
+ if (strncmp(instmode,"TTAG",4)){
+ FITS_close_file(infits, &status);
+ cf_verbose(1,"Skipped HIST file: %s\n",argv[optind]);
+ optind ++;
+ continue;
+ }
+
+ if (newfile) {
+ cf_verbose(3, "newfile is TRUE");
+ strcpy(gainfile, newfile);
+
+ if ((fptr = fopen(gainfile, "r")))
+ fclose(fptr);
+ else {
+ if (!detector)
+ cf_if_error("No detector segment specified: "
+ "1A, 1B, 2A, 2B");
+ detector[1] = toupper(detector[1]);
+ status = new_gainfile(gainfile, detector, effmjd,
+ newbinx, newbiny, truncate);
+ }
+ FITS_open_file(&gainfits, gainfile, READWRITE, &status);
+ FITS_movabs_hdu(gainfits, 1, &hdutype, &status);
+ detector = buffer;
+ FITS_read_key(gainfits, TSTRING, "DETECTOR", detector, NULL,
+ &status);
+ } else {
+ cf_verbose(3, "newfile is FALSE");
+ FITS_read_key(infits, TSTRING, "GAIN_MAP", gainfile, NULL,
+ &status);
+ detector = buffer;
+ FITS_read_key(infits, TSTRING, "DETECTOR", detector, NULL,
+ &status);
+
+ if ((fptr = fopen(cf_hist_file(gainfile), "r"))) {
+ fclose(fptr);
+ } else {
+ char keyword[FLEN_CARD]={'\0'}, segment[FLEN_CARD]={'\n'};
+ char linin[MAXCHARS]={'\0'}, filename[FLEN_CARD]={'\0'};
+ int interp;
+ float aftermjd;
+ FILE *master;
+ /*
+ * Get the effective MJD from the master calibration file
+ */
+ master = fopen(cf_parm_file(MASTER_CAL_FILE), "r");
+ if (master == NULL)
+ cf_if_error("Master calibration database file not found");
+ while (fgets(linin, MAXCHARS, master) != NULL) {
+ /*
+ * Check for comment lines
+ */
+ if ((linin[0] != '#') && (linin[0] != '\n')) {
+ sscanf(linin, "%4c%*2c%2c%*2c%13c%*2c%9f%*2c%1d",
+ keyword, segment, filename, &aftermjd, &interp);
+ if (strcmp(filename, gainfile) == 0) {
+ effmjd = aftermjd;
+ break;
+ }
+ }
+ }
+ fclose(master);
+
+ status = new_gainfile(cf_hist_file(gainfile), detector,
+ effmjd, newbinx, newbiny, truncate);
+ }
+ FITS_open_file(&gainfits, cf_hist_file(gainfile), READWRITE,
+ &status);
+ }
+
+ strncpy(rootname, argv[optind], ROOT_LENGTH);
+ FITS_read_key(infits, TSTRING, "DETECTOR", segment, NULL, &status);
+ if (strcmp(detector, segment) != 0) {
+ cf_if_warning("%s segment is not the same as %s",
+ argv[optind], gainfile);
+ }
+ else if (subtract) {
+ int index;
+ if ((index = _check_index(gainfits, rootname, &scale, &phamin,
+ &phamax, truncate))) {
+
+ int hdutype;
+ float exptime=0., acctime=0.;
+ double ncounts=0, accevnt=0;
+
+ FITS_read_key(infits, TFLOAT, "EXPTIME", &exptime, NULL,
+ &status);
+
+
+ ncounts = subtract_events(infits, xcolumn, ycolumn, scale, phamin,
+ phamax-truncate, gainfits);
+
+ FITS_movabs_hdu(gainfits, IMAGE_EXT, &hdutype, &status);
+ FITS_read_key(gainfits, TDOUBLE, "NEVENTS", &accevnt,
+ NULL, &status);
+ FITS_read_key(gainfits, TFLOAT, "EXPTIME", &acctime,
+ NULL, &status);
+ accevnt -= ncounts;
+ acctime -= exptime;
+ FITS_update_key(gainfits, TDOUBLE, "NEVENTS", &accevnt,
+ NULL, &status);
+ FITS_update_key(gainfits, TFLOAT, "EXPTIME", &acctime,
+ NULL, &status);
+
+ FITS_movabs_hdu(gainfits, INDEX_EXT-truncate, &hdutype, &status);
+ FITS_delete_rows(gainfits, index, 1, &status);
+ cf_verbose(1, "%s is subtracted from %s",
+ argv[optind], gainfile);
+ } else {
+ cf_if_warning("%s is not in %s.", argv[optind], gainfile);
+ }
+ } else {
+ cf_verbose(3, "Adding new data file.");
+ if (!force && _check_index(gainfits, rootname, &old_scale,
+ &old_phamin, &old_phamax, truncate)) {
+ cf_if_warning("%s is already in %s with scale = %f, "
+ "phamin = %d, phamax=%d.",
+ argv[optind], gainfile, old_scale, old_phamin, old_phamax);
+ } else {
+ char obsdate[FLEN_CARD]={'\0'}, dateobs[FLEN_CARD]={'\0'};
+ char timeobs[FLEN_CARD]={'\0'};
+
+ float exptime=0., acctime=0.;
+ double ncounts=0, accevnt=0;
+
+ FITS_read_key(infits, TSTRING, "DATEOBS", dateobs, NULL,
+ &status);
+ FITS_read_key(infits, TSTRING, "TIMEOBS", timeobs, NULL,
+ &status);
+ obsdate[0] = '\0';
+ sprintf(obsdate, "%s%s%s", dateobs, "T", timeobs);
+ FITS_read_key(infits, TFLOAT, "EXPTIME", &exptime, NULL,
+ &status);
+
+ if ((force) && _check_index(gainfits, rootname, &old_scale,
+ &old_phamin, &old_phamax, truncate))
+ cf_if_warning("FORCE MODE - %s is already in %s with scale = %f, "
+ "phamin = %d, phamax=%d.",
+ argv[optind], gainfile, old_scale, old_phamin, old_phamax);
+
+ ncounts = add_events(infits, xcolumn, ycolumn, scale, phamin,
+ phamax-truncate, gainfits);
+
+ FITS_movabs_hdu(gainfits, IMAGE_EXT, &hdutype, &status);
+ FITS_read_key(gainfits, TDOUBLE, "NEVENTS", &accevnt,
+ NULL, &status);
+ FITS_read_key(gainfits, TFLOAT, "EXPTIME", &acctime,
+ NULL, &status);
+ accevnt += ncounts;
+ acctime += exptime;
+ FITS_update_key(gainfits, TDOUBLE, "NEVENTS", &accevnt,
+ NULL, &status);
+ FITS_update_key(gainfits, TFLOAT, "EXPTIME", &acctime,
+ NULL, &status);
+
+ status = _update_index(gainfits, rootname, obsdate, ncounts,
+ exptime, scale, phamin, phamax, truncate);
+ cf_verbose(1, "%s is added to %s", argv[optind], gainfile);
+ }
+ }
+ FITS_close_file(gainfits, &status);
+ FITS_close_file(infits, &status);
+ optind ++;
+ }
+
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Finished processing");
+ return 0;
+}
diff --git a/src/fuv/cf_hist_init.c b/src/fuv/cf_hist_init.c
new file mode 100644
index 0000000..5829707
--- /dev/null
+++ b/src/fuv/cf_hist_init.c
@@ -0,0 +1,653 @@
+/*****************************************************************************
+ * Johns Hopkins University
+ * Center For Astrophysical Sciences
+ * FUSE
+ *****************************************************************************
+ *
+ * Synopsis: cf_hist_init input_file output_file
+ *
+ * Description: Reads a raw hist file and generates a list of the pixels
+ * that contain counts. This list has the same format
+ * as a TTAG list, so analysis can be carried out using
+ * many of the modules in the TTAG pipeline.
+ *
+ * The output file will have the same format as a TTAG
+ * file, with the following definitions:
+ *
+ * TIME FLOAT all equal to the midpoint of the exp.
+ * XRAW INT raw x position
+ * YRAW INT raw y position (midpoint of 8 pixel bin)
+ * PHA BYTE all 20 - indicating a good point
+ * WEIGHT FLOAT number of photons falling in that bin
+ * XFARF FLOAT x position in geometrically
+ * corrected frame
+ * YFARF FLOAT y position in geometrically
+ * corrected frame
+ * X FLOAT final x position (after correcting
+ * for motions)
+ * Y FLOAT final y position (after correcting
+ * for motions)
+ * CHANNEL CHAR aperture ID for the photon
+ * TIMEFLGS CHAR temporal flags - various bits set
+ * LOC_FLGS CHAR location flags - various bits set
+ * LAMBDA FLOAT wavelength for each photon
+ * ERGCM2 FLOAT flux calibration for each photon
+ *
+ * This routine calls cf_fuv_init to populate the calibration
+ * entries in the header.
+ *
+ * In addition, the routine creates a timeline table containing
+ * the orbital information, count rates, and high-voltage
+ * information to be used by later pipeline modules.
+ *
+ *
+ * Returns: Zero upon successful completion.
+ *
+ * History: 05/09/03 v1.1 rdr Adapted from cf_ttag_init
+ * 06/04/03 v1.4 wvd Add call to cf_proc_check
+ * 06/10/03 v1.5 rdr Scan for hot pixels and remove
+ * 06/11/03 v1.6 wvd Compress data with TSCALE and TZERO
+ * 06/11/03 v1.7 wvd Changed calfusettag.h to calfuse.h
+ * 07/17/03 v1.8 wvd Move various subroutines to libcf.
+ * 07/22/03 v1.9 wvd Read all extensions of raw hist file.
+ * Test for hot pixels only in active
+ * region and only if img[j]>15.
+ * 07/28/03 v1.10 wvd Close ELEC_CAL.
+ * 10/26/03 v1.11 wvd Increase morekeys to 38.
+ * 11/05/03 v1.12 wvd Added stdlib.h and unistd.h.
+ * 12/05/03 bjg Modified cf_make_pixel_list
+ * to unbin data in y.
+ * 12/10/03 v1.13 bjg Unbin is now optional
+ * Update keyword SPECBINY
+ * 03/03/04 v1.14 rdr Exit if no photons found
+ * 03/05/04 v1.15 rdr Warning if data is missing
+ * 03/09/04 v1.16 rdr Fix error in correcting hot pixels
+ * 03/26/04 v1.17 rdr Continue when there is no data
+ * 03/30/04 v1.20 wvd Clean up i/o.
+ * 03/31/04 v1.21 bjg functions return EXIT_SUCCESS instead
+ * of zero.
+ * >>>Include math.h for fmod<<<
+ * Type conversion for
+ * double fmod(double, double)
+ * 04/06/04 v1.22 bjg Remove unused variables
+ * Change format to match arg type
+ * in printf
+ * Remove static keyword in struct key
+ * definition
+ * Add braces in TSCAL and TZERO
+ * definitions
+ * Test for fill data
+ * 10/15/04 v1.23 bjg Fill Data Photon gets flag
+ * LOCATION_FILL only if it is
+ * in Active Aperture
+ * 12/24/2004 v1.24 bjg npix_max now takes binx and biny into
+ * account
+ * 01/24/2005 v1.25 wvd Modify cf_make_pixel_list to ignore
+ * duplicate data, which can sneak in
+ * with the stim pulses.
+ * 02/11/2005 v1.26 wvd Hot pixels and fill data cause the
+ * header keyword NEVENTS to be too
+ * large, so we correct/exclude these
+ * points and update NEVENTS.
+ * 03/10/2005 v1.27 wvd Move flagging of airglow photons to
+ * cf_screen_airglow.
+ * Initialize XFARF and YFARF with zero.
+ * 05/20/2005 v1.28 wvd Clean up i/o.
+ * 05/31/2006 v1.29 wvd OPUS has been updated, so we can set
+ * morekeys to 0.
+ * 06/12/2006 v1.30 wvd Call cf_verbose rather than
+ * cf_if_warning when an extension
+ * contains no photons.
+ * 02/23/2007 v1.31 wvd If a raw hist file contains no
+ * extensions, set EXP_STAT to -1.
+ * Make active-area limits inclusive.
+ * 11/07/2008 v1.32 wvd Clean up i/o.
+ *
+ ****************************************************************************/
+
+#include <stdlib.h>
+#include <string.h>
+#include <unistd.h>
+#include <math.h>
+
+#include "calfuse.h"
+
+struct key {
+ char keyword[FLEN_KEYWORD];
+ float value;
+ };
+
+
+/*****************************************************************************/
+
+static int
+cf_make_pixel_list(fitsfile *infits, fitsfile *outfits, char unbin) {
+
+/*****************************************************************************
+ *
+ * Description: procedure to take input data from the raw HIST data file and
+ * generate a table with 14 columns. This table lists
+ * all pixels that contain counts. The
+ * table is designed to mimic a TTAG list, so that it can be used
+ * with various modules in the TTAG pipeline
+ *
+ ****************************************************************************/
+
+ fitsfile *elecfits;
+ char fmt_byte[FLEN_CARD], fmt_float[FLEN_CARD], fmt_short[FLEN_CARD];
+ char *inpha, *dumarrb;
+ char elecfile[FLEN_CARD];
+ char keyword[FLEN_CARD], card[FLEN_CARD];
+ unsigned char *dumarrbu ;
+ int status=0, anynull=0, hdutype;
+ int colval;
+ int tfields = 14; /* output table will have 14 columns */
+ int active_l, active_r ;
+ long i, j, k, nevents, nfdp;
+ long frow=1, felem=1, nrows=1; /* output table will have 1 row */
+ float sum_weights;
+ float *intime, *dumarrf;
+
+ int sia_index[512], sia_temp[512]; /* Indices of SIA rectangles */
+ int sia_width = 2048, sia_height = 16; /* Dimensions of SIA recs */
+ int xx, yy; /* Pixel coordinates (unbinned) */
+
+ int x0, y0, xstart, ystart ;
+ int nx, ny, yndx, binx, biny ;
+ int nhot, nhot_total=0;
+ short *xpix, *ypix, *img, *fill_data_pix, nmask, nullval;
+ long npix, npix_max, npts ;
+ float *wtpix, exptime, mtime, nave;
+
+
+ char extname[]="TTAG DATA"; /* name of the extension containing
+ the binary table holding the data */
+
+ char *ttype[]={"TIME", "XRAW", "YRAW", "PHA", "WEIGHT", "XFARF",
+ "YFARF", "X", "Y", "CHANNEL", "TIMEFLGS",
+ "LOC_FLGS", "LAMBDA", "ERGCM2" };
+
+ char *tform[14]; /* will asign values when we find out
+ the number of elements in the data set */
+
+
+ char *tunit[]={"SECONDS", "PIXELS", "PIXELS", " ", " ", "PIXELS",
+ "PIXELS", "PIXELS", "PIXELS", "UNITLESS",
+ "UNITLESS", "UNITLESS", "ANGSTROMS","ERG CM^-2"};
+
+ struct key tscal[] = {{"TSCAL6", 0.25}, {"TSCAL7", 0.1},
+ {"TSCAL8", 0.25}, {"TSCAL9", 0.1}};
+ struct key tzero[] = {{"TZERO6", 8192.}, {"TZERO7", 0.},
+ {"TZERO8", 8192.}, {"TZERO9", 0.}};
+
+ /* Initialize the sia_index array. */
+ (void) memset((void *)sia_index, 0, sizeof(int)*512);
+
+ /*********************************************************************
+ Read information to be used in the screening analysis.
+ **********************************************************************/
+
+ /* Read the limits to the active area of the detector */
+ FITS_read_key(outfits, TSTRING, "ELEC_CAL", elecfile, NULL, &status);
+ FITS_open_file(&elecfits, cf_cal_file(elecfile), READONLY, &status);
+ FITS_read_key(elecfits, TINT, "ACTIVE_L", &active_l, NULL,
+ &status);
+ FITS_read_key(elecfits, TINT, "ACTIVE_R", &active_r, NULL,
+ &status);
+ FITS_close_file(elecfits, &status);
+ cf_verbose(3, "Limits to the active area: left=%d, right=%d",
+ active_l, active_r);
+
+ /* Determine the exposure time */
+ FITS_movabs_hdu(infits, 1, &hdutype, &status);
+ FITS_read_key(infits, TFLOAT, "EXPTIME", &exptime, NULL, &status);
+ cf_verbose(3, "exptime = %d ", cf_nint(exptime)) ;
+
+ /* Determine the binning factors in x and y */
+ FITS_read_key(infits, TINT, "SPECBINX", &binx, NULL, &status);
+ FITS_read_key(infits, TINT, "SPECBINY", &biny, NULL, &status);
+ cf_verbose(3, "binx=%d, biny=%d \n", binx, biny) ;
+
+
+ /***********************************************************************
+ Get information from the raw data file.
+ ************************************************************************/
+
+ /* Set up arrays to contain the positions of each pixel with counts
+ (xpix and ypix) and the number of counts in that pixel (wtpix) */
+ if (unbin) npix_max = (NXMAX * NYMAX) / binx ;
+ else npix_max = (NXMAX * NYMAX) / (binx * biny) ;
+ npix = -1 ;
+ xpix = (short *) cf_calloc(npix_max, sizeof(short)) ;
+ ypix = (short *) cf_calloc(npix_max, sizeof(short)) ;
+ wtpix = (float *) cf_calloc(npix_max, sizeof(float)) ;
+ fill_data_pix = (short *) cf_calloc(npix_max, sizeof(short)) ;
+
+ /*
+ * Loop over all extensions in the input file, and place contents
+ * of each extension into outbuff. Keep trying to read a new HDU
+ * in the input file until we hit the EOF.
+ *
+ * Note that the next line cannot use FITS_movrel_hdu because
+ * that routine automatically exits upon error.
+ */
+ for (i = 1; !(fits_movrel_hdu(infits, 1, &hdutype, &status)); i++) {
+ FITS_read_key(infits, TINT, "NAXIS1", &nx, NULL, &status);
+ FITS_read_key(infits, TINT, "NAXIS2", &ny, NULL, &status);
+ FITS_read_key(infits, TINT, "XORIGIN", &x0, NULL, &status);
+ FITS_read_key(infits, TINT, "YORIGIN", &y0, NULL, &status);
+
+ npts = nx * ny ;
+ cf_verbose(3, "Reading extension %d ", i) ;
+ cf_verbose(3, "nx=%d, ny=%d, npts=%d ", nx, ny, npts) ;
+
+ /* Initialize the sia_temp array */
+ (void) memset((void *)sia_temp, 0, sizeof(int)*512);
+
+ /*
+ * Allocate space for and read the input image
+ */
+ img = cf_malloc(sizeof(short) * npts);
+ FITS_read_img(infits, TSHORT, 1L, npts, &nullval,
+ img, &anynull, &status);
+
+ /* Specify the starting positions in x and y.
+ Place the point in the middle of the bin */
+ xstart = x0 * binx + binx / 2;
+ ystart = y0 * biny + biny / 2;
+ cf_verbose(3, "x0=%d, y0=%d, xstart=%d, ystart=%d",
+ x0, y0, xstart, ystart) ;
+
+ /* Go through the image and tabulate the pixels with counts. */
+ nhot=0 ;
+ for (j=0; j<npts; j++)
+ if (img[j] > 0) {
+
+ /* If this pixel falls in an SIA rectangle that has already
+ been written to the output data list, skip it and move to
+ the next SIA rectangle. */
+
+ xx = (x0 + (j % nx)) * binx;
+ yy = (y0 + (j / nx)) * biny;
+ k = xx / sia_width + yy / sia_height * 8;
+ if (k < 0 || k >= 512) continue;
+ if (sia_index[k] == 1) {
+ j += sia_width / binx - 1;
+ continue;
+ }
+ sia_temp[k] = 1;
+
+ /* Check for a hot pixel by comparing the value in the pixel with
+ those on either side. Ignore pixels outside of detector active
+ area or with fewer than 16 counts. */
+
+ if ((img[j] > 15) && (img[j] != FILL_DATA) &&
+ (xx >= active_l) && (xx <= active_r)) {
+ nave = (img[j-2] + img[j-1] + img[j+1] + img[j+2]) / 4. ;
+ if (nave > 32000.) nave=32000. ;
+ if (nave < 1.) nave=1. ;
+ if (img[j] > 8*nave) {
+ nhot++ ;
+ cf_verbose(3,"hot pixel at point %d: val=%d, ave=%f",
+ j, img[j], nave) ;
+ /* Repair hot pixel */
+ nmask=1 ;
+ while(nmask < 8*nave && nmask < 16000) nmask*=2 ;
+ nmask-=1 ;
+ img[j] = (img[j] & nmask) ;
+ cf_verbose(3,"corrected value = %d ",img[j]) ;
+ }
+ }
+
+ yndx = (int) ( j / nx ) ;
+
+ if (unbin){
+ for (k=0;k<biny;k++){
+ npix++ ;
+ xpix[npix] = (short) (xstart + (j - (yndx * nx)) * binx) ;
+ ypix[npix] = (short) (ystart + yndx*biny + k) ;
+ wtpix[npix] = (float) ( ((float) img[j]) / ((float) biny) ) ;
+
+ fill_data_pix[npix]=0;
+ if (img[j] == FILL_DATA) {
+ wtpix[npix]=0.0;
+ fill_data_pix[npix]=1;
+ }
+
+
+ }
+ }
+ else {
+ npix++ ;
+ xpix[npix] = (short) (xstart + (j - (yndx * nx)) * binx) ;
+ ypix[npix] = (short) (ystart + yndx*biny) ;
+ wtpix[npix] = (float) img[j];
+
+ fill_data_pix[npix]=0;
+ if (img[j] == FILL_DATA) {
+ wtpix[npix]=0.0;
+ fill_data_pix[npix]=1;
+ }
+ }
+
+
+ }
+
+ nhot_total += nhot;
+ cf_verbose(2, "%d hot pixels found in extn %d ",nhot, i) ;
+ cf_verbose(3, "total number of pixels with counts = %d \n", npix) ;
+ if (nhot > 10) cf_if_warning("More than 10 hot pixels found.") ;
+
+ free(img) ;
+
+ /* Copy new SIA indices to sia_index array. */
+ for (k = 0; k < 512; k++)
+ if (sia_temp[k]) sia_index[k] = 1;
+
+ } /* Finished loop over all extensions in infits */
+ /*
+ * Status upon completion of loop through extensions should be
+ * "end of file."
+ * If it is, reset value of status. If not, there was an
+ * unexpected error. Print and exit.
+ */
+ if( status == END_OF_FILE ) { /* !!! This might cause a problem */
+ status = 0;
+ } else {
+ cf_if_fits_error(status);
+ }
+
+ cf_verbose(3, "Finished reading data ") ;
+
+ if (npix < 0) {
+ npix = 1 ;
+ xpix[0]=0 ;
+ ypix[0]=0 ;
+ wtpix[0]=1. ;
+ fill_data_pix[0]=0;
+ }
+
+ /*************************************************************************
+ Output the analysis to a table
+ ***************************************************************************/
+
+ /* Generate the tform array */
+ sprintf(fmt_byte, "%ldB", npix);
+ sprintf(fmt_float, "%ldE", npix);
+ sprintf(fmt_short, "%ldI", npix);
+
+ tform[0] = fmt_float;
+ tform[1] = fmt_short;
+ tform[2] = fmt_short;
+ tform[3] = fmt_byte;
+ tform[4] = fmt_float;
+ for (i=5; i<9; i++)
+ tform[i] = fmt_short;
+ for ( ; i<12; i++)
+ tform[i] = fmt_byte;
+ tform[12] = fmt_float;
+ tform[13] = fmt_float;
+
+ /* Append a new empty binary table to the output file */
+
+ cf_verbose(3, "Creating table") ;
+ FITS_movabs_hdu(outfits, 1, &hdutype, &status);
+ FITS_create_tbl(outfits, BINARY_TBL, nrows, tfields, ttype, tform,
+ tunit, extname, &status);
+
+ /* Write TSCALE and TZERO entries to header */
+ for (i=0; i<4; i++) {
+ sprintf(keyword, "TUNIT%ld", i+6);
+ if (fits_read_keyword(outfits, keyword, card, NULL, &status))
+ cf_if_fits_error(status);
+ FITS_insert_key_flt(outfits, tscal[i].keyword, tscal[i].value, -2,
+ NULL, &status);
+ FITS_insert_key_flt(outfits, tzero[i].keyword, tzero[i].value, -4,
+ NULL, &status);
+ }
+
+ /* Allocate space for the new output columns */
+
+ cf_verbose(3, "Allocating space for arrays ") ;
+ dumarrf = (float *) cf_malloc(sizeof(float) * npix);
+ dumarrb = (char *) cf_malloc(sizeof(char) * npix);
+ dumarrbu = (unsigned char *) cf_malloc(sizeof(unsigned char) * npix);
+
+ cf_verbose(3, "Filling columns") ;
+ /* Specify a time array - assign the same time to all pixels */
+ intime = (float *) cf_calloc(npix, sizeof(float) ) ;
+ mtime = exptime/2. ;
+ for (i=0; i<npix; i++) intime[i]=mtime ;
+ FITS_write_col(outfits, TFLOAT, 1, frow, felem, npix, intime, &status);
+
+ /* Output X and Y positions */
+ FITS_write_col(outfits, TSHORT, 2, frow, felem, npix, xpix, &status);
+ FITS_write_col(outfits, TSHORT, 3, frow, felem, npix, ypix, &status);
+
+ /* Now do the PHA column - assume 20 */
+ inpha = (char *) cf_calloc(npix, sizeof(char) ) ;
+ for (i=0; i<npix; i++) inpha[i] = 20 ;
+ FITS_write_col(outfits, TBYTE, 4, frow, felem, npix, inpha, &status);
+
+ /* Fill weights column with counts per pixel */
+ FITS_write_col(outfits, TFLOAT, 5, frow, felem, npix, wtpix, &status);
+
+ /* Fill XFARF, YFARF, X and Y arrays with 0 */
+ for (i=0; i<npix; i++) dumarrf[i] = 0.;
+ for (colval=6; colval<10; colval++) {
+ FITS_write_col(outfits, TFLOAT, colval, frow, felem, npix,
+ dumarrf, &status);
+ }
+
+ /* Initialize the channel and temporal flag arrays with 0 - byte */
+ for (i=0;i<npix; i++)
+ dumarrb[i] = 0;
+ FITS_write_col(outfits, TBYTE, 10, frow, felem, npix, dumarrb, &status);
+ for (i=0;i<npix; i++)
+ dumarrbu[i] = 0;
+ FITS_write_col(outfits, TBYTE, 11, frow, felem, npix, dumarrbu, &status);
+
+ /* Flag events that lie outside of the active area or
+ in fill-data regions. */
+ for (i=0 ; i<npix ; i++) {
+ dumarrbu[i] = 0 ;
+
+ /* Flag photons that lie outside the active area. */
+ if (xpix[i] < active_l || xpix[i] > active_r)
+ dumarrbu[i] = dumarrbu[i] | LOCATION_SHLD ;
+
+ /* If inside, is it fill data? */
+ else if (fill_data_pix[i])
+ dumarrbu[i] = dumarrbu[i] | LOCATION_FILL ;
+ }
+
+ FITS_movabs_hdu(outfits, 2, &hdutype, &status);
+ FITS_write_col(outfits, TBYTE, 12, frow, felem, npix,
+ dumarrbu, &status);
+
+ /* Fill lambda and energy arrays with 0 */
+ FITS_write_col(outfits, TFLOAT, 13, frow, felem, npix, dumarrf, &status);
+ FITS_write_col(outfits, TFLOAT, 14, frow, felem, npix, dumarrf, &status);
+
+ /* If file contains hot pixels or fill data, recalculate NEVENTS. */
+ nfdp = 0;
+ sum_weights = 0.;
+ for (i=0; i<npix; i++) {
+ nfdp += fill_data_pix[i];
+ sum_weights += wtpix[i];
+ }
+ if (nfdp > 0 || nhot_total > 0) {
+ nevents = cf_nlong(sum_weights);
+ FITS_movabs_hdu(outfits, 1, &hdutype, &status);
+ FITS_update_key(outfits, TLONG, "NEVENTS", &nevents, NULL, &status);
+ cf_verbose(1, "Image contains bad pixels. Updating NEVENTS = %ld",
+ nevents);
+ }
+
+ free(intime);
+ free(xpix);
+ free(ypix);
+ free(wtpix) ;
+ free(inpha);
+ free(dumarrf);
+ free(dumarrb);
+
+ cf_verbose(3, "Leaving cf_make_pixel_list") ;
+
+ return EXIT_SUCCESS ;
+}
+
+
+/******************************************************************************
+ *
+ * CF_MAKE_GTI_TABLE
+ *
+ * Procedure to generate a table of GTI values and put it into the second
+ * extension of the output file. There is assumed to be only one GTI,
+ * starting at 0 and ending at the exposure time.
+ *
+ ******************************************************************************/
+
+int cf_make_gti_table (fitsfile *outfits) {
+
+ int status=0, nrows=1, hdutype, frow=1, felem=1 ;
+ float exptime ;
+ double start[1]={0.}, stop[1] ;
+
+ char extname[]="GTI";
+ int tfields=2; /* output table will have 2 columns */
+ char *ttype[]={"START", "STOP"};
+
+ char *tform[]={"1D","1D"};
+ char *tunit[]={"seconds", "seconds"};
+
+ cf_verbose(3, "Entering cf_make_gti_table ") ;
+
+ /* read the exposure time from the output file */
+ FITS_movabs_hdu(outfits, 1, &hdutype, &status);
+ FITS_read_key(outfits, TFLOAT, "EXPTIME", &exptime, NULL, &status);
+ stop[0] = (double) exptime ;
+
+
+ /* Append a new empty binary table to the output file */
+ FITS_movabs_hdu(outfits, 2, &hdutype, &status);
+ FITS_create_tbl(outfits, BINARY_TBL, nrows, tfields, ttype, tform,
+ tunit, extname, &status);
+
+ /* Write the data to the table */
+ FITS_write_col(outfits, TDOUBLE, 1, frow, felem, nrows, start, &status);
+ FITS_write_col(outfits, TDOUBLE, 2, frow, felem, nrows, stop, &status);
+
+ return EXIT_SUCCESS ;
+}
+
+/*****************************************************************************/
+
+
+int main(int argc, char *argv[])
+{
+ char CF_PRGM_ID[] = "cf_hist_init";
+ char CF_VER_NUM[] = "1.32";
+
+ fitsfile *infits, *outfits;
+
+ int morekeys=0; /* Change to zero when OPUS is updated. */
+ int status=0, optc;
+
+ int one=1, nextn=0;
+ long nphot ;
+ char unbin=0;
+
+ char opts[] = "hsv:";
+ char usage[] =
+ "Usage:\n"
+ " cf_hist_init [-h] [-s] [-v level] rawhistfile idf_file\n";
+ char option[] =
+ "Options:\n"
+ " -h: this help message\n"
+ " -s: to unbin the data in Y\n"
+ " -v: verbosity level (default is 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 's':
+ unbin=1;
+ break;
+ case 'v':
+ verbose_level = atoi(optarg);
+ break;
+ }
+ }
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+ if (argc != optind+2) {
+ printf("%s", usage);
+ cf_if_error("Incorrect number of arguments");
+ }
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin execution");
+
+ FITS_open_file(&infits, argv[optind++], READONLY, &status);
+
+ /* Exit if data are not in HIST mode. */
+ if ( cf_proc_check(infits, CF_PRGM_ID) ) {
+ FITS_close_file(infits, &status);
+ return (-1);
+ }
+
+ /* Issue warning message if data are missing. */
+ FITS_read_key(infits, TLONG, "NEVENTS", &nphot, NULL, &status);
+ FITS_read_key(infits, TINT, "NEXTEND", &nextn, NULL, &status);
+ cf_verbose(3,"Number of extensions to the file = %d ", nextn) ;
+ if (nphot < 1) cf_verbose(1, "%s contains no photons.", argv[optind-1]);
+ if (nextn < 1) cf_if_warning("%s contains no extensions.", argv[optind-1]);
+ else if (nextn < 4) cf_verbose(1, "%s contains only %d extensions.",
+ argv[optind-1], nextn) ;
+
+ /* Create output file. */
+ FITS_create_file(&outfits, argv[optind], &status);
+ /* Copy primary image header from the input file. */
+ FITS_copy_hdu(infits, outfits, morekeys, &status);
+
+ /* Update header keywords. */
+ FITS_update_key(outfits, TSTRING, "FILENAME", argv[optind], NULL, &status);
+ FITS_update_key(outfits, TSTRING, "FILETYPE",
+ "INTERMEDIATE DATA FILE", NULL, &status);
+
+ if (unbin) FITS_update_key(outfits, TINT, "SPECBINY", &one, NULL, &status);
+
+ /* If input file contains no extensions, set EXP_STAT to -1. */
+ if (nextn < 1) {
+ char comment[] = "Raw file contains no data.";
+ int eflag = -1;
+ FITS_update_key(outfits, TINT, "EXP_STAT", &eflag, comment, &status);
+ }
+
+ /* Populate the calibration keywords in the primary header. */
+ cf_fuv_init(outfits);
+
+ /* Fill in the first extension with the photon information */
+ cf_make_pixel_list(infits, outfits, unbin) ;
+
+ /* Generate a good-times table and put it in the second extension. */
+ cf_make_gti_table(outfits);
+
+ /* Append a table containing the orbital timeline */
+ cf_timeline(outfits);
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+
+ /* Update processing flags. */
+ cf_proc_update(outfits, CF_PRGM_ID, "COMPLETE");
+
+ FITS_close_file(infits, &status);
+ FITS_close_file(outfits, &status);
+
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Finished execution.");
+
+ return EXIT_SUCCESS ;
+}
diff --git a/src/fuv/cf_remove_motions.c b/src/fuv/cf_remove_motions.c
new file mode 100644
index 0000000..ec2b332
--- /dev/null
+++ b/src/fuv/cf_remove_motions.c
@@ -0,0 +1,316 @@
+/*****************************************************************************
+ * Johns Hopkins University
+ * Center For Astrophysical Sciences
+ * FUSE
+ *****************************************************************************
+ *
+ * Synopsis: cf_remove_motions options intermediate_file
+ *
+ * Description: Associates a photon with an aperture and corrects for grating,
+ * mirror, spacecraft (jitter) and focal plane assembly (FPA)
+ * motion. The position of the Y centroid is also determined
+ * both as a function of time (before correction for s/c motion)
+ * and for the entire exposure (after motion correction).
+ * The resulting coordinates represent the data obtained by a
+ * motionless instrument.
+ *
+ * By default, all corrections are performed. Command line
+ * options allow one or more corrections to be omitted.
+ *
+ * Arguments: input_file Intermediate data file
+ *
+ * Calibration files:
+ *
+ * Returns: 0 on successful completion
+ *
+ * Calls: cf_find_spectra, cf_calculate_ycent_motion,
+ * cf_grating_motion, cf_fpa_position, cf_mirror_motion,
+ * cf_satellite_jitter, cf_calculate_y_centroid
+ *
+ * History: 09/30/02 peb 1.1 Begin work
+ * 12/04/02 wvd 1.5 Install Y centroid routines.
+ * Modify call to cf_satellite_jitter
+ * 12/12/02 wvd 1.6 Change infits -> header throughout
+ * 03/04/03 peb 1.10 Minor code changes to remove gcc -Wall
+ * messages: unused j in main.
+ * Added unistd.h, so getopt is portable
+ * 03/05/03 wvd 1.11 Added cf_target_count_rate
+ * 03/10/03 peb 1.12 Added verbose_level, unistd.h for
+ * getopt portability, and changed
+ * locflags to unsigned char *.
+ * 03/13/03 rdr 1.13 corrected small error in tabulating
+ * count rate_sic
+ * 04/17/03 wvd 1.17 Add cf_find_spectra for initial
+ * assignment of photons to
+ * target apertures. Add last_call to
+ * cf_identify_channel, last_call and
+ * weight to cf_calculate_y_centroid
+ * 04/21/03 wvd 1.18 Pass tsunrise to cf_grating_motion,
+ * tsunset to cf_mirror_motion,
+ * channel to cf_satellite_jitter.
+ * 06/02/03 wvd 1.20 Don't modify count-rate arrays
+ * for HIST data.
+ * 06/11/03 wvd 1.21 Pass datatype to cf_read_col and
+ * cf_write_col.
+ * 08/01/03 wvd 1.23 Add cf_error_init after calls to
+ * subroutines.
+ * 08/06/03 wvd 1.24 Change channel to unsigned char.
+ * Remove GTI's from cf_satellite_jitter.
+ * 08/25/03 wvd 1.25 Change coltype from string to int in
+ * cf_read_col and cf_write_col.
+ * 09/18/03 wvd 1.26 Pass locflags to cf_identify_channel
+ * 10/21/03 wvd 1.27 For HIST data, don't write blank
+ * LiF and SiC count-rate arrays to IDF.
+ * 10/27/03 wvd 1.28 Pass pha array to cf_find_spectra.
+ * Delete last_call from argument list
+ * of cf_calculate_y_centroid.
+ * 10/31/03 wvd 1.29 Treat channel as unsigned char
+ * throughout.
+ * 04/06/04 bjg 1.30 Remove unused variables
+ * 05/03/04 wvd 1.31 Improve documentation.
+ * Move first call of cf_identify_channel
+ * from cf_find_spectra to this routine.
+ * 06/02/04 wvd 1.32 Copy/move cf_set_photon_flags(),
+ * cf_set_good_time_intervals(), and
+ * cf_modify_hist_times() from
+ * cf_screen_photons. Run them only if
+ * cf_satellite_jitter rejects some time
+ * AND if they were already run before.
+ * Pass timeflags to cf_find_spectra()
+ * and cf_calculate_y_centroid().
+ * Get exp_jitr from cf_satellite_jitter.
+ * 07/21/04 wvd 1.33 Add include file <string.h>
+ * 03/10/05 wvd 1.34 Change use of variable pad. Now add
+ * 10 to user_pad on first call to
+ * cf_identify_channel, 0 on second call.
+ * 03/22/05 wvd 1.35 Change TIME_SUNRISE and TIME_SUNSET
+ * from floats to shorts.
+ * 05/20/05 wvd 1.36 Clean up i/o.
+ * 09/09/05 wvd 1.37 Add -j to list of allowed options.
+ * 11/29/05 wvd 1.38 Separate cf_satellite_jitter into two
+ * routines. Screening is performed
+ * earlier, so we need not copy flags.
+ * 08/30/06 wvd 1.39 Add -a to list of allowed options.
+ * 12/29/06 wvd 1.40 Pass weights array and initial LiF
+ * and SiC counts arrays to
+ * cf_target_count_rate. Make arrays
+ * rate_lif and rate_sic floats. If
+ * they are too large to store as shorts,
+ * store as unsigned shorts.
+ * 02/02/07 wvd 1.41 If LiF or SiC count rate is too large
+ * for an unsigned int, set to zero.
+ *
+ ****************************************************************************/
+
+#include <unistd.h>
+#include <stdlib.h>
+#include <string.h>
+#include "calfuse.h"
+
+static char CF_PRGM_ID[] = "cf_remove_motions";
+static char CF_VER_NUM[] = "1.41";
+
+int
+main(int argc, char *argv[])
+{
+ unsigned char *channel=NULL, *timeflags=NULL, *locflags=NULL,
+ *statflag=NULL;
+ int grating_motion=1, fpa_position=1, mirror_motion=1, sat_jitter=1;
+ int last_call, set_times=FALSE, set_gtis=FALSE, airglow_centroid=FALSE;
+ int user_pad=0, status=0, set_rate=FALSE, optc;
+ long i, nevents, ntimes;
+ float *time=NULL, *ttime=NULL, *xfarf=NULL, *yfarf=NULL, *x, *y;
+ float *weight=NULL, *ycent_lif, *ycent_sic;
+ float max_rate_lif=0, max_rate_sic=0, *rate_lif=NULL, *rate_sic=NULL;
+ float tscale=1, tzero=32768;
+ short *tsunrise=NULL, *tsunset=NULL;
+ GTI gti;
+ fitsfile *header;
+
+ char opts[] = "hagfmjp:v:";
+ char usage[] =
+ "Usage:\n"
+ " cf_remove_motions [-hagfmj] [-p pad] [-v level] idf_file\n";
+ char option[] =
+ "Options:\n"
+ " -h: this help message\n"
+ " -v: verbose level (default is 1; 0 is silent)\n"
+ " -a: use airglow lines to determine spectral centroid\n"
+ " -p: aperture padding in pixels (default is 10)\n"
+ " -g: no grating motion correction\n"
+ " -f: no fpa position correction\n"
+ " -m: no mirror motion correction\n"
+ " -j: no satellite jitter correction\n";
+
+ verbose_level = 1;
+
+ while ((optc = getopt(argc, argv, opts)) != -1) {
+ switch(optc) {
+
+ case 'h':
+ printf("%s\n%s", usage, option);
+ return 0;
+ case 'v':
+ verbose_level = atoi(optarg);
+ break;
+ case 'p':
+ user_pad = atoi(optarg);
+ break;
+ case 'g':
+ grating_motion = 0;
+ break;
+ case 'f':
+ fpa_position = 0;
+ break;
+ case 'a':
+ airglow_centroid = TRUE;
+ break;
+ case 'm':
+ mirror_motion = 0;
+ break;
+ case 'j':
+ sat_jitter = 0;
+ break;
+ }
+ }
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+
+ if (argc <= optind) {
+ printf("%s", usage);
+ cf_if_error("Incorrect number of command-line arguments");
+ }
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");
+
+ FITS_open_file(&header, argv[optind], READWRITE, &status);
+
+ FITS_movabs_hdu(header, 2, NULL, &status);
+ nevents = cf_read_col(header, TFLOAT, "TIME", (void **) &time);
+ nevents = cf_read_col(header, TFLOAT, "WEIGHT", (void **) &weight);
+ nevents = cf_read_col(header, TFLOAT, "XFARF", (void **) &xfarf);
+ nevents = cf_read_col(header, TFLOAT, "YFARF", (void **) &yfarf);
+ nevents = cf_read_col(header, TBYTE, "TIMEFLGS", (void **) &timeflags);
+ nevents = cf_read_col(header, TBYTE, "LOC_FLGS", (void **) &locflags);
+ channel = cf_calloc(nevents, sizeof(unsigned char));
+
+ FITS_movabs_hdu(header, 4, NULL, &status) ;
+ ntimes = cf_read_col(header, TFLOAT, "TIME", (void **) &ttime) ;
+ ntimes = cf_read_col(header, TBYTE, "STATUS_FLAGS", (void **) &statflag);
+ ntimes = cf_read_col(header, TSHORT, "TIME_SUNRISE", (void **) &tsunrise) ;
+ ntimes = cf_read_col(header, TSHORT, "TIME_SUNSET", (void **) &tsunset) ;
+ ntimes = cf_read_col(header, TFLOAT, "LIF_CNT_RATE", (void **) &rate_lif) ;
+ ntimes = cf_read_col(header, TFLOAT, "SIC_CNT_RATE", (void **) &rate_sic) ;
+ ycent_lif=(float *) cf_calloc(ntimes, sizeof(float)) ;
+ ycent_sic=(float *) cf_calloc(ntimes, sizeof(float)) ;
+
+ FITS_movabs_hdu(header, 1, NULL, &status);
+ x = xfarf;
+ y = yfarf;
+
+ /* Find all six spectra and compute centroids. */
+ cf_find_spectra(header, nevents, weight, x, y, channel, timeflags,
+ locflags, airglow_centroid);
+
+ /* Initial assignment of channel numbers to photons */
+ cf_identify_channel(header, nevents, xfarf, yfarf, channel, locflags,
+ user_pad+10, last_call=FALSE);
+
+ /* Track Y centroids of moving spectra. */
+ cf_calculate_ycent_motion(header, nevents, time, y, channel, locflags,
+ ntimes, ttime, ycent_lif, ycent_sic);
+
+ /* Correct for detector motion during exposure. */
+ if (grating_motion)
+ cf_grating_motion(header, nevents, time, x, y, channel,
+ ntimes, ttime, tsunrise);
+
+ /* Shift spectra in X to account for position of the FPA. */
+ if (fpa_position)
+ cf_fpa_position(header, nevents, x, channel);
+
+ /* Correct for mirror motion during exposure. */
+ if (mirror_motion)
+ cf_mirror_motion(header, nevents, time, x, y, channel,
+ ntimes, ttime, tsunset);
+
+ /* Correct for spacecraft motion during exposure. */
+ if (sat_jitter)
+ cf_satellite_jitter(header, nevents, time, x, y, channel,
+ ntimes, ttime, statflag);
+
+ /* Compute centroids of motion-corrected spectra. */
+ cf_calculate_y_centroid(header, nevents, weight, x, y, channel,
+ timeflags, locflags);
+
+ /* Final assignment of channel numbers to photons */
+ cf_identify_channel(header, nevents, x, y, channel, locflags,
+ user_pad, last_call=TRUE);
+
+ /* Compute count rates for TTAG target spectra. */
+ set_rate = !(cf_target_count_rate(header, nevents, time, weight,
+ channel, locflags, ntimes, ttime, rate_lif, rate_sic));
+
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+
+ cf_verbose(3, "Writing X, Y, CHANNEL to IDF");
+ FITS_movabs_hdu(header, 2, NULL, &status);
+ if (set_times)
+ cf_write_col(header, TFLOAT, "TIME", (void *) time, nevents);
+ cf_write_col(header, TFLOAT, "X", (void *) x, nevents);
+ cf_write_col(header, TFLOAT, "Y", (void *) y, nevents);
+ cf_write_col(header, TBYTE, "CHANNEL", (void *) channel, nevents);
+
+ if (set_gtis) {
+ cf_verbose(3, "Writing good-time intervals to IDF");
+ FITS_movabs_hdu(header, 3, NULL, &status);
+ cf_write_col(header, TDOUBLE, "START", (void *) gti.start, gti.ntimes);
+ cf_write_col(header, TDOUBLE, "STOP", (void *) gti.stop, gti.ntimes);
+ free(gti.stop);
+ free(gti.start);
+ }
+
+ FITS_movabs_hdu(header, 4, NULL, &status);
+ cf_verbose(3, "Writing timeline information to IDF");
+ if (set_rate) {
+ for (i = 0; i < ntimes; i++) {
+ if (rate_lif[i] > 65535) rate_lif[i] = 0;
+ if (rate_sic[i] > 65535) rate_sic[i] = 0;
+ if (max_rate_lif < rate_lif[i]) max_rate_lif = rate_lif[i];
+ if (max_rate_sic < rate_sic[i]) max_rate_sic = rate_sic[i];
+ }
+ cf_verbose(2, "max_rate_lif = %.0f", max_rate_lif);
+ cf_verbose(2, "max_rate_sic = %.0f", max_rate_sic);
+ if (max_rate_lif > 32767) {
+ FITS_update_key(header, TFLOAT, "TSCAL10", &tscale, NULL, &status);
+ FITS_update_key(header, TFLOAT, "TZERO10", &tzero, NULL, &status);
+ fits_set_tscale(header, 10, (double)tscale, (double)tzero, &status);
+ }
+ if (max_rate_sic > 32767) {
+ FITS_update_key(header, TFLOAT, "TSCAL11", &tscale, NULL, &status);
+ FITS_update_key(header, TFLOAT, "TZERO11", &tzero, NULL, &status);
+ fits_set_tscale(header, 11, (double)tscale, (double)tzero, &status);
+ }
+ FITS_write_col(header, TFLOAT, 10, 1, 1, ntimes, rate_lif, &status);
+ FITS_write_col(header, TFLOAT, 11, 1, 1, ntimes, rate_sic, &status);
+ }
+ cf_write_col(header, TFLOAT, "YCENT_LIF", (void *) ycent_lif, ntimes);
+ cf_write_col(header, TFLOAT, "YCENT_SIC", (void *) ycent_sic, ntimes);
+ cf_verbose(3, "IDF updates complete");
+
+ FITS_close_file(header, &status);
+
+ free(time);
+ free(weight);
+ free(xfarf);
+ free(yfarf);
+ free(locflags);
+ free(timeflags);
+ free(channel);
+ free(ttime);
+ free(rate_lif);
+ free(rate_sic);
+ free(ycent_lif);
+ free(ycent_sic);
+
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done processing");
+ return status;
+}
diff --git a/src/fuv/cf_screen_photons.c b/src/fuv/cf_screen_photons.c
new file mode 100644
index 0000000..bd82bba
--- /dev/null
+++ b/src/fuv/cf_screen_photons.c
@@ -0,0 +1,298 @@
+/*****************************************************************************
+ * Johns Hopkins University
+ * Center For Astrophysical Sciences
+ * FUSE
+ *****************************************************************************
+ *
+ * Synopsis: cf_screen_photons options intermediate_file
+ *
+ * Description: Filters events in the intermediate data file (IDF) for limb
+ * angle, SAA crossing, high voltage changes, bursts, and PHA.
+ * Sets the event flags and GTIs. Data are modified in place.
+ *
+ * By default all corrections are performed. Command line
+ * options allow one or more corrections to be omitted.
+ *
+ * Arguments: input_file FARF-corrected intermediate data file
+ *
+ * Calibration files:
+ *
+ * Returns: 0 on successful completion
+ *
+ * Calls: cf_screen_limb_angle, cf_screen_saa, cf_screen_high_voltage,
+ * cf_screen_burst, cf_screen_jitter,
+ * cf_set_user_gtis, cf_set_photon_flags, cf_modify_hist_times,
+ * cf_screen_pulse_height, cf_screen_airglow, cf_screen_bad_pixels
+ *
+ * History: 11/05/02 1.1 peb Begin work
+ * 11/13/02 1.2 peb Added screening functions to program
+ * 11/14/02 1.3 peb Added burst-screening function
+ * 11/18/02 1.4 peb Corrected logic dealing with reading,
+ * writing, and freeing of GTI struct.
+ * Corrected some grammatical errors.
+ * 12/10/02 1.5 rdr changed call to cf_screen_bursts
+ * 12/18/02 1.6 rdr updated names of some columns in
+ * the calls to the timeline extension
+ * 12/20/02 1.7 wvd Change flags to unsigned char
+ * 03/10/03 1.9 peb Added verbose_level and added unistd.h
+ * for getopt portability
+ * 05/10/03 1.10 wvd Pass locflag to cf_set_photon_flags
+ * 05/30/03 1.11 wvd Pass weight to cf_set_photon_flags
+ * 06/11/03 1.12 wvd Treat HV as array of shorts.
+ * Pass datatype to cf_read_col and
+ * cf_write_col.
+ * 08/25/03 1.13 wvd Change coltype from string to int in
+ * cf_read_col and cf_write_col.
+ * 09/10/03 1.14 wvd Change background array to type short.
+ * Add cf_set_user_gtis().
+ * 10/02/03 1.15 wvd Pass locflag array to
+ * cf_screen_pulse_height; write to IDF
+ * 02/10/04 1.16 wvd Add cf_screen_fifo_overflow()
+ * 06/02/04 1.17 wvd Add cf_modify_hist_times()
+ * 06/03/04 1.18 wvd Read XFARF and YFARF, not X and Y.
+ * If INSTMODE = TTAG, don't bother to
+ * call cf_modify_hist_times().
+ * 02/02/05 1.19 wvd Pass rate_aic to cf_screen_burst
+ * 03/10/05 1.20 wvd Add cf_screen_airglow()
+ * 03/30/05 1.21 wvd Don't change TIME_COR to SKIPPED
+ * for TTAG data. Leave as OMIT.
+ * 05/20/05 1.22 wvd Clean up i/o.
+ * 06/03/05 1.23 wvd Fix typo in CF_VER_NUM.
+ * 09/09/05 1.24 wvd Update list of allowed options.
+ * 11/09/05 1.25 wvd Change argument list for
+ * cf_screen_fifo_overflow.
+ * 11/22/05 1.26 wvd Add cf_screen_jitter and
+ * cf_screen_bad_pixels.
+ * 01/24/06 1.27 wvd Move cf_screen_airglow before
+ * cf_screen_burst
+ * 12/29/06 1.28 wvd Move cf_screen_fifo_overflow to
+ * cf_convert_to_farf.
+ * 07/18/08 1.29 wvd Write STATUS_FLAGS to IDF before
+ * calling cf_set_photon_flags, which
+ * may change the array.
+ *
+ ****************************************************************************/
+
+#include <unistd.h>
+#include <stdlib.h>
+#include <string.h>
+#include "calfuse.h"
+
+int
+main(int argc, char *argv[])
+{
+ char CF_PRGM_ID[] = "cf_screen_photons";
+ char CF_VER_NUM[] = "1.29";
+
+ char instmode[FLEN_VALUE];
+ unsigned char *pha=NULL, *timeflag=NULL, *statflag=NULL, *locflag=NULL;
+ int limb_screening=1, saa_screening=1, voltage_screening=1;
+ int burst_screening=1, set_flags=1, airglow_screening=1, pha_screening=1;
+ int modify_times=1, set_gtis=1, set_user_gtis=1;
+ int bad_pixel_screening=1, set_times=FALSE, jitter_screening=1;
+ int status=0, optc;
+ long nevents, nseconds;
+ short *background=NULL, *rate_lif=NULL, *rate_sic=NULL, *voltage=NULL;
+ float *time=NULL, *x=NULL, *y=NULL, *weight=NULL;
+ float *timeline=NULL, *limb=NULL, *longitude=NULL, *latitude=NULL;
+ float *rate_aic=NULL;
+ GTI gti;
+ fitsfile *header;
+
+ char opts[] = "halstbjufimdpv:";
+ char usage[] =
+ "Usage:\n"
+ " cf_screen_photons [-halstbjufimdp] [-v level] idf_file\n";
+ char option[] =
+ "Options:\n"
+ " -h: this help message\n"
+ " -v: verbosity level (default is 1; 0 is silent)\n"
+ " -a: no airglow screening\n"
+ " -l: no limb screening\n"
+ " -s: no SAA screening\n"
+ " -t: no high voltage screening\n"
+ " -b: no burst screening\n"
+ " -j: no jitter screening\n"
+ " -u: no user-defined good-time intervals\n"
+ " -f: do not set photon event flags\n"
+ " -i: do not update GTI arrays\n"
+ " -m: do not modify HIST photon times\n"
+ " -d: no bad-pixel screening\n"
+ " -p: no PHA screening\n";
+
+ verbose_level = 1;
+
+ while ((optc = getopt(argc, argv, opts)) != -1) {
+ switch(optc) {
+
+ case 'h':
+ printf("%s\n%s", usage, option);
+ return 0;
+ case 'v':
+ verbose_level = atoi(optarg);
+ break;
+ case 'l':
+ limb_screening = 0;
+ break;
+ case 't':
+ voltage_screening = 0;
+ break;
+ case 's':
+ saa_screening = 0;
+ break;
+ case 'b':
+ burst_screening = 0;
+ break;
+ case 'j':
+ jitter_screening = 0;
+ break;
+ case 'u':
+ set_user_gtis = 0;
+ break;
+ case 'f':
+ set_flags = 0;
+ break;
+ case 'i':
+ set_gtis = 0;
+ break;
+ case 'm':
+ modify_times = 0;
+ break;
+ case 'a':
+ airglow_screening = 0;
+ break;
+ case 'd':
+ bad_pixel_screening = 0;
+ break;
+ case 'p':
+ pha_screening = 0;
+ break;
+ }
+ }
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+ if (argc <= optind) {
+ printf("%s", usage);
+ cf_if_error("Incorrect number of program arguments");
+ }
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");
+
+ FITS_open_file(&header, argv[optind], READWRITE, &status);
+
+ FITS_movabs_hdu(header, 2, NULL, &status);
+ nevents = cf_read_col(header, TFLOAT, "TIME", (void **) &time);
+ nevents = cf_read_col(header, TBYTE, "PHA", (void **) &pha);
+ nevents = cf_read_col(header, TFLOAT, "WEIGHT", (void **) &weight);
+ nevents = cf_read_col(header, TFLOAT, "XFARF", (void **) &x);
+ nevents = cf_read_col(header, TFLOAT, "YFARF", (void **) &y);
+ nevents = cf_read_col(header, TBYTE, "TIMEFLGS", (void **) &timeflag);
+ nevents = cf_read_col(header, TBYTE, "LOC_FLGS", (void **) &locflag);
+
+ if (burst_screening) {
+ FITS_movabs_hdu(header, 3, NULL, &status);
+ gti.ntimes = cf_read_col(header, TDOUBLE, "START",(void **) &gti.start);
+ gti.ntimes = cf_read_col(header, TDOUBLE, "STOP", (void **) &gti.stop);
+ }
+ FITS_movabs_hdu(header, 4, NULL, &status);
+ nseconds = cf_read_col(header, TFLOAT, "TIME", (void **) &timeline);
+ nseconds = cf_read_col(header, TBYTE, "STATUS_FLAGS", (void **) &statflag);
+ nseconds = cf_read_col(header, TFLOAT, "LIMB_ANGLE", (void **) &limb);
+ nseconds = cf_read_col(header, TFLOAT, "LONGITUDE", (void **) &longitude);
+ nseconds = cf_read_col(header, TFLOAT, "LATITUDE", (void **) &latitude);
+ nseconds = cf_read_col(header, TSHORT, "HIGH_VOLTAGE", (void **) &voltage);
+ nseconds = cf_read_col(header, TSHORT, "LIF_CNT_RATE", (void **) &rate_lif);
+ nseconds = cf_read_col(header, TSHORT, "SIC_CNT_RATE", (void **) &rate_sic);
+ nseconds = cf_read_col(header, TFLOAT, "AIC_CNT_RATE", (void **) &rate_aic);
+ background = (short *) cf_calloc(nseconds, sizeof(short)) ;
+
+ /* If INSTMODE = TTAG, don't bother to call cf_modify_hist_times. */
+ FITS_movabs_hdu(header, 1, NULL, &status);
+ FITS_read_key(header, TSTRING, "INSTMODE", instmode, NULL, &status);
+ if (!strncmp(instmode, "TTAG", 4)) modify_times = FALSE;
+
+ if (airglow_screening)
+ cf_screen_airglow(header, nevents, x, y, locflag);
+
+ if (limb_screening)
+ cf_screen_limb_angle(header, nseconds, statflag, limb);
+
+ if (saa_screening)
+ cf_screen_saa(header, nseconds, statflag, longitude, latitude);
+
+ if (voltage_screening)
+ cf_screen_high_voltage(header, nseconds, statflag, voltage);
+
+ if (burst_screening) {
+ cf_screen_burst(header, nevents, time, x, y, locflag, &gti,
+ nseconds, timeline, statflag, rate_aic, background);
+ free(gti.stop);
+ free(gti.start);
+ }
+
+ if (jitter_screening)
+ cf_screen_jitter(header, nseconds, timeline, statflag);
+
+ if (set_user_gtis)
+ cf_set_user_gtis(header, nseconds, timeline, statflag);
+
+ /* Write STATUS_FLAGS array to IDF before calling cf_set_photon_flags. */
+ FITS_movabs_hdu(header, 4, NULL, &status);
+ cf_write_col(header, TBYTE, "STATUS_FLAGS", (void *) statflag, nseconds);
+ FITS_movabs_hdu(header, 1, NULL, &status);
+
+ if (set_flags)
+ cf_set_photon_flags(header, nevents, time, weight, timeflag,
+ locflag, nseconds, timeline, statflag);
+
+ if (set_gtis)
+ cf_set_good_time_intervals(header, nseconds, timeline,
+ statflag, &gti);
+
+ if (modify_times)
+ set_times = !(cf_modify_hist_times(header, nevents, time, &gti));
+
+ if (bad_pixel_screening)
+ cf_screen_bad_pixels(header, nevents, x, y, locflag);
+
+ if (pha_screening)
+ cf_screen_pulse_height(header, nevents, pha, locflag);
+
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+
+ FITS_movabs_hdu(header, 2, NULL, &status);
+ if (set_times)
+ cf_write_col(header, TFLOAT, "TIME", (void *) time, nevents);
+ cf_write_col(header, TBYTE, "TIMEFLGS", (void *) timeflag, nevents);
+ cf_write_col(header, TBYTE, "LOC_FLGS", (void *) locflag, nevents);
+
+ if (set_gtis) {
+ FITS_movabs_hdu(header, 3, NULL, &status);
+ cf_write_col(header, TDOUBLE, "START", (void *) gti.start, gti.ntimes);
+ cf_write_col(header, TDOUBLE, "STOP", (void *) gti.stop, gti.ntimes);
+ free(gti.stop);
+ free(gti.start);
+ }
+
+ if (burst_screening) {
+ FITS_movabs_hdu(header, 4, NULL, &status);
+ cf_write_col(header, TSHORT, "BKGD_CNT_RATE", (void *) background,
+ nseconds);
+ }
+
+ FITS_close_file(header, &status);
+
+ free(background);
+ free(voltage);
+ free(latitude);
+ free(longitude);
+ free(limb);
+ free(statflag);
+ free(timeline);
+ free(timeflag);
+ free(pha);
+ free(y);
+ free(x);
+ free(time);
+
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done processing");
+ return 0;
+}
diff --git a/src/fuv/cf_ttag_init.c b/src/fuv/cf_ttag_init.c
new file mode 100644
index 0000000..ec11f19
--- /dev/null
+++ b/src/fuv/cf_ttag_init.c
@@ -0,0 +1,594 @@
+/*****************************************************************************
+ * Johns Hopkins University
+ * Center For Astrophysical Sciences
+ * FUSE
+ *****************************************************************************
+ *
+ * Synopsis: cf_ttag_init input_file output_file
+ *
+ * Description: Reads a raw ttag file, which is a FITS binary table
+ * consisting of four columns:
+ *
+ * TIME FLOAT
+ * X INT
+ * Y INT
+ * PHA BYTE
+ *
+ * This routine will convert the raw ttag file to another binary
+ * table with 13 columns the values for most of these columns
+ * will be added by various modules of the pipeline:
+ *
+ * TIME FLOAT time of arrival for photon
+ * XRAW INT raw x position
+ * YRAW INT raw y position
+ * PHA BYTE pha value
+ * WEIGHT FLOAT effective area and flux calibration
+ * XFARF FLOAT x position in geometrically
+ * corrected frame
+ * YFARF FLOAT y position in geometrically
+ * corrected frame
+ * X FLOAT final x position (after correcting
+ * for motions)
+ * Y FLOAT final y position (after correcting
+ * for motions)
+ * CHANNEL CHAR aperture ID for the photon
+ * TIMEFLGS CHAR temporal flags - various bits set
+ * LOC_FLGS CHAR location flags - various bits set
+ * LAMBDA FLOAT wavelength for each photon
+ * ERGCM2 FLOAT flux calibration for each photon
+ *
+ * This routine calls cf_fuv_init to set up the calibration
+ * entries in the header.
+ *
+ * In addition, the routine creates a timeline containing
+ * orbital information, count rates, and high-voltage
+ * information which is used in later pipeline modules.
+ *
+ * Turns out that OPUS good-time intervals may exclude SAA
+ * passes that remain in the photon-event list. This means
+ * that there may be photons in the photon-event list whose
+ * times are not listed in the timeline table. To catch
+ * them, we set all photon status flags to TEMPORAL_HV.
+ *
+ *
+ * Returns: Zero upon successful completion
+ *
+ * History: 03/30/98 gak Begin work
+ * 03/18/99 emm Modified header.
+ * 04/06/99 emm Added check on nrows to be no
+ * larger than nevents.
+ * 04/15/99 barrett Deleted FF and DY columns,
+ * added FITS wrapper functions,
+ * tidied code.
+ * 04/22/99 emm Program now returns int.
+ * 04/23/99 emm Fixed two bugs causing segmentation
+ * faults (pha was being written as
+ * long, at definition of outdx was
+ * given as 1B instead of 1E.) Added
+ * check on number of arguments.
+ * 04/29/99 emm Added check on nrows to make sure
+ * it is smaller than nevents.
+ * 04/30/99 peb Added cf_malloc functions.
+ * 06/07/99 peb Added reporting of version number.
+ * 18/07/99 1.12 jm Populate APER_ACT keyword
+ * 10/28/99 1.14 emm Added DY column.
+ * 12/02/99 1.15 emm Added corrections to some keywords
+ * removed APER_ACT keyword update
+ * 12/21/00 1.16 peb Fixed type mismatches in cfitsio
+ * calls
+ * 05/16/01 1.17 rdr added night/day flag column (dnflg)
+ * 08/20/01 1.18 wvd Check keywords SPECBINX and SPECBINY;
+ * if 0, set to 1
+ * 10/19/01 1.19 wvd Move correction of SPECBINX, etc.
+ * from cf_ttag_image.
+ * 04/03/02 1.20 wvd Clean up keywords in photon-event list.
+ * 08/02/02 1.21 rdr - changed the format of the binary
+ * tables to improve speed
+ * - added new columns to the table
+ * 09/15/02 1.22 rdr added extension containing timeline
+ * 11/07/02 1.23 peb Modified columns for event and
+ * timeline flags. Removed temporary
+ * keywords. Changed include files to
+ * use calfusettag.h
+ * 11/13/02 1.24 rdr Added ycentl and ycents columns to
+ * timeline.
+ * Added section to flag geocoronal photons
+ * and photons in the inactive part of the
+ * detector. Moved processing of the
+ * photon list into a separate subroutine.
+ * Fixed problem of program crashing when
+ * needed columns are not in the hskp file
+ * 12/06/02 1.25 rdr removed variables which were not being
+ * used
+ * 12/12/02 1.5 wvd Installed program.
+ * 12/18/02 1.6 rdr use AIRG calibration file for airglow
+ * locations
+ * 12/20/02 1.7 rdr - check for negative count rates and
+ * correct
+ * - convert status flags to unsigned char
+ * - flag photons near start and end of
+ * GTI
+ * 01/24/03 1.8 rdr - correct situation where hskp file
+ * has blank columns
+ * 02/11/03 1.9 rdr Added ERGCM2S column to photon table
+ * 03/04/03 1.10 rdr removed YFARFC column and deleted some
+ * unused variables
+ * 03/19/03 1.12 peb Added command-line options and
+ * verbose_level parameter
+ * Changed all printf calls to cf_verbose
+ * calls with verbose level == 1
+ * 03/25/03 1.13 peb Changed a cf_verbose to cf_timestamp.
+ * 04/01/03 1.15 wvd Call cf_error_init after cf_timeline.
+ * 04/04/03 1.16 wvd Update FILETYPE keyword. Set verbose
+ * level of each call to cf_verbose.
+ * 04/07/03 1.17 wvd Update FILENAME keyword.
+ * 04/07/03 1.18 wvd Add FPADXLIF and FPADXSIC to header.
+ * 04/09/03 1.19 rdr Correct error in determining the name
+ * of the housekeeping file
+ * 04/29/03 1.21 wvd Allow possiblity that housekeeping
+ * filename is in lower-case letters.
+ * 05/06/03 1.22 rdr Add daynight keyword, bkg_max and
+ * bkg_min keywords, and changed the
+ * column ergcm2s to ergcm2
+ * 05/10/03 1.23 wvd Modify internal copy of GTI's.
+ * Make number of timeline entries more
+ * nearly equal to EXPTIME. Discard
+ * photons with bad time markers.
+ * 06/03/03 1.24 wvd Add call to cf_proc_check.
+ * Keep trying to get GTI's right.
+ * 06/09/03 1.25 wvd Compress data with TSCALE and TZERO.
+ * 06/11/03 1.26 wvd Change calfusettag.h to calfuse.h
+ * 07/07/03 1.27 rdr Correct initialization of dn_flag_last
+ * when searching for time_sunrise and
+ * time_sunset
+ * 07/17/03 1.28 wvd Move various subroutines to library.
+ * 07/23/03 1.29 wvd Close ELEC_CAL file.
+ * Delete goto statement.
+ * 09/15/03 1.31 wvd Set temporal flags to TEMPORAL_HV.
+ * 10/02/03 1.32 wvd Delete reference to LOCATION_GTI flag
+ * 10/26/03 1.33 wvd Increase to morekeys to 38.
+ * 11/05/03 1.34 wvd Added stdlib.h and unistd.h
+ * 03/03/04 1.35 rdr Exit if no photons found.
+ * 04/06/04 1.37 bjg Change formats to match arg types
+ * in printf
+ * Functions return EXIT_SUCCESS
+ * Remove static keyword in struct key
+ * definition
+ * Add braces in TSCAL and TZERO
+ * definitions
+ * When no events are present in
+ * the raw file, create a photon list
+ * with a dummy event.
+ * 03/10/05 1.38 wvd Move airglow screening to
+ * cf_screen_airglow.
+ * Initialize XFARF and YFARF to zero.
+ * 05/20/05 1.39 wvd Clean up i/o.
+ * 05/31/06 1.40 wvd OPUS has been updated, so we can set
+ * morekeys to 0.
+ * 06/12/06 1.41 wvd Change cf_if_warning to cf_verbose.
+ * 08/21/07 1.42 wvd Add -e flag to set HKEXISTS to YES.
+ *
+ ****************************************************************************/
+
+#include <stdlib.h>
+#include <unistd.h>
+#include "calfuse.h"
+
+struct key {
+ char keyword[FLEN_KEYWORD];
+ float value;
+ };
+
+static int
+cf_init_null_photon_list(fitsfile *outfits) {
+
+ int status=0;
+
+ long i;
+
+ char keyword[FLEN_CARD], card[FLEN_CARD];
+
+ char extname[]="TTAG DATA"; /* Name of this extension */
+
+ char *ttype[]={"TIME", "XRAW", "YRAW", "PHA", "WEIGHT", "XFARF",
+ "YFARF", "X", "Y", "CHANNEL", "TIMEFLGS",
+ "LOC_FLGS", "LAMBDA", "ERGCM2" };
+
+ char *tform[14]={"1E","1I","1I","1B","1E","1I","1I",
+ "1I","1I","1B","1B","1B","1E","1E"};
+
+ char *tunit[]={"SECONDS", "PIXELS", "PIXELS", "UNITLESS", "UNITLESS",
+ "PIXELS", "PIXELS", "PIXELS", "PIXELS", "UNITLESS",
+ "UNITLESS", "UNITLESS", "ANGSTROMS", "ERG CM^-2"};
+
+ struct key tscal[] = {{"TSCAL6", 0.25}, {"TSCAL7", 0.1},
+ {"TSCAL8", 0.25}, {"TSCAL9", 0.1}};
+ struct key tzero[] = {{"TZERO6", 8192.}, {"TZERO7", 0.},
+ {"TZERO8", 8192.}, {"TZERO9", 0.}};
+
+ float hdu2_time[1] = {0.0};
+ short hdu2_xraw[1] = {0};
+ short hdu2_yraw[1] = {0};
+ char hdu2_pha[1] = {0};
+ float hdu2_weight[1] = {0.0};
+ float hdu2_xfarf[1] = {0.0};
+ float hdu2_yfarf[1] = {0.0};
+ float hdu2_x[1] = {0.0};
+ float hdu2_y[1] = {0.0};
+ char hdu2_channel[1] = {0};
+ char hdu2_timeflgs[1] = {TEMPORAL_HV};
+ char hdu2_loc_flgs[1] = {LOCATION_SHLD};
+ float hdu2_lambda[1] = {0.0};
+ float hdu2_ergcm2s[1] = {0.0};
+
+ cf_verbose(3, "Entering cf_init_null_photon_list.");
+
+ /* Append a new empty binary table to the output file */
+ FITS_create_tbl(outfits, BINARY_TBL, 1, 14, ttype, tform,
+ tunit, extname, &status);
+
+ /* Write TSCALE and TZERO entries to header */
+ for (i=0; i<4; i++) {
+ sprintf(keyword, "TUNIT%ld", i+6);
+ if (fits_read_keyword(outfits, keyword, card, NULL, &status))
+ cf_if_fits_error(status);
+ FITS_insert_key_flt(outfits, tscal[i].keyword, tscal[i].value,
+ -2, NULL, &status);
+ FITS_insert_key_flt(outfits, tzero[i].keyword, tzero[i].value,
+ -4, NULL, &status);
+ }
+
+
+
+
+ FITS_write_col(outfits, TFLOAT, 1, 1, 1, 1, hdu2_time , &status);
+ FITS_write_col(outfits, TSHORT, 2, 1, 1, 1, hdu2_xraw , &status);
+ FITS_write_col(outfits, TSHORT, 3, 1, 1, 1, hdu2_yraw , &status);
+ FITS_write_col(outfits, TBYTE , 4, 1, 1, 1, hdu2_pha , &status);
+ FITS_write_col(outfits, TFLOAT, 5, 1, 1, 1, hdu2_weight , &status);
+ FITS_write_col(outfits, TFLOAT, 6, 1, 1, 1, hdu2_xfarf , &status);
+ FITS_write_col(outfits, TFLOAT, 7, 1, 1, 1, hdu2_yfarf , &status);
+ FITS_write_col(outfits, TFLOAT, 8, 1, 1, 1, hdu2_x , &status);
+ FITS_write_col(outfits, TFLOAT, 9, 1, 1, 1, hdu2_y , &status);
+ FITS_write_col(outfits, TBYTE , 10, 1, 1, 1, hdu2_channel , &status);
+ FITS_write_col(outfits, TBYTE , 11, 1, 1, 1, hdu2_timeflgs , &status);
+ FITS_write_col(outfits, TBYTE , 12, 1, 1, 1, hdu2_loc_flgs , &status);
+ FITS_write_col(outfits, TFLOAT, 13, 1, 1, 1, hdu2_lambda , &status);
+ FITS_write_col(outfits, TFLOAT, 14, 1, 1, 1, hdu2_ergcm2s , &status);
+
+ cf_verbose(3, "Exiting cf_init_null_photon_list.");
+ return EXIT_SUCCESS;
+}
+
+/*****************************************************************************/
+
+static int
+cf_init_photon_list(fitsfile *infits, fitsfile *outfits) {
+
+/*****************************************************************************
+ *
+ * Description: procedure to take input data from the raw data file and
+ * generate a table with 14 columns. These columns
+ * will be modified by later pipeline modules.
+ *
+ * ARGUMENTS:
+ * infits fitsfile name of the input raw data file
+ * outfits fitsfile name of the output intermediate data
+ * file
+ * INPUT DATA:
+ * TIME FLOAT time of arrival for each photon
+ * XRAW INT x position of photon
+ * YRAW INT y position of photon
+ * PHA BYTE pha value of the photon
+ *
+ * OUTPUT DATA:
+ * TIME FLOAT time of arrival for photon
+ * XRAW INT raw x position
+ * YRAW INT raw y position
+ * PHA BYTE pha value
+ * WEIGHT FLOAT effective area and flux calibration
+ * XFARF FLOAT x position in geometrically
+ * corrected frame
+ * YFARF FLOAT y position in geometrically
+ * corrected frame
+ * X FLOAT final x position (after correcting
+ * for motions)
+ * Y FLOAT final y position (after correcting
+ * for motions)
+ * CHANNEL CHAR aperture ID for the photon
+ * TIMEFLGS CHAR temporal flags - various bits set
+ * LOC_FLGS CHAR location flags - various bits set
+ * LAMBDA FLOAT wavelength for each photon
+ * ERGCM2 FLOAT flux cal for each photon
+ *
+ ****************************************************************************/
+
+ fitsfile *elecfits;
+ char fmt_byte[FLEN_CARD], fmt_float[FLEN_CARD], fmt_short[FLEN_CARD];
+ char *inpha;
+ char elecfile[FLEN_FILENAME];
+ char keyword[FLEN_CARD], card[FLEN_CARD];
+ unsigned char *channel, *location, *temporal;
+ short *inx, *iny ;
+ int status=0, intnull=0, anynull, hdutype;
+ int tcol, xcol, ycol, phacol, colval;
+ int tfields = 14; /* output table will have 14 columns */
+ int active_l, active_r;
+ long i;
+ long frow=1, felem=1, nevents, nrows=1; /* output table will have 1 row */
+ float *intime, *dumarrf;
+
+ char extname[]="TTAG DATA"; /* Name of this extension */
+
+ char *ttype[]={"TIME", "XRAW", "YRAW", "PHA", "WEIGHT", "XFARF",
+ "YFARF", "X", "Y", "CHANNEL", "TIMEFLGS",
+ "LOC_FLGS", "LAMBDA", "ERGCM2" };
+
+ char *tform[14]; /* We'll assign values when we know
+ the number of elements in the data set. */
+
+
+ char *tunit[]={"SECONDS", "PIXELS", "PIXELS", "UNITLESS", "UNITLESS",
+ "PIXELS", "PIXELS", "PIXELS", "PIXELS", "UNITLESS",
+ "UNITLESS", "UNITLESS", "ANGSTROMS", "ERG CM^-2"};
+
+ struct key tscal[] = {{"TSCAL6", 0.25}, {"TSCAL7", 0.1},
+ {"TSCAL8", 0.25}, {"TSCAL9", 0.1}};
+ struct key tzero[] = {{"TZERO6", 8192.}, {"TZERO7", 0.},
+ {"TZERO8", 8192.}, {"TZERO9", 0.}};
+
+ cf_verbose(3, "Entering cf_init_photon_list.");
+
+ /* Move to the extension of the input file containing the data. */
+ FITS_movabs_hdu(infits, 2, &hdutype, &status);
+
+ /* Read the limits to the active area of the detector. */
+ FITS_read_key(outfits, TSTRING, "ELEC_CAL", elecfile, NULL, &status);
+ FITS_open_file(&elecfits, cf_cal_file(elecfile), READONLY, &status);
+ FITS_read_key(elecfits, TINT, "ACTIVE_L", &active_l, NULL, &status);
+ FITS_read_key(elecfits, TINT, "ACTIVE_R", &active_r, NULL, &status);
+ FITS_close_file(elecfits, &status);
+ cf_verbose(3, "Limits to the active area: left=%d, right=%d",
+ active_l, active_r);
+
+ /*
+ * Get the number of events in the input file.
+ */
+ FITS_read_key(infits, TLONG, "NAXIS2", &nevents, NULL, &status);
+ cf_verbose(3, "Raw data file contains %ld entries", nevents);
+
+ /* Generate the tform array */
+ sprintf(fmt_byte, "%ldB", nevents);
+ sprintf(fmt_float, "%ldE", nevents);
+ sprintf(fmt_short, "%ldI", nevents);
+
+ tform[0] = fmt_float;
+ tform[1] = fmt_short;
+ tform[2] = fmt_short;
+ tform[3] = fmt_byte;
+ tform[4] = fmt_float;
+ for (i=5; i<9; i++)
+ tform[i] = fmt_short;
+ for ( ; i<12; i++)
+ tform[i] = fmt_byte;
+ tform[12] = fmt_float;
+ tform[13] = fmt_float;
+
+ /* Append a new empty binary table to the output file */
+ FITS_create_tbl(outfits, BINARY_TBL, nrows, tfields, ttype, tform,
+ tunit, extname, &status);
+
+ /* Write TSCALE and TZERO entries to header */
+ for (i=0; i<4; i++) {
+ sprintf(keyword, "TUNIT%ld", i+6);
+ if (fits_read_keyword(outfits, keyword, card, NULL, &status))
+ cf_if_fits_error(status);
+ FITS_insert_key_flt(outfits, tscal[i].keyword, tscal[i].value,
+ -2, NULL, &status);
+ FITS_insert_key_flt(outfits, tzero[i].keyword, tzero[i].value,
+ -4, NULL, &status);
+ }
+
+ /*
+ * Allocate space for the input times, coords, and pha.
+ */
+ intime = (float *) cf_malloc(sizeof(float) * nevents);
+ inx = (short *) cf_malloc(sizeof(short) * nevents);
+ iny = (short *) cf_malloc(sizeof(short) * nevents);
+ inpha = (char *) cf_malloc(sizeof(char) * nevents);
+ /*
+ * Allocate space for the new output columns.
+ */
+ dumarrf = (float *) cf_malloc(sizeof(float) * nevents);
+ channel = (unsigned char *) cf_calloc(nevents, sizeof(unsigned char));
+ location = (unsigned char *) cf_calloc(nevents, sizeof(unsigned char));
+ temporal = (unsigned char *) cf_calloc(nevents, sizeof(unsigned char));
+
+ /*
+ * Read the data from the input file and write to the output
+ * file in a new format
+ *
+ * read the times
+ */
+ FITS_get_colnum(infits, TRUE, "TIME", &tcol, &status);
+ FITS_read_col(infits, TFLOAT, tcol, frow, felem, nevents, &intnull,
+ intime, &anynull, &status);
+ FITS_write_col(outfits, TFLOAT, 1, frow, felem, nevents, intime, &status);
+
+ /* process X and Y positions */
+ FITS_get_colnum(infits, TRUE, "X", &xcol, &status);
+ FITS_read_col(infits, TSHORT, xcol, frow, felem, nevents, &intnull,
+ inx, &anynull, &status);
+ FITS_write_col(outfits, TSHORT, 2, frow, felem, nevents, inx, &status);
+
+ FITS_get_colnum(infits, TRUE, "Y", &ycol, &status);
+ FITS_read_col(infits, TSHORT, ycol, frow, felem, nevents, &intnull,
+ iny, &anynull, &status);
+ FITS_write_col(outfits, TSHORT, 3, frow, felem, nevents, iny, &status);
+
+ /* now do the PHA column */
+ FITS_get_colnum(infits, TRUE, "PHA", &phacol, &status);
+ FITS_read_col(infits, TBYTE, phacol, frow, felem, nevents, &intnull,
+ inpha, &anynull, &status);
+ FITS_write_col(outfits, TBYTE, 4, frow, felem, nevents, inpha, &status);
+
+ /* fill weights with a dummy float - all equal to 1. */
+ for (i=0; i<nevents; i++)
+ dumarrf[i] = 1.;
+ FITS_write_col(outfits, TFLOAT, 5, frow, felem, nevents, dumarrf, &status);
+
+ /* fill XFARF, YFARF, X and Y arrays with 0. */
+ for (i=0; i<nevents; i++)
+ dumarrf[i] = 0.;
+ for (colval=6; colval<10; colval++) {
+ FITS_write_col(outfits, TFLOAT, colval, frow, felem, nevents,
+ dumarrf, &status);
+ }
+
+ /* Default value of temporal array is TEMPORAL_HV */
+ for (i=0; i<nevents; i++)
+ temporal[i] = TEMPORAL_HV;
+ FITS_write_col(outfits, TBYTE, 10, frow, felem, nevents, channel, &status);
+ FITS_write_col(outfits, TBYTE, 11, frow, felem, nevents, temporal, &status);
+
+ /* Set flags in the screening array for photons outside of the active
+ region of the detector. */
+ for (i=0; i<nevents; i++)
+ if (inx[i] < active_l || inx[i] > active_r)
+ location[i] = location[i] | LOCATION_SHLD;
+
+ /* Write location flags to output file. */
+ FITS_movabs_hdu(outfits, 2, &hdutype, &status);
+ FITS_write_col(outfits, TBYTE, 12, frow, felem, nevents,
+ location, &status);
+
+ /* Fill lambda and energy arrays with 0. */
+ FITS_write_col(outfits, TFLOAT, 13, frow, felem, nevents, dumarrf,
+ &status);
+ FITS_write_col(outfits, TFLOAT, 14, frow, felem, nevents, dumarrf,
+ &status);
+
+ free(intime);
+ free(inx);
+ free(iny);
+ free(inpha);
+ free(dumarrf);
+ free(channel);
+ free(location);
+ free(temporal);
+
+ cf_verbose(3, "Exiting cf_init_photon_list.");
+ return EXIT_SUCCESS;
+}
+
+
+/*****************************************************************************/
+
+int main(int argc, char *argv[])
+
+/*****************************************************************************/
+{
+ char CF_PRGM_ID[] = "cf_ttag_init";
+ char CF_VER_NUM[] = "1.42";
+
+ fitsfile *infits, *outfits;
+ int morekeys=0;
+ int modify_hkexists=FALSE;
+ int status=0, hdutype, optc;
+ long nphot ;
+
+ char opts[] = "ehv:";
+ char usage[] =
+ "Usage:\n"
+ " cf_ttag_init [-he] [-v level] rawttagfile idf_file\n";
+ char option[] =
+ "Options:\n"
+ " -h: this help message\n"
+ " -e: set header keyword HKEXISTS to YES in IDF file\n"
+ " -v: verbosity level (default is 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 'e':
+ modify_hkexists = TRUE;
+ break;
+ case 'v':
+ verbose_level = atoi(optarg);
+ break;
+ }
+ }
+
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+ if (argc != optind+2) {
+ printf("%s", usage);
+ cf_if_error("Incorrect number of arguments");
+ }
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin execution");
+
+ FITS_open_file(&infits, argv[optind++], READONLY, &status);
+
+ /* Exit if data are not in TTAG mode. */
+ if ( cf_proc_check(infits, CF_PRGM_ID) ) {
+ FITS_close_file(infits, &status);
+ return (-1);
+ }
+
+ /* Issue warning if there are no events in the data set. */
+ FITS_read_key(infits, TLONG, "NEVENTS", &nphot, NULL, &status);
+ cf_verbose(3,"Number of events = %d ",nphot) ;
+ if (nphot < 1) {
+ cf_verbose(1, "No photons found in the data.") ;
+ }
+
+
+ /* Create output file. */
+ FITS_create_file(&outfits, argv[optind], &status);
+ /* Copy primary image header from the input file. */
+ FITS_copy_hdu(infits, outfits, morekeys, &status);
+
+ /* Update header keywords. */
+ FITS_update_key(outfits, TSTRING, "FILENAME", argv[optind], NULL, &status);
+ FITS_update_key(outfits, TSTRING, "FILETYPE",
+ "INTERMEDIATE DATA FILE", NULL, &status);
+ if (modify_hkexists)
+ FITS_update_key(outfits, TSTRING, "HKEXISTS", "YES", NULL, &status);
+
+ /* Populate the calibration keywords in the primary header. */
+ cf_fuv_init(outfits);
+
+ /* Write background limits to file header */
+ cf_set_background_limits(outfits);
+
+ /* Fill in the first extension with the photon information */
+ if (nphot<1) {
+ cf_init_null_photon_list(outfits);
+ }
+ else {
+ cf_init_photon_list(infits, outfits);
+ }
+
+ /* Copy good-time intervals from input to output. */
+ cf_verbose(3, "Copying good-time intervals from input to output.");
+ FITS_movabs_hdu(infits, 3, &hdutype, &status);
+ FITS_create_hdu(outfits, &status);
+ FITS_copy_hdu(infits, outfits, 0, &status);
+
+ /* Append a table containing the orbital timeline */
+ cf_timeline(outfits);
+ cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
+
+ /* Update processing flags. */
+ cf_proc_update(outfits, CF_PRGM_ID, "COMPLETE");
+
+ FITS_close_file(infits, &status);
+ FITS_close_file(outfits, &status);
+
+ cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Finished execution.");
+ return EXIT_SUCCESS;
+}
generated by cgit (git 2.34.1) at 2025-09-20 13:03:36 -0400