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