Initial commit

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;
+}