aboutsummaryrefslogtreecommitdiff
path: root/src/fuv/cf_remove_motions.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/fuv/cf_remove_motions.c')
-rw-r--r--src/fuv/cf_remove_motions.c316
1 files changed, 316 insertions, 0 deletions
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;
+}
generated by cgit (git 2.34.1) at 2025-09-20 05:11:45 -0400