rm funtools for update

1 files changed, 0 insertions, 1693 deletions

diff --git a/funtools/funcnts.c b/funtools/funcnts.c
deleted file mode 100644
index 6da09ed..0000000
--- a/funtools/funcnts.c
+++ /dev/null
@@ -1,1693 +0,0 @@
-/*
- *	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);
-}