diff options
Diffstat (limited to 'funtools/funcnts.c')
-rw-r--r-- | funtools/funcnts.c | 1693 |
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); +} |