/*
* 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);
}