summaryrefslogtreecommitdiffstats
path: root/funtools/funcnts.c
diff options
context:
space:
mode:
Diffstat (limited to 'funtools/funcnts.c')
-rw-r--r--funtools/funcnts.c1693
1 files changed, 1693 insertions, 0 deletions
diff --git a/funtools/funcnts.c b/funtools/funcnts.c
new file mode 100644
index 0000000..6da09ed
--- /dev/null
+++ b/funtools/funcnts.c
@@ -0,0 +1,1693 @@
+/*
+ * Copyright (c) 1999-2003 Smithsonian Astrophysical Observatory
+ */
+
+#include <math.h>
+#include <funtools.h>
+#include <NaN.h>
+#include <filter.h>
+#include <swap.h>
+#include <word.h>
+#include <xalloc.h>
+
+#define MAXROW 8192
+static int maxrow=MAXROW;
+
+/* max number of intervals we will process in one pass through table data */
+#define MAXINTV 1024
+
+/* type of data we process */
+#define SRC 0
+#define BKG 1
+/* this must match number of types above */
+#define NTYPE 2
+
+/* how corrections (time and exposure) are applied to source/bkgd */
+#define NONE 0
+#define SCORR 1
+#define BCORR 2
+#define BOTH 3
+
+/* types of background */
+#define BKG_VAL 1
+#define BKG_ALL 2
+#define BKG_EACH 3
+
+/* types of image -> exposure pixel conversion algorithms */
+#define EXP_RATIO 1
+#define EXP_WCS 2
+
+#define ARCSEC_PER_DEG 3600.0
+#define ARCSEC_PER_DEGSQ (ARCSEC_PER_DEG*ARCSEC_PER_DEG)
+
+extern char *optarg;
+extern int optind;
+
+#ifdef ANSI_FUNC
+static void
+usage (char *fname)
+#else
+static void usage(fname)
+ char *fname;
+#endif
+{
+ fprintf(stderr,
+ "usage: %s <switches> sname [sreg] [bname breg|breg|cnts]\n",
+ fname);
+ fprintf(stderr, "optional switches:\n");
+ fprintf(stderr, " -e \"source_exposure[;bkgd_exposure]\" # exp matches data\n");
+ fprintf(stderr, " -w \"source_exposure[;bkgd_exposure]\" # WCS method\n");
+ fprintf(stderr, "\t\t# source (bkgd) FITS exposure image\n");
+ fprintf(stderr, " -t \"source_timecorr[;bkgd_timecorr]\"\n");
+ fprintf(stderr, "\t\t# source (bkgd) time correction value or header parameter name\n");
+ fprintf(stderr, " -g\t\t# output using nice g format\n");
+ fprintf(stderr, " -G\t\t# output using %%.14g format (maximum precision)\n");
+ fprintf(stderr, " -i \"[column;]int1;int2...\" # column-based intervals\n");
+ fprintf(stderr, " -m\t\t# match individual source and bkgd regions\n");
+ fprintf(stderr, " -p\t\t# output in pixels, even if wcs is present\n");
+ fprintf(stderr, " -r\t\t# output inner/outer radii (and angles) for annuli (and pandas)\n");
+ fprintf(stderr, " -s\t\t# output summed values\n");
+ fprintf(stderr, " -v \"scol[;bcol]\" # src and bkgd value columns for tables\n");
+ fprintf(stderr, " -T\t\t# output in starbase/rdb table format\n");
+ fprintf(stderr, " -z\t\t# include regions with zero area in output\n");
+ fprintf(stderr, "\n(version: %s)\n", FUN_VERSION);
+ exit(1);
+}
+
+#ifdef ANSI_FUNC
+static double
+DConvert (char *buf, int type, int n)
+#else
+static double
+DConvert(buf, type, n)
+ char *buf;
+ int type;
+ int n;
+#endif
+{
+ int ival;
+ double dval=0.0;
+
+ switch(type){
+ case 'X':
+ if( n == 16 )
+ dval = (double)*(unsigned short *)buf;
+ else if( n == 32 )
+ dval = (double)*(unsigned int *)buf;
+ else
+ dval = (double)*(unsigned char *)buf;
+ break;
+ case 'B':
+ dval = (double)*(unsigned char *)buf;
+ break;
+ case 'I':
+ dval = (double)*(short *)buf;
+ break;
+ case 'U':
+ dval = (double)*(unsigned short *)buf;
+ break;
+ case 'J':
+ dval = (double)*(int *)buf;
+ break;
+ case 'K':
+#if HAVE_LONG_LONG == 0
+ gerror(stderr,
+ "64-bit data support not built (long long not available)\n");
+#endif
+ break;
+ case 'V':
+ dval = (double)*(unsigned int *)buf;
+ break;
+ case 'E':
+ dval = (double)*(float *)buf;
+ break;
+ case 'D':
+ dval = *(double *)buf;
+ break;
+ case 'L':
+ ival = (int)*(unsigned char *)buf;
+ if( !ival || (ival == 'F') || (ival == 'f') )
+ dval = 0.0;
+ else
+ dval = 1.0;
+ break;
+ default:
+ dval = 0.0;
+ break;
+ }
+ return(dval);
+}
+
+#ifdef ANSI_FUNC
+int
+main (int argc, char **argv)
+#else
+int
+main(argc, argv)
+ int argc;
+ char **argv;
+#endif
+{
+ int c;
+ int i, j, k, v;
+ int args;
+ int bkarea=0;
+ int bktype=BKG_VAL;
+ int dosum=0;
+ int dog=0;
+ int domatch=0;
+ int doradang=0;
+ int dobkgderr=1;
+ int doexp=0;
+ int dotim=0;
+ int dopixels=0;
+ int dozero=0;
+ int evimage=0;
+ int cold=' ';
+ int nintv=1;
+ int endian[NTYPE], type[NTYPE], rawsize[NTYPE], rowsize[NTYPE];
+ int x0[NTYPE], x1[NTYPE], y0[NTYPE], y1[NTYPE];
+ int dim1[4], dim2[4], block[4];
+ int *area[NTYPE];
+ int *savearea[NTYPE];
+ int nmask[NTYPE]={0,0};
+ int nreg[NTYPE]={0,0};
+ int valtypes[NTYPE]={0,0};
+ int valoffsets[NTYPE]={0,0};
+ int valns[NTYPE]={0,0};
+ int exptrans[NTYPE]={EXP_RATIO,EXP_RATIO};
+ double bkexp=0.0;
+ double dpp[NTYPE]={-1.0,-1.0};
+ double timecorr[NTYPE]={1.0,1.0};
+ double *cnts[MAXINTV][NTYPE];
+ double *savecnts[MAXINTV][NTYPE];
+ double *bncnts[MAXINTV];
+ double *bnerr[MAXINTV];
+ double *bscnts[MAXINTV];
+ double *bserr[MAXINTV];
+ double *exp[NTYPE];
+ double *saveexp[NTYPE];
+ double dppnorm=1.0;
+ double bkval=0.0;
+ char mode[SZ_LINE];
+ char tbuf[SZ_LINE];
+ char tbuf2[SZ_LINE];
+ char *intvs[MAXINTV];
+ char *s, *t, *u;
+ char *ebuf;
+ char *eptr;
+ char *fmt=NULL;
+ char *expstr=NULL;
+ char *intvstr=NULL;
+ char *valstr=NULL;
+ char *timestr=NULL;
+ char *radang=NULL;
+ char *cradang=NULL;
+ char *tradang=NULL;
+ char *region[NTYPE]={NULL,NULL};
+ char *filtstr[NTYPE]={NULL,NULL};
+ char *bincols[NTYPE]={NULL,NULL};
+ char *valname[2]={NULL,NULL};
+ char *name[4]={NULL,NULL,NULL, NULL};
+ Fun fun[4]={NULL,NULL,NULL,NULL};
+ FITSHead header[NTYPE]={NULL,NULL};
+ FilterMask masks[NTYPE]={NULL,NULL};
+ Filter filter[NTYPE]={NULL,NULL};
+ struct WorldCoor *wcs[4]={NULL,NULL,NULL,NULL};
+
+ /* clear array elelemts to mark them as unallocated */
+ memset(cnts, 0, MAXINTV*NTYPE*sizeof(double *));
+ memset(savecnts, 0, MAXINTV*NTYPE*sizeof(double *));
+ memset(bncnts, 0, NTYPE*sizeof(double *));
+ memset(bnerr, 0, NTYPE*sizeof(double *));
+ memset(bscnts, 0, NTYPE*sizeof(double *));
+ memset(bserr, 0, NTYPE*sizeof(double *));
+ memset(exp, 0, NTYPE*sizeof(double *));
+ memset(saveexp, 0, NTYPE*sizeof(double *));
+ memset(area, 0, NTYPE*sizeof(int *));
+ memset(savearea, 0, NTYPE*sizeof(int *));
+ memset(intvs, 0, MAXINTV*sizeof(char *));
+
+ /* exit on gio errors */
+ if( !getenv("GERROR") )
+ setgerror(2);
+
+ /* we want the args in the same order in which they arrived, and
+ gnu getopt sometimes changes things without this */
+ putenv("POSIXLY_CORRECT=true");
+
+ /* process switch arguments */
+ while ((c = getopt(argc, argv, "e:gGhi:mprst:v:w:zEIT")) != -1){
+ switch(c){
+ case 'e':
+ expstr = optarg;
+ break;
+ case 'g':
+ dog = 1;
+ break;
+ case 'G':
+ dog = 2;
+ break;
+ case 'h':
+ usage(argv[0]);
+ break;
+ case 'i':
+ intvstr = optarg;
+ break;
+ case 'm':
+ domatch = 1;
+ break;
+ case 'p':
+ dopixels = 1;
+ break;
+ case 'r':
+ doradang = 1;
+ break;
+ case 's':
+ dosum = 1;
+ break;
+ case 't':
+ timestr = optarg;
+ break;
+ case 'v':
+ valstr = optarg;
+ break;
+ case 'w':
+ putenv("FUNCNTS_EXPTRANS=wcs");
+ expstr = optarg;
+ break;
+ case 'z':
+ dozero = 1;
+ break;
+ case 'E':
+ evimage = 0;
+ break;
+ case 'I':
+ evimage = 1;
+ break;
+ case 'T':
+ cold = '\t';
+ break;
+ }
+ }
+
+ /* get maxrow,if user-specified */
+ if( (s=getenv("FUN_MAXROW")) != NULL )
+ maxrow = atoi(s);
+
+ /* check for required arguments */
+ args = argc - optind;
+ if( args < 1 ) usage(argv[0]);
+
+ /* arg 1: source file name */
+ name[SRC] = xstrdup(argv[optind+0]);
+
+ /* arg 2: source region */
+ if( (args == 1) || !*(region[SRC] = argv[optind+1]) )
+ region[SRC] = "field()";
+
+ /* arg 3: background region */
+ if( args >= 3 ){
+ if( args == 3 ){
+ name[BKG] = NULL;
+ region[BKG] = argv[optind+2];
+ }
+ else{
+ name[BKG] = xstrdup(argv[optind+2]);
+ region[BKG] = argv[optind+3];
+ }
+ /* check for constant numeric value -- background counts */
+ bkval = strtod(region[BKG], &s);
+ /* if we did not get a valid numeric constant, its a region */
+ if( region[BKG] && *region[BKG] && (s == region[BKG]) ){
+ bktype = BKG_ALL;
+ bkval = 0.0;
+ }
+ else{
+ /* can't have background file and background value */
+ if( args == 4 ) usage(argv[0]);
+ bktype = BKG_VAL;
+ region[BKG] = NULL;
+ }
+ }
+
+ /* get value column names, if necessary */
+ if( valstr ){
+ int ip=0;
+ int i=0;
+ newdtable(";");
+ for(i=0; i<NTYPE; i++){
+ if( word(valstr, tbuf, &ip) && *tbuf ){
+ if( strcasecmp(tbuf, "$none") )
+ valname[i] = xstrdup(tbuf);
+ else
+ valname[i] = xstrdup("-");
+ }
+ }
+ freedtable();
+ }
+
+ /* open the source FITS file */
+ if( !(fun[SRC] = FunOpen(name[SRC], "r", NULL)) )
+ gerror(stderr, "can't FunOpen source file (or find extension): %s\n",
+ name[SRC]);
+ /* get required information from funtools structure */
+ FunInfoGet(fun[SRC], FUN_ENDIAN, &endian[SRC],
+ FUN_TYPE, &type[SRC], FUN_HEADER, &header[SRC],
+ FUN_SECT_X0, &x0[SRC], FUN_SECT_X1, &x1[SRC],
+ FUN_SECT_Y0, &y0[SRC], FUN_SECT_Y1, &y1[SRC],
+ FUN_SECT_DIM1, &dim1[SRC], FUN_SECT_DIM2, &dim2[SRC],
+ FUN_SECT_BLOCK, &block[SRC], FUN_WCS, &wcs[SRC],
+ FUN_RAWSIZE, &rawsize[SRC], FUN_BINCOLS, &bincols[SRC],
+ FUN_ROWSIZE, &rowsize[SRC],
+ 0);
+ /* set dim2 to 1 for a 1D image */
+ if( dim2[SRC] == 0 ) dim2[SRC] = 1;
+ /* make up filter mode string using source bincols */
+ strcpy(mode, "type=image");
+ /* add the binning key */
+ if( bincols[SRC] ){
+ strcat(mode, ",");
+ strcat(mode, bincols[SRC]);
+ }
+
+ /* open the source region */
+ filter[SRC] = FilterOpen(header[SRC], region[SRC], mode);
+ if( filter[SRC] && (filter[SRC] != NOFILTER) ){
+ /* retrieve region mask segments */
+ nmask[SRC] = FilterImage(filter[SRC],
+ x0[SRC], x1[SRC], y0[SRC], y1[SRC],
+ block[SRC], &masks[SRC], &nreg[SRC]);
+
+ /* filtstr[SRC] = xstrdup(FilterString(filter[SRC])); */
+ /* for display, we have to add # comment chars after each \n */
+ t = FilterString(filter[SRC]);
+ filtstr[SRC] = calloc(strlen(t)*3, sizeof(char));
+ for(u=filtstr[SRC]; *t; t++){
+ *u++ = *t;
+ if( *t == '\n' ){
+ *u++ = '#';
+ *u++ = ' ';
+ }
+ }
+ filtstr[SRC] = realloc(filtstr[SRC], strlen(filtstr[SRC])+1);
+ }
+ if( filter[SRC] ) FilterClose(filter[SRC]);
+
+ /* make sure we have something to do */
+ if( nreg[SRC] == 0 )
+ gerror(stderr, "no valid regions included in '%s'\n", region[SRC]);
+
+ /* gather up intervals, if necessary (must be done before memory alloc) */
+ if( intvstr ){
+ int ip=0;
+ int ntok=0;
+ newdtable(";");
+ nintv = 0;
+ *tbuf2 = '\0';
+ while( word(intvstr, tbuf, &ip) && *tbuf ){
+ /* check for default column as first token */
+ if( !ntok && !strpbrk(tbuf, "=><&|") ){
+ strcpy(tbuf2, tbuf);
+ strcat(tbuf2, "=");
+ }
+ else{
+ if( nintv >= MAXINTV ){
+ gwarning(stderr, "Too many intervals; ignoring: %s\n", tbuf);
+ continue;
+ }
+ intvs[nintv] = (char *)calloc(SZ_LINE, sizeof(char));
+ if( *tbuf2 ) strcat(intvs[nintv], tbuf2);
+ strcat(intvs[nintv], tbuf);
+ nintv++;
+ }
+ ntok++;
+ }
+ freedtable();
+ }
+
+ /* if a value name was specified, make sure its a calid column */
+ if( valname[SRC] && (*valname[SRC] != '-') ){
+ if( !FunColumnLookup(fun[SRC], valname[SRC], 0, NULL,
+ &valtypes[SRC], NULL,
+ &valoffsets[SRC], NULL, NULL) ){
+ gerror(stderr, "value column does not exist: %s\n", valname[SRC]);
+ }
+ }
+
+ /* allocate space for results */
+ for(i=0; i<nintv; i++){
+ cnts[i][SRC] = (double *)xcalloc(nreg[SRC]+1, sizeof(double));
+ savecnts[i][SRC] = (double *)xcalloc(nreg[SRC]+1, sizeof(double));
+ bncnts[i] = (double *)xcalloc(nreg[SRC]+1, sizeof(double));
+ bnerr[i] = (double *)xcalloc(nreg[SRC]+1, sizeof(double));
+ bscnts[i] = (double *)xcalloc(nreg[SRC]+1, sizeof(double));
+ bserr[i] = (double *)xcalloc(nreg[SRC]+1, sizeof(double));
+ }
+ area[SRC] = (int *)xcalloc(nreg[SRC]+1, sizeof(int));
+ exp[SRC] = (double *)xcalloc(nreg[SRC]+1, sizeof(double));
+ savearea[SRC] = (int *)xcalloc(nreg[SRC]+1, sizeof(int));
+ saveexp[SRC] = (double *)xcalloc(nreg[SRC]+1, sizeof(double));
+ /* get radii and angle string, if needed */
+ if( doradang ){
+ if( (radang = FilterRadAng()) )
+ radang = xstrdup(radang);
+ }
+
+ /* process background region, if necessary */
+ if( region[BKG] ){
+ /* open the background FITS file */
+ if( name[BKG] ){
+ if( !(fun[BKG] = FunOpen(name[BKG], "r", NULL)) )
+ gerror(stderr,
+ "can't FunOpen background file (or find extension): %s\n",
+ name[BKG]);
+ /* get required information from funtools structure */
+ FunInfoGet(fun[BKG], FUN_ENDIAN, &endian[BKG],
+ FUN_TYPE, &type[BKG], FUN_HEADER, &header[BKG],
+ FUN_SECT_X0, &x0[BKG], FUN_SECT_X1, &x1[BKG],
+ FUN_SECT_Y0, &y0[BKG], FUN_SECT_Y1, &y1[BKG],
+ FUN_SECT_DIM1, &dim1[BKG], FUN_SECT_DIM2, &dim2[BKG],
+ FUN_SECT_BLOCK, &block[BKG], FUN_WCS, &wcs[BKG],
+ FUN_RAWSIZE, &rawsize[BKG], FUN_BINCOLS, &bincols[BKG],
+ FUN_ROWSIZE, &rowsize[BKG],
+ 0);
+ /* set dim2 to 1 for a 1D image */
+ if( dim2[BKG] == 0 ) dim2[BKG] = 1;
+ /* make up filter mode string using back bincols */
+ strcpy(mode, "type=image");
+ /* add the binning key */
+ if( bincols[BKG] ){
+ strcat(mode, ",");
+ strcat(mode, bincols[BKG]);
+ }
+ }
+ /* use values from the source */
+ else{
+ fun[BKG] = fun[SRC];
+ endian[BKG] = endian[SRC];
+ type[BKG] = type[SRC]; header[BKG] = header[SRC];
+ x0[BKG] = x0[SRC]; x1[BKG] = x1[SRC];
+ y0[BKG] = y0[SRC]; y1[BKG] = y1[SRC];
+ dim1[BKG] = dim1[SRC]; dim2[BKG] = dim2[SRC];
+ block[BKG] = block[SRC];
+ rawsize[BKG] = rawsize[SRC];
+ rowsize[BKG] = rowsize[SRC];
+ bincols[BKG] = bincols[SRC];
+ }
+
+ /* open the background region */
+ filter[BKG] = FilterOpen(header[BKG], region[BKG], mode);
+ if( filter[BKG] && (filter[BKG] != NOFILTER) ){
+ /* retrieve region mask segments */
+ nmask[BKG] = FilterImage(filter[BKG],
+ x0[BKG], x1[BKG], y0[BKG], y1[BKG],
+ block[BKG], &masks[BKG], &nreg[BKG]);
+ /* see if we want to, and can, match source and bkgd regions */
+ if( domatch && (nreg[SRC] == nreg[BKG]) )
+ bktype = BKG_EACH;
+ /* filtstr[BKG] = xstrdup(FilterString(filter[BKG])); */
+ /* for display, we have to add # comment chars after each \n */
+ t = FilterString(filter[BKG]);
+ filtstr[BKG] = calloc(strlen(t)*3, sizeof(char));
+ for(u=filtstr[BKG]; *t; t++){
+ *u++ = *t;
+ if( *t == '\n' ){
+ *u++ = '#';
+ *u++ = ' ';
+ }
+ }
+ filtstr[BKG] = realloc(filtstr[BKG], strlen(filtstr[BKG])+1);
+ }
+ if( filter[BKG] ) FilterClose(filter[BKG]);
+
+ /* if a value name was specified, make sure its a calid column */
+ if( valname[BKG] ){
+ if( *valname[BKG] != '-' ){
+ if( !FunColumnLookup(fun[BKG], valname[BKG], 0, NULL,
+ &valtypes[BKG], NULL,
+ &valoffsets[BKG], NULL, NULL) ){
+ gerror(stderr, "value column does not exist: %s\n", valname[BKG]);
+ }
+ }
+ }
+ else{
+ if( valname[SRC] && (*valname[SRC] != '-') ){
+ valname[BKG] = xstrdup(valname[SRC]);
+ valtypes[BKG] = valtypes[SRC];
+ valoffsets[BKG] = valoffsets[SRC];
+ }
+ }
+
+ /* allocate space for results */
+ for(i=0; i<nintv; i++){
+ cnts[i][BKG] = (double *)xcalloc(nreg[BKG]+1, sizeof(double));
+ savecnts[i][BKG] = (double *)xcalloc(nreg[BKG]+1, sizeof(double));
+ }
+ area[BKG] = (int *)xcalloc(nreg[BKG]+1, sizeof(int));
+ exp[BKG] = (double *)xcalloc(nreg[BKG]+1, sizeof(double));
+ savearea[BKG] = (int *)xcalloc(nreg[BKG]+1, sizeof(int));
+ saveexp[BKG] = (double *)xcalloc(nreg[BKG]+1, sizeof(double));
+ }
+
+ /* look for degrees/pixel in source and background files -- we will
+ use these to normalize background area, if both are present */
+ if( name[BKG] ){
+ for(i=0; i<2; i++){
+ if( wcs[i] && iswcs(wcs[i]) ){
+ if( !wcs[i]->coorflip )
+ dpp[i] = ABS(wcs[i]->cdelt[0]) * block[i];
+ else
+ dpp[i] = ABS(wcs[i]->cdelt[1]) * block[i];
+ }
+ }
+ /* note that BOTH must be present or we do not do this normalization */
+ if( (dpp[SRC] > 0.0) && (dpp[BKG] > 0.0) )
+ dppnorm = (dpp[SRC]/dpp[BKG]) * (dpp[SRC]/dpp[BKG]);
+ }
+ /* get degrees/pixel for source; not used in norm, but users want to know */
+ else{
+ if( wcs[SRC] && iswcs(wcs[SRC]) ){
+ if( !wcs[SRC]->coorflip )
+ dpp[SRC] = ABS(wcs[SRC]->cdelt[0]) * block[SRC];
+ else
+ dpp[SRC] = ABS(wcs[SRC]->cdelt[1]) * block[SRC];
+ }
+ }
+
+ /* open exposure file(s), if necessary */
+ if( expstr ){
+ int ip=0;
+ int i=0;
+ int got=0;
+ newdtable(";");
+ for(i=0; i<NTYPE; i++){
+ if( word(expstr, tbuf, &ip) && *tbuf ){
+ name[NTYPE+i] = xstrdup(tbuf);
+ if( !(fun[NTYPE+i] = FunOpen(name[NTYPE+i], "r", NULL)) )
+ gerror(stderr, "can't FunOpen exp file: %s\n", name[NTYPE+i]);
+ FunInfoGet(fun[NTYPE+i],
+ FUN_SECT_DIM1, &dim1[NTYPE+i],
+ FUN_SECT_DIM2, &dim2[NTYPE+i],
+ FUN_SECT_BLOCK, &block[NTYPE+i], FUN_WCS, &wcs[NTYPE+i], 0);
+ /* set dim2 to 1 for a 1D image */
+ if( dim2[NTYPE+i] == 0 ) dim2[NTYPE+i] = 1;
+ doexp |= (i+1);
+ got++;
+ }
+ }
+ freedtable();
+ /* if we have wcs for the data and exposure, and the user asks for it,
+ use wcs for conversion */
+ if( (s=getenv("FUNCNTS_EXPTRANS")) && !strcasecmp(s, "wcs") ){
+ for(i=0; i<got; i++){
+ if( wcs[i] && wcs[NTYPE+i] )
+ exptrans[i] = EXP_WCS;
+ else
+ gerror(stderr,
+ "no WCS present for WCS-based exposure", name[NTYPE+i]);
+ }
+ }
+ }
+
+ /* if we have source exposure file but no background exposure file,
+ and if background data comes from source file , then use the source
+ exposure file for bkgd exposure as well */
+ if( (doexp == SCORR) && !name[BKG] && (bktype != BKG_VAL) ){
+ doexp |= BCORR;
+ fun[3] = fun[2];
+ dim1[3] = dim1[2];
+ dim2[3] = dim2[2];
+ block[3] = block[2];
+ }
+
+ /* get time correction values */
+ if( timestr ){
+ int ip=0;
+ int i=0;
+ int got;
+ double dval;
+ char *t;
+ newdtable(";");
+ for(i=0; i<2; i++){
+ if( word(timestr, tbuf, &ip) && *tbuf ){
+ dval = strtod(tbuf, &t);
+ if( t != tbuf ){
+ timecorr[i] = dval;
+ }
+ else{
+ timecorr[i] = FunParamGetd(fun[i], tbuf, 0, 1.0, &got);
+ if( !got ){
+ /* try pure upper case */
+ cluc(tbuf);
+ timecorr[i] = FunParamGetd(fun[i], tbuf, 0, 1.0, &got);
+ if( !got ){
+ gerror(stderr, "can't find time correction parameter: %s\n",
+ tbuf);
+ }
+ }
+ }
+ dotim |= (i+1);
+ }
+ }
+ freedtable();
+ }
+ /* if we have source time correction but no background time correction,
+ and if background data comes from source file , then use the source
+ correction for bkgd as well */
+ if( (dotim == SCORR) && !name[BKG] && (bktype != BKG_VAL) ){
+ dotim |= BCORR;
+ timecorr[BKG] = timecorr[SRC];
+ }
+
+ /* process separate source/background files */
+ if( name[BKG] ){
+ for(k=0; k<NTYPE; k++){
+ /* process counts in all regions */
+ switch(type[k]){
+ case FUN_IMAGE:
+ case FUN_ARRAY:
+ {
+ int y, lasty;
+ double *dbuf=NULL;
+ /* allocate a row buffer */
+ dbuf = xmalloc(dim1[k] * sizeof(double));
+ /* seed with impossible value so we load first line */
+ lasty = -1;
+ for(i=0; i<nmask[k]; i++){
+ y = masks[k][i].y;
+ if( y != lasty ){
+ if( !FunImageRowGet(fun[k], dbuf, y, y, "bitpix=-64") )
+ gerror(stderr, "can't FunImageRowGet: %d %s\n", y, name[k]);
+ lasty = y;
+ }
+ area[k][masks[k][i].region-1] +=
+ masks[k][i].xstop - masks[k][i].xstart + 1;
+ for(j=masks[k][i].xstart-1; j<=masks[k][i].xstop-1; j++){
+ if( !isnand(dbuf[j]) ){
+ cnts[0][k][masks[k][i].region-1] += dbuf[j];
+ }
+ }
+ }
+ /* free temp buffers */
+ if( dbuf ) xfree(dbuf);
+ }
+ break;
+ case FUN_TABLE:
+ case FUN_EVENTS:
+ if( evimage ){
+ int *iptr;
+ int *ibuf=NULL;
+ /* read data as ints */
+ if( !(ibuf = FunImageGet(fun[k], NULL, "bitpix=32")) )
+ gerror(stderr, "can't FunImageGet: %s\n", name[k]);
+ /* get source counts */
+ for(i=0; i<nmask[k]; i++){
+ iptr = &(ibuf[(masks[k][i].y-1)*dim1[k]]);
+ for(j=masks[k][i].xstart-1; j<=masks[k][i].xstop-1; j++){
+ cnts[0][k][masks[k][i].region-1] += iptr[j];
+ }
+ }
+ /* free up space */
+ if( ibuf ) xfree(ibuf);
+ }
+ else{
+ int *rbuf;
+ int got;
+ char *rawbuf;
+ Filter efilter[MAXINTV];
+ /* make sure we will have a valid image section in which to filter */
+ if( (x1[k]-x0[k] <= 0) || (y1[k]-y0[k] <= 0) )
+ gerror(stderr,
+ "invalid or zero image dimensions(s) for table (%s)\n",
+ bincols[k] ? bincols[k] : "invalid bincols?");
+ /* make up the mode string */
+ snprintf(mode, SZ_LINE,
+ "type=events,convert=%s,evsect=\"%d %d %d %d %d\"",
+ (endian[k] == is_bigendian()) ? "false" : "true",
+ x0[k], x1[k], y0[k], y1[k], block[k]);
+ /* add columns */
+ if( bincols[k] ){
+ strcat(mode, ",");
+ strcat(mode, bincols[k]);
+ }
+ /* open filters for all intervals */
+ for(v=0; v<nintv; v++){
+ strcpy(tbuf, region[k]);
+ if( intvs[v] ){
+ strcat(tbuf, "&&(");
+ strcat(tbuf, intvs[v]);
+ strcat(tbuf, ")");
+ }
+ efilter[v] = FilterOpen(header[k], tbuf, mode);
+ }
+ /* allocate region value buffer */
+ rbuf = xmalloc(maxrow*sizeof(int));
+ /* extract events */
+ while( (ebuf = FunTableRowGet(fun[k], NULL, maxrow, NULL, &got)) ){
+ /* get pointer to raw buffer */
+ FunInfoGet(fun[k], FUN_RAWBUF, &rawbuf, 0);
+ /* process all intervals in one pass through this data */
+ for(v=0; v<nintv; v++){
+ /* get events which pass the region filter */
+ if( efilter[v] && (efilter[v] != NOFILTER) &&
+ FilterEvents(efilter[v], rawbuf, rawsize[k], got, rbuf)){
+ /* loop through events, process those which are in a region */
+ for(i=0; i<got; i++){
+ eptr = (char *)(ebuf+(rowsize[k]*i)+valoffsets[k]);
+ if( rbuf[i] > 0 ){
+ if( valname[k] && (*valname[k] != '-') )
+ cnts[v][k][rbuf[i]-1] +=
+ DConvert(eptr, valtypes[k], valns[k]);
+ else
+ cnts[v][k][rbuf[i]-1] += 1;
+ }
+ }
+ }
+ }
+ /* free for next read */
+ if( ebuf ) xfree(ebuf);
+ }
+ /* Done with region values */
+ if( rbuf ) xfree(rbuf);
+ /* close filters */
+ for(v=0; v<nintv; v++){
+ if( efilter[v] && (efilter[v] != NOFILTER) ){
+ FilterClose(efilter[v]);
+ }
+ }
+ }
+ /* get area */
+ for(i=0; i<nmask[k]; i++){
+ area[k][masks[k][i].region-1] +=
+ masks[k][i].xstop - masks[k][i].xstart + 1;
+ }
+ break;
+ }
+ }
+ }
+ /* same file for source and background:
+ this code is optimized so that we only traverse the file once
+ */
+ else{
+ /* process counts in all regions */
+ switch(type[SRC]){
+ case FUN_IMAGE:
+ case FUN_ARRAY:
+ {
+ int state;
+ int y, lasty;
+ double *dbuf=NULL;
+ /* allocate a row buffer */
+ dbuf = xmalloc(dim1[SRC] * sizeof(double));
+ /* seed with impossible value so we load first line */
+ lasty = -1;
+ for(i=0, j=0; i<nmask[SRC] || j<nmask[BKG]; ){
+ if( (i<nmask[SRC]) && (j<nmask[BKG]) ) state = 3;
+ else if( i<nmask[SRC] ) state = 1;
+ else if( j<nmask[BKG] ) state = 2;
+ else break;
+ switch(state){
+ case 1:
+ sline:
+ y = masks[SRC][i].y;
+ if( y != lasty ){
+ if( !FunImageRowGet(fun[SRC], dbuf, y, y, "bitpix=-64") )
+ gerror(stderr, "can't FunImageRowGet: %d %s\n", y, name[SRC]);
+ lasty = y;
+ }
+ area[SRC][masks[SRC][i].region-1] +=
+ masks[SRC][i].xstop - masks[SRC][i].xstart + 1;
+ for(k=masks[SRC][i].xstart-1; k<=masks[SRC][i].xstop-1; k++){
+ if( !isnand(dbuf[k]) ){
+ cnts[0][SRC][masks[SRC][i].region-1] += dbuf[k];
+ }
+ }
+ i++;
+ break;
+ case 2:
+ bline:
+ y = masks[BKG][j].y;
+ if( y != lasty ){
+ if( !FunImageRowGet(fun[SRC], dbuf, y, y, "bitpix=-64") )
+ gerror(stderr, "can't FunImageRowGet: %d %s\n", y, name[SRC]);
+ lasty = y;
+ }
+ area[BKG][masks[BKG][j].region-1] +=
+ masks[BKG][j].xstop - masks[BKG][j].xstart + 1;
+ for(k=masks[BKG][j].xstart-1; k<=masks[BKG][j].xstop-1; k++){
+ if( !isnand(dbuf[k]) ){
+ cnts[0][BKG][masks[BKG][j].region-1] += dbuf[k];
+ }
+ }
+ j++;
+ break;
+ case 3:
+ if( masks[SRC][i].y <= masks[BKG][j].y ){
+ goto sline;
+ }
+ goto bline;
+ }
+ }
+ if( dbuf ) xfree(dbuf);
+ }
+ break;
+ case FUN_TABLE:
+ case FUN_EVENTS:
+ if( evimage ){
+ int *iptr;
+ int *ibuf=NULL;
+ /* read data as ints */
+ if( !(ibuf = FunImageGet(fun[SRC], NULL, "bitpix=32")) )
+ gerror(stderr, "can't FunImageGet: %s\n", name[SRC]);
+ /* get source counts */
+ for(i=0; i<nmask[SRC]; i++){
+ iptr = &(ibuf[(masks[SRC][i].y-1)*dim1[SRC]]);
+ for(j=masks[SRC][i].xstart-1; j<=masks[SRC][i].xstop-1; j++){
+ cnts[0][SRC][masks[SRC][i].region-1] += iptr[j];
+ }
+ }
+ /* get background counts and area, if necessary */
+ if( bktype != BKG_VAL ){
+ for(i=0; i<nmask[BKG]; i++){
+ iptr = &(ibuf[(masks[BKG][i].y-1)*dim1[SRC]]);
+ for(j=masks[BKG][i].xstart-1; j<=masks[BKG][i].xstop-1; j++){
+ cnts[0][BKG][masks[BKG][i].region-1] += iptr[j];
+ }
+ }
+ }
+ /* free up space */
+ if( ibuf ) xfree(ibuf);
+ }
+ else{
+ int *rbuf;
+ int got;
+ char *rawbuf;
+ Filter efilter[MAXINTV][NTYPE];
+ /* make sure we will have a valid image section in which to filter */
+ if( (x1[SRC]-x0[SRC] <= 0) || (y1[SRC]-y0[SRC] <= 0) )
+ gerror(stderr,
+ "invalid or zero image dimensions(s) for table (%s)\n",
+ bincols[SRC] ? bincols[SRC] : "invalid bincols?");
+ /* open new filters to filter events through regions */
+ snprintf(mode, SZ_LINE,
+ "type=events,convert=%s,evsect=\"%d %d %d %d %d\"",
+ (endian[SRC] == is_bigendian()) ? "false" : "true",
+ x0[SRC], x1[SRC], y0[SRC], y1[SRC], block[SRC]);
+ /* add columns */
+ if( bincols[SRC] ){
+ strcat(mode, ",");
+ strcat(mode, bincols[SRC]);
+ }
+ /* open filters for all intervals */
+ for(v=0; v<nintv; v++){
+ for(i=0; i<NTYPE; i++){
+ if( region[i] ){
+ strcpy(tbuf, region[i]);
+ if( intvs[v] ){
+ strcat(tbuf, "&&(");
+ strcat(tbuf, intvs[v]);
+ strcat(tbuf, ")");
+ }
+ efilter[v][i] = FilterOpen(header[i], tbuf, mode);
+ }
+ else{
+ efilter[v][i] = NULL;
+ }
+ }
+ }
+ /* allocate region value buffer */
+ rbuf = xmalloc(maxrow*sizeof(int));
+ /* extract and filter events */
+ while( (ebuf = FunTableRowGet(fun[SRC], NULL, maxrow, NULL, &got)) ){
+ /* get pointer to raw buffer */
+ FunInfoGet(fun[SRC], FUN_RAWBUF, &rawbuf, 0);
+ /* process all intervals in one pass through this data */
+ for(v=0; v<nintv; v++){
+ /* get events which pass the region filter */
+ for(i=0; i<NTYPE; i++){
+ if( efilter[v][i] && (efilter[v][i] != NOFILTER) &&
+ FilterEvents(efilter[v][i], rawbuf, rawsize[i], got, rbuf)){
+ /* count events which are in a region */
+ for(j=0; j<got; j++){
+ eptr = (char *)(ebuf+(rowsize[i]*j)+valoffsets[i]);
+ if( rbuf[j] > 0 ){
+ if( valname[i] && (*valname[i] != '-') )
+ cnts[v][i][rbuf[j]-1] +=
+ DConvert(eptr, valtypes[i], valns[i]);
+ else
+ cnts[v][i][rbuf[j]-1] += 1;
+ }
+ }
+ }
+ }
+ }
+ /* free for next read */
+ if( ebuf ) xfree(ebuf);
+ }
+ /* done with region values */
+ if( rbuf ) xfree(rbuf);
+ /* close filters */
+ for(v=0; v<nintv; v++)
+ for(i=0; i<NTYPE; i++){
+ if( efilter[v][i] && (efilter[v][i] != NOFILTER) ){
+ FilterClose(efilter[v][i]);
+ }
+ }
+ }
+ /* get source area */
+ for(i=0; i<nmask[SRC]; i++){
+ area[SRC][masks[SRC][i].region-1] +=
+ masks[SRC][i].xstop - masks[SRC][i].xstart + 1;
+ }
+ /* get background area, if necessary */
+ if( bktype != BKG_VAL ){
+ for(i=0; i<nmask[BKG]; i++){
+ area[BKG][masks[BKG][i].region-1] +=
+ masks[BKG][i].xstop - masks[BKG][i].xstart + 1;
+ }
+ }
+ break;
+ }
+ }
+
+ /* accumulate averge exposure for each region */
+ if( doexp ){
+ for(k=0; k<2; k++){
+ if( fun[NTYPE+k] ){
+ int ex, ey, lastey;
+ double d1, d2;
+ double *exbuf=NULL;
+ char *tname="unknown";
+ exbuf = xmalloc(dim1[NTYPE+k] * sizeof(double));
+ if( name[NTYPE+k] )
+ tname = name[NTYPE+k];
+ else if( (k>0) && name[NTYPE+k-1] )
+ tname = name[NTYPE+k-1];
+ switch(exptrans[k]){
+ case EXP_RATIO:
+ lastey = -1;
+ d1 = ((double)dim1[NTYPE+k]/(double)dim1[k]);
+ d2 = ((double)dim2[NTYPE+k]/(double)dim2[k]);
+ for(i=0; i<nmask[k]; i++){
+ ey = ((masks[k][i].y - 1) * d2) + 1;
+ if( ey != lastey ){
+ if( !FunImageRowGet(fun[NTYPE+k], exbuf, ey, ey, "bitpix=-64") )
+ gerror(stderr, "can't FunImageRowGet (exp file): %d %s\n",
+ ey, tname);
+ lastey = ey;
+ }
+ for(j=masks[k][i].xstart; j<=masks[k][i].xstop; j++){
+ ex = (j-1) * d1 + 1;
+ if( ex < 1 ) ex = 1;
+ if( ex > dim1[NTYPE+k] ) ex = dim1[NTYPE+k];
+ exp[k][masks[k][i].region-1] += exbuf[ex-1];
+ }
+ }
+ break;
+ case EXP_WCS:
+ lastey = -1;
+ for(i=0; i<nmask[k]; i++){
+ for(j=masks[k][i].xstart; j<=masks[k][i].xstop; j++){
+ double dval1, dval2;
+ double dex, dey;
+ int offscl;
+ /* convert data image pixels to ra/dec using wcs */
+ pix2wcs(wcs[k], (double)j,(double)masks[k][i].y, &dval1,&dval2);
+ /* convert ra/dec to exp image pixels using wcs */
+ wcs2pix(wcs[NTYPE+k], dval1, dval2, &dex, &dey, &offscl);
+ ex = (int)(dex+0.5); ey = (int)(dey+0.5);
+ if( ex < 1 ) ex = 1;
+ if( ex > dim1[NTYPE+k] ) ex = dim1[NTYPE+k];
+ if( ey < 1 ) ey = 1;
+ if( ey > dim2[NTYPE+k] ) ey = dim2[NTYPE+k];
+ if( ey != lastey ){
+ if( !FunImageRowGet(fun[NTYPE+k], exbuf, ey, ey, "bitpix=-64") )
+ gerror(stderr, "can't FunImageRowGet (exp file): %d %s\n",
+ ey, tname);
+ lastey = ey;
+ }
+ exp[k][masks[k][i].region-1] += exbuf[ex-1];
+ }
+ }
+ break;
+ default:
+ gerror(stderr, "unknown exposure conversion type");
+ }
+ if( exbuf ) xfree(exbuf);
+ }
+ /* since we already checked for source exposure, we must be missing
+ bkgd exposure. This is OK, unless we have a bkgd file. In that case,
+ we need to warn the user and make an assumption about errors */
+ else if( (k==1) && name[BKG] ){
+ gwarning(stderr,
+ "No exposure file was specified for the background file.\nWe therefore assume that the bkgd already has been corrected for\nexposure and that the error associated with each bkgd pixel is 0.\n");
+ dobkgderr = 0;
+ }
+ }
+ }
+
+ /* reset the interval index for next stage of proessing*/
+ v = 0;
+ /* come back here if we are processing multiple intervals */
+intvagain:
+ /* for data-based background, check validity of background area */
+ if( bktype != BKG_VAL ){
+ bkval = 0.0;
+ bkarea = 0;
+ bkexp = 0.0;
+ for(i=0; i<nreg[BKG]; i++){
+ bkval += cnts[v][BKG][i];
+ bkarea += area[BKG][i];
+ bkexp += exp[BKG][i];
+ }
+ /* if background area is 0, that is bad */
+ if( bkarea == 0 ){
+ gerror(stdout, "background has zero area\n");
+ }
+ }
+
+ /* display source header information */
+ fprintf(stdout, "# source\n");
+ fprintf(stdout, "# data_file:\t\t%s\n", name[SRC]);
+ if( intvs[v] )
+ fprintf(stdout, "# interval:\t\t%s\n", intvs[v]);
+ if( valname[SRC] && (*valname[SRC] != '-') )
+ fprintf(stdout, "# value column:\t%s\n", valname[SRC]);
+ if( dpp[SRC] > 0.0 )
+ fprintf(stdout, "# arcsec/pixel:\t%g\n", dpp[SRC]*ARCSEC_PER_DEG);
+ if( doexp & SCORR )
+ fprintf(stdout, "# exp_correction:\t%s\n", name[NTYPE+SRC]);
+ if( dotim & SCORR )
+ fprintf(stdout, "# time_correction:\t%g\n", timecorr[SRC]);
+
+ /* display bkgd header information */
+ fprintf(stdout, "# background\n");
+ if( name[BKG] ){
+ fprintf(stdout, "# data_file:\t\t%s\n", name[BKG]);
+ if( intvs[v] )
+ fprintf(stdout, "# interval:\t\t%s\n", intvs[v]);
+ if( valname[BKG] && (*valname[BKG] != '-') )
+ fprintf(stdout, "# value column:\t%s\n", valname[BKG]);
+ if( dpp[BKG] > 0.0 )
+ fprintf(stdout, "# arcsec/pixel:\t%g\n", dpp[BKG]*ARCSEC_PER_DEG);
+ if( doexp & BCORR )
+ fprintf(stdout, "# exp_correction:\t%s\n", name[NTYPE+BKG]);
+ if( dotim & BCORR )
+ fprintf(stdout, "# time_correction:\t%g\n", timecorr[BKG]);
+ if( dppnorm != 1.0 )
+ fprintf(stdout, "# wcs area norm factor:\t%g/%g (source/bkgd))\n",
+ dpp[SRC],dpp[BKG]);
+ }
+ else if( bktype != BKG_VAL ){
+ fprintf(stdout, "# data_file:\t\t%s\n", name[SRC]);
+ if( valname[SRC] && valname[BKG] && strcmp(valname[SRC], valname[BKG]) &&
+ (*valname[BKG] != '-') )
+ fprintf(stdout, "# value_column:\t%s\n", valname[BKG]);
+ }
+ else
+ fprintf(stdout, "# constant_value:\t%.6f\n", bkval);
+
+ /* dislay table information */
+ if( (dpp[SRC] > 0.0) && !dopixels )
+ s = "arcsec";
+ else
+ s = "pixel";
+ fprintf(stdout, "# column units\n");
+ fprintf(stdout, "# area:\t\t%s**2\n", s);
+ fprintf(stdout, "# surf_bri:\t\tcnts/%s**2%s%s\n",
+ s, (dotim & SCORR)? "/sec" : "", (doexp & SCORR)? "/expval" : "");
+ fprintf(stdout, "# surf_err:\t\tcnts/%s**2%s%s\n",
+ s, (dotim & SCORR)? "/sec" : "", (doexp & SCORR)? "/expval" : "");
+ if( doradang ){
+ fprintf(stdout, "# radii:\t\t%ss\n", s);
+ fprintf(stdout, "# angles:\t\tdegrees\n");
+ }
+
+ /* come back here if we also are outputting summed results */
+sumagain:
+ /* if we need to display sums, so the sum now, but save unsummed values,
+ because we will have to display those as well */
+ switch(dosum){
+ case 1:
+ memcpy(savecnts[v][SRC], cnts[v][SRC], (nreg[SRC]+1)*sizeof(double));
+ memcpy(savearea[SRC], area[SRC], (nreg[SRC]+1)*sizeof(int));
+ for(i=1; i<nreg[SRC]; i++){
+ cnts[v][SRC][i] += cnts[v][SRC][i-1];
+ area[SRC][i] += area[SRC][i-1];
+ }
+ if( bktype != BKG_VAL ){
+ memcpy(savecnts[v][BKG], cnts[v][BKG],
+ (nreg[BKG]+1)*sizeof(double));
+ memcpy(savearea[BKG], area[BKG], (nreg[BKG]+1)*sizeof(int));
+ for(i=1; i<nreg[BKG]; i++){
+ cnts[v][BKG][i] += cnts[v][BKG][i-1];
+ area[BKG][i] += area[BKG][i-1];
+ }
+ }
+ break;
+ case 2:
+ memcpy(cnts[v][SRC], savecnts[v][SRC], (nreg[SRC]+1)*sizeof(double));
+ memcpy(area[SRC], savearea[SRC], (nreg[SRC]+1)*sizeof(int));
+ if( bktype != BKG_VAL ){
+ memcpy(cnts[v][BKG], savecnts[v][BKG],
+ (nreg[BKG]+1)*sizeof(double));
+ memcpy(area[BKG], savearea[BKG], (nreg[BKG]+1)*sizeof(int));
+ }
+ break;
+ default:
+ break;
+ }
+
+ /* process the background */
+ switch(bktype){
+ /* constant background is counts/pixel */
+ case BKG_VAL:
+ for(i=0; i<nreg[SRC]; i++){
+ bncnts[v][i] = (bkval * area[SRC][i] * timecorr[SRC]);
+ bscnts[v][i] = cnts[v][SRC][i] - bncnts[v][i];
+ bserr[v][i] = sqrt(cnts[v][SRC][i]);
+ bnerr[v][i] = 0.0;
+ }
+ break;
+ case BKG_ALL:
+ switch( dosum ){
+ case 1:
+ bkval = cnts[v][BKG][nreg[BKG]-1];
+ bkarea = area[BKG][nreg[BKG]-1];
+ bkexp = exp[BKG][nreg[BKG]-1];
+ break;
+ case 2:
+ default:
+ bkval = 0.0;
+ bkarea = 0;
+ bkexp = 0.0;
+ /* get total background and background area */
+ for(i=0; i<nreg[BKG]; i++){
+ bkval += cnts[v][BKG][i];
+ bkarea += area[BKG][i];
+ bkexp += exp[BKG][i];
+ }
+ break;
+ }
+ {
+ double tempnorm=0.0;
+ /* subtract entire normalized background from each source region */
+ for(i=0; i<nreg[SRC]; i++){
+ /* get area and exposure normalization */
+ switch(doexp){
+ case NONE:
+ tempnorm = (double)area[SRC][i] / bkarea;
+ break;
+ case SCORR:
+ tempnorm = exp[SRC][i] / bkarea;
+ break;
+ case BCORR:
+ tempnorm = (double)area[SRC][i] / bkexp;
+ break;
+ case BOTH:
+ tempnorm = exp[SRC][i] / bkexp;
+ break;
+ }
+ /* add time normalization */
+ tempnorm *= (timecorr[SRC]/timecorr[BKG]);
+ /* add normalization due to different pixels sizes in src and bkg */
+ tempnorm *= dppnorm;
+ bncnts[v][i] = (bkval * tempnorm);
+ bscnts[v][i] = cnts[v][SRC][i] - bncnts[v][i];
+ bserr[v][i] = sqrt(cnts[v][SRC][i] + (tempnorm*tempnorm*bkval));
+ if( dobkgderr )
+ bnerr[v][i] = sqrt(bkval) * tempnorm;
+ else
+ bnerr[v][i] = 0.0;
+ }
+ }
+ break;
+ case BKG_EACH:
+ {
+ double tempnorm=0.0;
+ /* subtract each background from each source region */
+ for(i=0; i<nreg[SRC]; i++){
+ /* get area and exposure normalization */
+ switch(doexp){
+ case NONE:
+ tempnorm = (double)area[SRC][i] / (double)area[BKG][i];
+ break;
+ case SCORR:
+ tempnorm = exp[SRC][i] / (double)area[BKG][i];
+ break;
+ case BCORR:
+ tempnorm = (double)area[SRC][i] / (double)exp[BKG][i];
+ break;
+ case BOTH:
+ tempnorm = exp[SRC][i] / (double)exp[BKG][i];
+ break;
+ }
+ /* add time normalization */
+ tempnorm *= (timecorr[SRC]/timecorr[BKG]);
+ /* add normalization due to different pixels sizes in src and bkg */
+ tempnorm *= dppnorm;
+ bncnts[v][i] = (cnts[v][BKG][i] * tempnorm);
+ bscnts[v][i] = cnts[v][SRC][i] - bncnts[v][i];
+ bserr[v][i] = sqrt(cnts[v][SRC][i] +
+ (tempnorm * tempnorm * cnts[v][BKG][i]));
+ if( dobkgderr )
+ bnerr[v][i] = sqrt(cnts[v][BKG][i]) * tempnorm;
+ else
+ bnerr[v][i] = 0.0;
+ }
+ }
+ break;
+ }
+
+ /* display results */
+ /* display the main output table */
+ switch( dosum ){
+ case 1:
+ fprintf(stdout, "\n");
+ fprintf(stdout, "# summed background-subtracted results\n");
+ fprintf(stdout, "upto%c net_counts%c error", cold, cold);
+ break;
+ case 2:
+ default:
+ if( cold == '\t' ) fprintf(stdout, "\f");
+ fprintf(stdout, "\n");
+ fprintf(stdout, "# background-subtracted results\n");
+ fprintf(stdout, " reg%c net_counts%c error", cold, cold);
+ break;
+ }
+ fprintf(stdout, "%c background%c berror", cold, cold);
+ fprintf(stdout, "%c area%c surf_bri%c surf_err", cold, cold, cold);
+ if( doradang )
+ fprintf(stdout, "%c radius1%c radius2%c angle1%c angle2",
+ cold, cold, cold, cold);
+ fprintf(stdout, "\n");
+ fprintf(stdout, "----%c------------%c---------", cold, cold);
+ fprintf(stdout, "%c------------%c---------", cold, cold);
+ fprintf(stdout, "%c---------%c---------%c---------", cold, cold, cold);
+ if( doradang )
+ fprintf(stdout, "%c---------%c---------%c---------%c---------",
+ cold, cold, cold, cold);
+ fprintf(stdout, "\n");
+ if( radang )
+ newdtable(",");
+ cradang = radang;
+ for(i=0; i<nreg[SRC]; i++){
+ /* get next line from radii/angle string */
+ if( cradang ){
+ tradang = (char *)strchr(cradang, '\n');
+ if( tradang )
+ *tradang = '\0';
+ }
+ if( area[SRC][i] ){
+ double cntsperarea;
+ double errperarea;
+ double areasq;
+ if( doexp & SCORR ){
+ cntsperarea = bscnts[v][i]/(exp[SRC][i]*timecorr[SRC]);
+ errperarea = bserr[v][i]/(exp[SRC][i]*timecorr[SRC]);
+ areasq = area[SRC][i];
+ }
+ else{
+ cntsperarea = bscnts[v][i]/(area[SRC][i]*timecorr[SRC]);
+ errperarea = bserr[v][i]/(area[SRC][i]*timecorr[SRC]);
+ areasq = area[SRC][i];
+ }
+ /* if we know how to convert to cnts/pix**2 to cnts/arcsec**2, do it */
+ if( !dopixels && (dpp[SRC] > 0.0) ){
+ cntsperarea = (cntsperarea / (dpp[SRC]*dpp[SRC])) / ARCSEC_PER_DEGSQ;
+ errperarea = (errperarea / (dpp[SRC]*dpp[SRC])) / ARCSEC_PER_DEGSQ;
+ areasq = (areasq * (dpp[SRC]*dpp[SRC])) * ARCSEC_PER_DEGSQ;
+ }
+ /* get correctly precisioned format statement */
+ switch(dog){
+ case 0:
+ fmt = "%4d%c%12.3f%c%9.3f%c%12.3f%c%9.3f%c%9.2f%c%9.3f%c%9.3f";
+ break;
+ case 1:
+ fmt = "%4d%c%12.3g%c%9.3g%c%12.3g%c%9.3g%c%9.2g%c%9.3g%c%9.3g";
+ break;
+ case 2:
+ fmt = "%4d%c%.14g%c%.14g%c%.14g%c%.14g%c%.14g%c%.14g%c%.14g";
+ break;
+ }
+ fprintf(stdout, fmt,
+ i+1, cold,
+ bscnts[v][i], cold, bserr[v][i], cold,
+ bncnts[v][i], cold, bnerr[v][i], cold,
+ areasq, cold, cntsperarea, cold, errperarea);
+ /* display values from this line of radii/angles */
+ if( doradang && cradang ){
+ int ip=0;
+ double dval;
+ for(j=0; j<4; j++){
+ if( word(cradang, tbuf, &ip) && strcmp(tbuf, "NA") ){
+ dval = strtod(tbuf, NULL);
+ if( (j<2) && !dopixels && (dpp[SRC]>0.0) )
+ fprintf(stdout, "%c%9.3f", cold, (dval*dpp[SRC]*ARCSEC_PER_DEG));
+ else
+ fprintf(stdout, "%c%9.3f", cold, dval);
+ }
+ else{
+ fprintf(stdout, "%c%9.9s", cold, "NA");
+ }
+ }
+ }
+ /* new-line at end */
+ fprintf(stdout, "\n");
+ }
+ /* might have to display zero area pixels */
+ else if( dozero ){
+ /* get correctly precisioned format statement */
+ switch(dog){
+ case 0:
+ fmt = "%4d%c%12.3f%c%9.3f%c%12.3f%c%9.3f%c%9.2f%c%9.3f%c%9.3f";
+ break;
+ case 1:
+ fmt = "%4d%c%12.3g%c%9.3g%c%12.3g%c%9.3g%c%9.2g%c%9.3g%c%9.3g";
+ break;
+ case 2:
+ fmt = "%4d%c%.14g%c%.14g%c%.14g%c%.14g%c%.14g%c%.14g%c%.14g";
+ break;
+ }
+ fprintf(stdout, fmt, i+1, cold, 0.0, cold, 0.0, cold, 0.0, cold, 0.0,
+ cold, 0.0, cold, 0.0, cold, 0.0);
+ /* add the correct radii and angle info, to make plotting easier */
+ if( doradang && cradang ){
+ int ip=0;
+ double dval;
+ for(j=0; j<4; j++){
+ if( word(cradang, tbuf, &ip) && strcmp(tbuf, "NA") ){
+ dval = strtod(tbuf, NULL);
+ if( (j<2) && !dopixels && (dpp[SRC]>0.0) )
+ fprintf(stdout, "%c%9.3f", cold, (dval*dpp[SRC]*ARCSEC_PER_DEG));
+ else
+ fprintf(stdout, "%c%9.3f", cold, dval);
+ }
+ else{
+ fprintf(stdout, "%c%9.9s", cold, "NA");
+ }
+ }
+ }
+ /* new-line at end */
+ fprintf(stdout, "\n");
+ }
+ /* bump to next line of radii/angles */
+ if( tradang ){
+ cradang = tradang+1;
+ /* put back the cr in case we pass through again */
+ *tradang = '\n';
+ }
+ }
+ if( radang )
+ freedtable();
+ fprintf(stdout, "\n");
+ fflush(stdout);
+
+ /* if we just summed, go back and display unsummed values */
+ switch(dosum){
+ case 1:
+ dosum = 2;
+ goto sumagain;
+ /* set flag back to 1 in case we have multiple intervals to process */
+ case 2:
+ dosum = 1;
+ break;
+ default:
+ break;
+ }
+
+ /* display raw source counts */
+ if( dosum ){
+ int tarea=0;
+ double tcnts=0;
+ if( cold == '\t' ) fprintf(stdout, "\f");
+ fprintf(stdout, "\n");
+ /* display source info */
+ if( filtstr[SRC] ){
+ fprintf(stdout, "# source_region(s):\n");
+ fprintf(stdout, "# %s\n\n", filtstr[SRC]);
+ }
+ fprintf(stdout, "# summed_source_data\n");
+ fprintf(stdout,
+ " reg%c counts%c pixels%c sumcnts%c sumpix\n",
+ cold, cold, cold, cold);
+ fprintf(stdout,
+ "----%c------------%c---------%c------------%c---------\n",
+ cold, cold, cold, cold);
+ for(i=0; i<nreg[SRC]; i++){
+ tcnts += cnts[v][SRC][i];
+ tarea += area[SRC][i];
+ /* get correctly precisioned format statement */
+ switch(dog){
+ case 0:
+ fmt = "%4d%c%12.3f%c%9d%c%12.3f%c%9d\n";
+ break;
+ case 1:
+ fmt = "%4d%c%12.3g%c%9d%c%12.3g%c%9d\n";
+ break;
+ case 2:
+ fmt = "%4d%c%.14g%c%9d%c%.14g%c%9d\n";
+ break;
+ }
+ fprintf(stdout, fmt,
+ i+1, cold,
+ cnts[v][SRC][i], cold, area[SRC][i], cold,
+ tcnts, cold, tarea);
+ }
+ } else{
+ if( cold == '\t' ) fprintf(stdout, "\f");
+ fprintf(stdout, "\n");
+ /* display source info */
+ if( filtstr[SRC] ){
+ fprintf(stdout, "# source_region(s):\n");
+ fprintf(stdout, "# %s\n\n", filtstr[SRC]);
+ }
+ fprintf(stdout, "# source_data\n");
+ fprintf(stdout, " reg%c counts%c pixels", cold, cold);
+ if( doexp & SCORR )
+ fprintf(stdout, "%c avg_exp", cold);
+ fprintf(stdout, "\n");
+ fprintf(stdout, "----%c------------%c---------", cold, cold);
+ if( doexp & SCORR )
+ fprintf(stdout, "%c---------", cold);
+ fprintf(stdout, "\n");
+ for(i=0; i<nreg[SRC]; i++){
+ /* get correctly precisioned format statement */
+ switch(dog){
+ case 0:
+ fmt = "%4d%c%12.3f%c%9d";
+ break;
+ case 1:
+ fmt = "%4d%c%12.3g%c%9d";
+ break;
+ case 2:
+ fmt = "%4d%c%.14g%c%9d";
+ break;
+ }
+ fprintf(stdout, fmt, i+1, cold, cnts[v][SRC][i], cold, area[SRC][i]);
+
+ if( doexp & SCORR ){
+ /* get correctly precisioned format statement */
+ switch(dog){
+ case 0:
+ fmt = "%c%9.3f";
+ break;
+ case 1:
+ fmt = "%c%9.3g";
+ break;
+ case 2:
+ fmt = "%c%.14g";
+ break;
+ }
+ if( area[SRC][i] > 0 )
+ fprintf(stdout, fmt, cold, exp[SRC][i]/area[SRC][i]);
+ else
+ fprintf(stdout, fmt, cold, 0.0);
+ }
+ fprintf(stdout, "\n");
+ }
+ }
+ fprintf(stdout, "\n");
+ fflush(stdout);
+
+ /* display raw background info */
+ switch(bktype){
+ case BKG_VAL:
+ break;
+ case BKG_ALL:
+ if( cold == '\t' ) fprintf(stdout, "\f");
+ fprintf(stdout, "\n");
+ if( filtstr[BKG] ){
+ fprintf(stdout, "# background_region(s)\n");
+ fprintf(stdout, "# %s\n\n", filtstr[BKG]);
+ }
+ fprintf(stdout, "# background_data\n");
+ fprintf(stdout, " reg%c counts%c pixels", cold, cold);
+ if( doexp & BCORR )
+ fprintf(stdout, "%c avg_exp", cold);
+ fprintf(stdout, "\n");
+ fprintf(stdout, "----%c------------%c---------", cold, cold);
+ if( doexp & BCORR )
+ fprintf(stdout, "%c---------", cold);
+ fprintf(stdout, "\n");
+ /* get correctly precisioned format statement */
+ switch(dog){
+ case 0:
+ fmt = "%s%c%12.3f%c%9d";
+ break;
+ case 1:
+ fmt = "%s%c%12.3g%c%9d";
+ break;
+ case 2:
+ fmt = "%s%c%.14g%c%9d";
+ break;
+ }
+ fprintf(stdout, fmt, "all ", cold, bkval, cold, bkarea);
+ if( doexp & BCORR ){
+ /* get correctly precisioned format statement */
+ switch(dog){
+ case 0:
+ fmt = "%c%9.3f";
+ break;
+ case 1:
+ fmt = "%c%9.3g";
+ break;
+ case 2:
+ fmt = "%c%.14g";
+ break;
+ }
+ fprintf(stdout, fmt, cold, bkexp/bkarea);
+ }
+ fprintf(stdout, "\n");
+ break;
+ case BKG_EACH:
+ if( dosum ){
+ int tarea=0;
+ double tcnts=0;
+ if( cold == '\t' ) fprintf(stdout, "\f");
+ fprintf(stdout, "\n");
+ if( filtstr[BKG] ){
+ fprintf(stdout, "# background_region(s)\n");
+ fprintf(stdout, "# %s\n\n", filtstr[BKG]);
+ }
+ fprintf(stdout, "# summed_background_data\n");
+ fprintf(stdout,
+ " reg%c counts%c pixels%c sumcnts%c sumpix\n",
+ cold, cold, cold, cold);
+ fprintf(stdout,
+ "----%c------------%c---------%c------------%c---------\n",
+ cold, cold, cold, cold);
+ for(i=0; i<nreg[BKG]; i++){
+ tcnts += cnts[v][BKG][i];
+ tarea += area[BKG][i];
+ /* get correctly precisioned format statement */
+ switch(dog){
+ case 0:
+ fmt = "%4d%c%12.3f%c%9d%c%12.3f%c%9d\n";
+ break;
+ case 1:
+ fmt = "%4d%c%12.3g%c%9d%c%12.3g%c%9d\n";
+ break;
+ case 2:
+ fmt = "%4d%c%.14g%c%9d%c%.14g%c%9d\n";
+ break;
+ }
+ fprintf(stdout, fmt, i+1, cold,
+ cnts[v][BKG][i], cold, area[BKG][i], cold,
+ tcnts, cold, tarea);
+ }
+ } else{
+ if( cold == '\t' ) fprintf(stdout, "\f");
+ fprintf(stdout, "\n");
+ if( filtstr[BKG] ){
+ fprintf(stdout, "# background_region(s)\n");
+ fprintf(stdout, "# %s\n\n", filtstr[BKG]);
+ }
+ fprintf(stdout, "# background_data\n");
+ fprintf(stdout, " reg%c counts%c pixels", cold, cold);
+ if( doexp & BCORR )
+ fprintf(stdout, "%c avg_exp", cold);
+ fprintf(stdout, "\n");
+ fprintf(stdout, "----%c------------%c---------", cold, cold);
+ if( doexp & BCORR )
+ fprintf(stdout, "%c---------", cold);
+ fprintf(stdout, "\n");
+ for(i=0; i<nreg[BKG]; i++){
+ /* get correctly precisioned format statement */
+ switch(dog){
+ case 0:
+ fmt = "%4d%c%12.3f%c%9d";
+ break;
+ case 1:
+ fmt = "%4d%c%12.3g%c%9d";
+ break;
+ case 2:
+ fmt = "%4d%c%.14g%c%9d";
+ break;
+ }
+ fprintf(stdout, fmt, i+1, cold,
+ cnts[v][BKG][i], cold, area[BKG][i]);
+ if( doexp & BCORR ){
+ /* get correctly precisioned format statement */
+ switch(dog){
+ case 0:
+ fmt = "%c%9.3f";
+ break;
+ case 1:
+ fmt = "%c%9.3g";
+ break;
+ case 2:
+ fmt = "%c%.14g";
+ break;
+ }
+ if( area[BKG][i] > 0 )
+ fprintf(stdout, fmt, cold, exp[BKG][i]/area[BKG][i]);
+ else
+ fprintf(stdout, fmt, cold, 0.0);
+ }
+ fprintf(stdout, "\n");
+ }
+ }
+ break;
+ }
+ fprintf(stdout, "\n");
+ fflush(stdout);
+
+ /* process the next interval, if necessary */
+ if( ++v < nintv ){
+ fprintf(stdout, "\014\n");
+ goto intvagain;
+ }
+
+ /* cleanup */
+ for(i=0; i<nintv; i++){
+ if( bncnts[i] ) xfree(bncnts[i]);
+ if( bnerr[i] ) xfree(bnerr[i]);
+ if( bscnts[i] ) xfree(bscnts[i]);
+ if( bserr[i] ) xfree(bserr[i]);
+ for(j=0; j<2; j++){
+ if( cnts[i][j] ) xfree(cnts[i][j]);
+ if( savecnts[i][j] ) xfree(savecnts[i][j]);
+ }
+ if( intvs[i] ) xfree(intvs[i]);
+ }
+ for(i=0; i<2; i++){
+ if( masks[i] ) xfree(masks[i]);
+ if( area[i] ) xfree(area[i]);
+ if( exp[i] ) xfree(exp[i]);
+ if( savearea[i] ) xfree(savearea[i]);
+ if( saveexp[i] ) xfree(saveexp[i]);
+ if( filtstr[i] ) xfree(filtstr[i]);
+ if( valname[i] ) xfree(valname[i]);
+ }
+ if( radang ) xfree(radang);
+ for(i=0; i<4; i++){
+ if( name[i] ){
+ xfree(name[i]);
+ if( fun[i] ) FunClose(fun[i]);
+ }
+ }
+
+ return(0);
+}
generated by cgit v0.12 at 2024-08-28 23:51:56 (GMT)