From d54fe7c1f704a63824c5bfa0ece65245572e9b27 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 4 Mar 2015 21:21:30 -0500 Subject: Initial commit --- src/analysis/idf_combine.c | 945 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 945 insertions(+) create mode 100644 src/analysis/idf_combine.c (limited to 'src/analysis/idf_combine.c') diff --git a/src/analysis/idf_combine.c b/src/analysis/idf_combine.c new file mode 100644 index 0000000..3f1a678 --- /dev/null +++ b/src/analysis/idf_combine.c @@ -0,0 +1,945 @@ + +/***************************************************************************** + * Johns Hopkins University + * Center For Astrophysical Sciences + * FUSE + ***************************************************************************** + * + * Synopsis: idf_combine outfile file1 file2 file3 file4 file5 ... + * + * Description: Reads multiple contiguous idf files and concatenates them + * into a single file. + * + * History: 08/11/03 bjg Begin work from ttag_combine.c + * sort the input files in time order + * 08/12/03 bjg v1.0 First usable version + * 08/13/03 bjg Some minor changes + * 08/14/03 bjg Fixed some important keywords + * in the main HDU + * 08/20/03 bjg v1.1 Changed the way the binary tables + * are read, created and written + * (same as cf_ttag_init) + * 09/09/03 bjg v1.2 Fixed some keywords + * 10/29/03 bjg Changed the scaling format for + * output timeline to conform the + * new IDF formatting. + * 11/17/03 bjg v1.3 Now warns users if the Focal + * Plane Assembly has moved between + * exposures. + * 11/22/03 bjg Now prints the name of the + * included files in the main + * header. + * 12/10/03 bjg Updated NSPEC keyword and + * SPECxxx, WOFFLxxx, WOFFSxxx keywords + * Removed path in filenames written + * into header. + * 02/25/04 bjg Bug fix + * 03/25/04 bjg Added TZERO and TSCALE values for + * AIC and FEC countrates + * 04/05/04 bjg Remove unused variables + * Change formats to match arg types + * in printf + * Remove static keyword in struct key + * definition + * Add braces in TSCAL and TZERO + * definitions + * 06/08/04 bjg Changes to handle EXP_STAT keyword + * 08/13/04 bjg v1.4 Recalculates Y centroid on demand + * 10/06/04 bjg v1.5 Added option -b to store + * ORBITAL_VEL in a float + * 10/26/04 bjg v1.6 Cosmetic change + * 01/28/2005 bjg v1.7 Fixed crash within CFITSIO with + * with IDF files having empty tables + * 02/16/2005 wvd v1.8 Double TSCALE and TZERO + * for the AIC_CNT_RATE array. + * 03/22/2005 wvd v1.9 Read and write TIME_SUNRISE + * and TIME_SUNSET as shorts. + * 04/15/2005 tc v1.10 Update PLANTIME keyword + * 06/08/2005 tc v1.11 Comment out warning if FPA has + * shifted. Pipeline corrects for this. + * 09/19/2005 wvd v1.12 Don't even read FPA positions. + * 05/22/2006 wvd v1.13 Change most header keyword + * types from float to long. + * 06/22/2006 wvd v1.14 Add -z flag to force creation of + * an output file, even if empty. + * 11/10/2006 wvd v1.15 Add comment fields to WOFFS and + * WOFFL keywords. + * 08/27/2007 bot v1.16 Changed long for int l.138 + * 07/25/2008 wvd v1.17 Set output EXP_STAT keyword to lowest + * non-negative value among input files. + * 08/22/2008 wvd v1.18 If -a flag is set, include files with + * positive values of EXP_STAT, but not + * those with negative values. + * + ****************************************************************************/ + + +#include +#include +#include +#include +#include +#include "calfuse.h" + +typedef char filename[FLEN_CARD]; + +struct key { + char keyword[FLEN_KEYWORD]; + float value; +}; + +static char CF_PRGM_ID[]= "idf_combine"; +static char CF_VER_NUM[]= "1.18"; + +int main(int argc,char *argv[]){ + + char date[FLEN_CARD]={'\0'}; + char rootname[FLEN_CARD]; + char comment[FLEN_COMMENT]; + + int felem_hdu2,felem_hdu4,frow_hdu3; + + char stime[FLEN_CARD],keyword[FLEN_CARD],card[FLEN_CARD]; + + char fmt_byte[FLEN_CARD],fmt_float[FLEN_CARD],fmt_short[FLEN_CARD]; + + char *corrected_filename; + char *string_pointer; + char string0[100], string1[100]; + + time_t vtime; + + fitsfile *infits,*outfits; + + filename *filelist; + double *expstartlist; + double *expendlist; + long *neventslist; + long *nrecordslist; + long *ngtislist; + + double delta_t; + + float zero=0.0; + + filename tempstring; + double tempdouble; + long templong; + float tempfloat; + char tempchar; + + double minexpstart; + int minindex; + + int nfiles; + + int intnull=0,anynull; + int ncol; + + int status=0; + int hdutype=0; + int tref = 0; + + int n2; + + int exp_stat, exp_stat_out=99; + + long nevents=0, nrecords=0, ngtis=0; + long n_real_events=0; + + long i,j; + + float *floatlist; + short *shortlist; + char *charlist; + double *doublelist; + + float exptime=0, rawtime=0; + long neventscreened=0, neventscreenedpha=0, plantime=0, timehv=0; + long timescreened=0, timesaa=0, timelowlimbangle=0, timeburst=0, + timejitter=0,timenight=0; + + /* float fpasx0, fpasx, fpalx0, fpalx, dfpasx, dfpalx ; */ + /* int fpa_split = FALSE; */ + + int hdu2_tfields=14; + + char hdu2_extname[]="TTAG DATA"; /* Name of this extension */ + + char *hdu2_ttype[]={"TIME", "XRAW", "YRAW", "PHA", "WEIGHT", "XFARF", + "YFARF", "X", "Y", "CHANNEL", "TIMEFLGS", + "LOC_FLGS", "LAMBDA", "ERGCM2" }; + + char *hdu2_tform[14]; /* We'll assign values when we know + the number of elements in the data set. */ + + + char *hdu2_tunit[]={"SECONDS", "PIXELS", "PIXELS", "UNITLESS", "UNITLESS", + "PIXELS", "PIXELS", "PIXELS", "PIXELS", "UNITLESS", + "UNITLESS", "UNITLESS", "ANGSTROMS", "ERG CM^-2"}; + + struct key hdu2_tscal[] = {{"TSCAL6", 0.25}, {"TSCAL7", 0.1}, + {"TSCAL8", 0.25}, {"TSCAL9", 0.1}}; + + struct key hdu2_tzero[] = {{"TZERO6", 8192.}, {"TZERO7", 0.}, + {"TZERO8", 8192.}, {"TZERO9", 0.}}; + + + + + int hdu3_tfields=2; + + char hdu3_extname[]="GTI"; /* Name of this extension */ + + char *hdu3_ttype[]={"START", "STOP" }; + + char *hdu3_tform[2]={"1D", "1D" }; /* We'll assign values when we know + the number of elements in the data set. */ + + + char *hdu3_tunit[]={"seconds", "seconds"}; + + struct key hdu4_tscal[16]; + struct key hdu4_tzero[16]; + + char hdu4_extname[]="TIMELINE"; + int hdu4_tfields=16; /* output table will have 16 columns */ + char *hdu4_ttype[]={"TIME", "STATUS_FLAGS", "TIME_SUNRISE", "TIME_SUNSET", + "LIMB_ANGLE", "LONGITUDE", "LATITUDE", "ORBITAL_VEL", + "HIGH_VOLTAGE", "LIF_CNT_RATE", "SIC_CNT_RATE", + "FEC_CNT_RATE", "AIC_CNT_RATE", "BKGD_CNT_RATE", + "YCENT_LIF","YCENT_SIC"}; + + char *hdu4_tform[16]; /* we will define tform later, when the number + of photons is known */ + + char *hdu4_tunit[]={"seconds", "unitless", "seconds", "seconds", "degrees", + "degrees", "degrees", "km/s", "unitless", "counts/sec", + "counts/sec", "counts/sec", "counts/sec", "counts/sec", + "pixels","pixels" }; + + + int ignore_exp_stat=0; + int do_ycent=0; + int big_vel=0; + int do_empty=FALSE, empty_output=FALSE; + + float *x=NULL, *y=NULL, *weight=NULL; + unsigned char *channel=NULL, *timeflags=NULL, *locflags=NULL; + + int optc; + + char opts[] = "hacbv:z"; + char usage[] = + "Usage:\n" + " idf_combine [-ahbcz] [-v level] output_idf_file input_idf_files\n"; + char option[] = + "Options:\n" + " -h: this help message\n" + " -v: verbosity level (=1; 0 is silent)\n" + " -a: ignore EXP_STAT keyword \n" + " -c: recalculates Y centroids (use target events) \n" + " -b: store ORBITAL_VEL in a float \n" + " -z: if no good data, generate empty output file\n"; + + verbose_level = 1; + + /* 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 'a': + ignore_exp_stat=1; + break; + case 'v': + verbose_level = atoi(optarg); + break; + case 'c': + do_ycent=1; + break; + case 'b': + big_vel=1; + break; + case 'z': + do_empty=TRUE; + break; + } + } + + cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr); + + if (argc < optind+2) { + printf("%s", usage); + return 1; + } + + + cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Started execution."); + + /* get and display time */ + vtime = time(NULL) ; + strcpy(stime,ctime(&vtime)); + + + nfiles=argc-optind-1; + + filelist = (filename *)cf_calloc(nfiles, sizeof(filename)); + expstartlist = (double *)cf_calloc(nfiles, sizeof(double)); + expendlist = (double *)cf_calloc(nfiles, sizeof(double)); + neventslist = (long *)cf_calloc(nfiles, sizeof(long)); + nrecordslist = (long *)cf_calloc(nfiles, sizeof(long)); + ngtislist = (long *)cf_calloc(nfiles, sizeof(long)); + + if (!ignore_exp_stat) + cf_verbose(1,"Will include only files with EXP_STAT=0") ; + else + cf_verbose(1,"Will include only files with EXP_STAT >= 0") ; + + cf_verbose(1,"GETTING INFORMATION ON INPUT FILES") ; + + n2=0; + + for (i=0; i exp_stat && exp_stat >= 0) exp_stat_out = exp_stat; + + /* Skip files with bad values of EXP_STAT. */ + if (!ignore_exp_stat && exp_stat != 0) { + cf_if_warning("File %s rejected. EXP_STAT not equal to zero. Use -a flag to override.", + argv[optind+i+1]) ; + FITS_close_file(infits,&status); + continue; + } + /* Never include files with negative values of EXP_STAT. */ + else if (exp_stat < 0) { + cf_if_warning("File %s rejected. EXP_STAT less than zero.", argv[optind+i+1]) ; + FITS_close_file(infits,&status); + continue; + } + + FITS_read_key(infits,TFLOAT,"EXPTIME",&(tempfloat),NULL,&status); + exptime+=tempfloat; + FITS_read_key(infits,TLONG,"NEVENTS",&templong,NULL,&status); + n_real_events+=templong; + FITS_read_key(infits,TFLOAT,"RAWTIME",&(tempfloat),NULL,&status); + rawtime+=tempfloat; + FITS_read_key(infits,TLONG,"PLANTIME",&(templong),NULL,&status); + plantime+=tempfloat; + + FITS_read_key(infits,TLONG,"NBADEVNT",&templong,NULL,&status); + neventscreened+=templong; + FITS_read_key(infits,TLONG,"NBADPHA",&templong,NULL,&status); + neventscreenedpha+=templong; + + FITS_read_key(infits,TLONG,"EXP_BAD",&templong,NULL,&status); + timescreened+=templong; + FITS_read_key(infits,TLONG,"EXP_BRST",&templong,NULL,&status); + timeburst+=templong; + FITS_read_key(infits,TLONG,"EXP_HV",&templong,NULL,&status); + timehv+=templong; + FITS_read_key(infits,TLONG,"EXP_JITR",&templong,NULL,&status); + timejitter+=templong; + FITS_read_key(infits,TLONG,"EXP_LIM",&templong,NULL,&status); + timelowlimbangle+=templong; + FITS_read_key(infits,TLONG,"EXP_SAA",&templong,NULL,&status); + timesaa+=templong; + FITS_read_key(infits,TLONG,"EXPNIGHT",&templong,NULL,&status); + timenight+=templong; + + FITS_movabs_hdu(infits,2,&hdutype,&status); + FITS_read_key(infits,TSTRING,"TFORM1",&tempstring,NULL,&status); + sscanf(tempstring,"%ld%c",&neventslist[n2],&tempchar); + nevents+=neventslist[n2]; + + FITS_movabs_hdu(infits,3,&hdutype,&status); + FITS_read_key(infits,TLONG,"NAXIS2",&(ngtislist[n2]),NULL,&status); + ngtis+=ngtislist[n2]; + + FITS_movabs_hdu(infits,4,&hdutype,&status); + FITS_read_key(infits,TSTRING,"TFORM1",&tempstring,NULL,&status); + sscanf(tempstring,"%ld%c",&nrecordslist[n2],&tempchar); + nrecords+=nrecordslist[n2]; + + FITS_close_file(infits,&status); + + n2++; + } + + nfiles=n2; + + if (nfiles==0){ + if (do_empty) { + cf_verbose(1,"No files to combine. Generating empty output file."); + do_ycent = FALSE; + empty_output = TRUE; + expendlist[0]=expstartlist[0]; + nfiles = 1; + } + else { + cf_verbose(1,"No files to combine. Exiting."); + cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Finished execution."); + exit(0); + } + } + + cf_verbose(1,"SORTING INPUT FILES IN TIME ORDER") ; + + for (i=0; i0){ + floatlist=(float *)malloc(neventslist[i]*sizeof(float)); + shortlist=(short *)malloc(neventslist[i]*sizeof(short)); + charlist=(char *)malloc(neventslist[i]*sizeof(char)); + + FITS_get_colnum(infits, TRUE, "TIME", &ncol, &status); + FITS_read_col(infits, TFLOAT, ncol, 1, 1, neventslist[i], &intnull, + floatlist, &anynull, &status); + for (j=0;j0){ + doublelist=(double *)malloc(ngtislist[i]*sizeof(double)); + + FITS_get_colnum(infits, TRUE, "START", &ncol, &status); + FITS_read_col(infits, TDOUBLE, ncol, 1, 1, ngtislist[i], &intnull, + doublelist, &anynull, &status); + for (j=0;j0){ + floatlist=(float *)malloc(nrecordslist[i]*sizeof(float)); + shortlist=(short *)malloc(nrecordslist[i]*sizeof(short)); + charlist=(char *)malloc(nrecordslist[i]*sizeof(char)); + + + FITS_get_colnum(infits, TRUE, "TIME", &ncol, &status); + FITS_read_col(infits, TFLOAT, ncol, 1, 1, nrecordslist[i], &intnull, + floatlist, &anynull, &status); + for (j=0;j