/*****************************************************************************
******************************************************************************
**<AUTO>
**
** FILE:	tclMerge.cc
**
**<HTML>
**	This file contains TCL functions to merge, unmerge, and generate
**	merging statistics.
**</HTML>
**</AUTO>
**
**
** ENVIRONMENT:
**	C++.
**
** REQUIRED PRODUCTS:
**	photoSch
**
******************************************************************************
******************************************************************************
*/

#include <arObject.h>
#include <arObjectList.h>
#include <arFramesPipeRun.h>
#include <arRun.h>
#include <arProgram.h>
#include <arCCDProgram.h>
#include <arAstromCalib.h>
#include <arFinalPhotoCalib.h>
#include <arProduct.h>
#include <arParamValue.h>
#include "tsMerge.h"
#include <taCosmic.h>
#include <dervish.h>
extern "C" {
#include <slalib.h>
#include <cpgplot.h>
}
#include "arObjectStatus.h"
#include <taCalibObj.h>

#include <arSegment.h>
#include <arTargetType.h>
#include <arTarget.h>

const int nFilters = 5;             /* Number of filters */
const int nColors = 4;              /* Number of colors */

/*============================================================================
**<AUTO EXTRACT>
**
** TCL VERB: xSetObjectStatusDontCall
**
**<HTML>
**	Set the object status for all objects in a single Frames Pipeline Run,
**	and merge objects in the overlap regions between adjacent fields.
**	This consists of three steps, performed in the following order:
**	<OL>
**	<LI> Set each object's status (i.e., set/don't set the GOOD flag)
**	based on the values of the object flags.
**	<LI> Merge objects lying in the overlap regions between adjacent
**	fields.
**	<LI> Set/don't set the OK_RUN and DUPLICATE (for GOOD objects only),
**	determining whether this object is a primary detection for this run or
**	not.  This includes resolving the matched pairs of objects in the
**	overlap region, requiring that exactly one object in each matched pair
**	has its status set to OK_RUN.
**	</OL>
**	The following object flags are set by this routine:
**	<DL>
**	<DT> SET
**	<DD> All objects will have this flag set, indicating the objects have
**	been processed by "setObjectStatus".
**	<DT> GOOD
**	<DD> If set, then this object has been properly measured, as
**	determined by the object flags, to the extent
**	that it MAY be used in the primary survey, depending on the status of
**	its other flags.  If not set, then the object has NOT been properly
**	measured, and it may NOT be used in the primary survey.
**	</DL>
**	For each GOOD object, the following two flags MAY be set:
**	<DL>
**	<DT> OK_RUN
**	<DD> This is the primary detection of this object for this Frames
**	Pipeline Run, and MAY be used in the primary survey, depending on the
**	status of its other flags.
**	<DT> DUPLICATE
**	<DD> This object has a duplicate detection in an adjacent field
**	(thus both objects lay in the overlap region between the two fields).
**	A link will be made in the database to the duplicate detection.
**	Such pairs of objects are useful for quality
**	analysis internal to this run of the Frames Pipeline (i.e., the
**	measured parameters should be all but exactly equal for each object
**	in a pair).  For all such pairs, exactly one object in the pair
**	will also have the OK_RUN flag set.
**	</DL>
**	Objects match if their separation is less than
**	or equal to the specified match 'radius'.  If more than one object in
**	an adjacent field matches a given object, then only the closest match
**	is merged.  This returns a CHAIN of TSMATCHSTATs, with one set of
**	statistics for each field in the run.	The returned statistics give
**	the row and column statistics in pixels, and the uncalibrated PSF
**	magnitude statistics in magnitudes.  Only objects classified STAR,
**	GALAXY, or UNK are matched.
**<P>
**	This verb should NOT be called directly.  The verb "setObjectStatus"
**	should be used; it calls this verb.
**</HTML>
**
**</AUTO>
**============================================================================
*/
static char setObjectStatusDontCallCmd[] = "xSetObjectStatusDontCall";
static int setObjectStatusDontCallFlg = FTCL_ARGV_NO_LEFTOVERS;
static ftclArgvInfo setObjectStatusDontCallTbl[] = {
   {NULL, FTCL_ARGV_HELP, NULL, NULL,
    "Set object status and self merge a Frames Pipeline Run\n"},
   {"<run>", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of Frames Pipeline run"},
   {"<matchRadius>", FTCL_ARGV_DOUBLE, NULL, NULL,
    "Matching 'radius' (pixels)"},
   {"<nSigma>", FTCL_ARGV_DOUBLE, NULL, NULL,
    "Compile matching statistics for nSigma detections (PSF, in r') only"},
   {NULL, FTCL_ARGV_END, NULL, NULL, NULL}
};

static int
tclSetObjectStatusDontCall (ClientData clientData,
			    Tcl_Interp *interp,
			    int argc,
			    char **argv)
{
  ooHandle(arFramesPipeRun) *fp = NULL;
  char *framesName = NULL;
  double matchRadius, nSigma;
  int status;

  // Parse the arguments
  setObjectStatusDontCallTbl[1].dst = &framesName;
  setObjectStatusDontCallTbl[2].dst = &matchRadius;
  setObjectStatusDontCallTbl[3].dst = &nSigma;
  if ((status = shTclParseArgv(interp, &argc, argv, setObjectStatusDontCallTbl,
			       setObjectStatusDontCallFlg,
			       setObjectStatusDontCallCmd))
      != FTCL_ARGV_SUCCESS)
    return status;
   
  // Get handle to frames pipeline run
  if (shTclAddrGetFromName (interp, framesName, (void **)&fp,
			    "ooHandle(arFramesPipeRun)") != TCL_OK)
    {
      Tcl_SetResult(interp,
		    "setObjectStatusDontCall: error getting handle to frames pipeline run",
		    TCL_VOLATILE);
      return TCL_ERROR;
    }
   
  // Verify that this run has not previously been through setObjectStatus
  if ((*fp)->setObjectStatus().isValid() == oocTrue)
    {
      Tcl_SetResult(interp, 
		    "setObjectStatusDontCall: run was previously processed by setObjectStatus",
		    TCL_VOLATILE);
      return TCL_ERROR;
    }

  // Get number of rows per frame.
  int camRow = 1;  // Any photometric CCD row will do
  int rowsPerField =
    (*fp)->imagingRun()->program()->ccd(camRow,(*fp)->camCol())->nRows();

  // HARDWIRED: number of overlap rows per frame.
  int overlapRows = 128;

  // Fetch the range of fields
  uint16 field0 = (*fp)->field0();
  uint16 nFields = (*fp)->nFields();

  // Allocate a CHAIN to return the matching statistics for each field.
  CHAIN *ch = shChainNewBlock("TSMATCHSTATS", nFields);
  if (ch == NULL)
    {
      Tcl_SetResult(interp, 
		    "setObjectStatusDontCall: couldn't allocate stats chain",
		    TCL_VOLATILE);
      return TCL_ERROR;
    }

  // Allocate pointers for some dynamic structures
  TSMATCHSTATS *statsNew = NULL;
  TSMATCHSTATS *statsOld = NULL;
  ooTVArray(ooHandle(arObject)) *handlesNew = NULL;
  ooTVArray(ooHandle(arObject)) *handlesOld = NULL;
  TSOBJ *arrayNew = NULL;
  TSOBJ *arrayOld = NULL;

  // For each field ...
  CURSOR_T cursor = shChainCursorNew(ch);
  for (int32 field = field0; field < field0 + nFields; field++)
    {
      printf("%4d", field);
      fflush(stdout);
      // Save the pointers for the previous field
      if (field != field0)
	{
	  arrayOld = arrayNew;
	  handlesOld = handlesNew;
	  statsOld = statsNew;
	}

      // Get the stats structure for this field
      if ((statsNew = (TSMATCHSTATS*) shChainWalk(ch, cursor, NEXT)) == NULL)
	{
	  // Error.  Clean up and return an error.
	  shChainCursorDel(ch, cursor);
	  taGenericChainDestroy(ch);
	  if (handlesOld != NULL) delete handlesOld;
	  if (arrayOld != NULL) shFree(arrayOld);
	  char fieldString[10];
	  sprintf(fieldString, "%d", field);
	  Tcl_AppendResult(interp, setObjectStatusDontCallCmd,
			   ": error getting stats in field ", fieldString,
			   NULL);
	  return TCL_ERROR;
	}
      tsMatchStatsInit(statsNew);
      statsNew->field = field;

      // Read in the objects for this field, setting the SET, GOOD, and OK_RUN
      // flags as appropriate.  Returns a pair of arrays of GOOD
      // objects in the field which MAY match objects in adjacent fields.
      // Each element in the "arrayNew" array contains an index into the
      // "handlesNew" array plus the object position in column and "scan row"
      // coordinates (i.e., row increases throughout the run).
      if (tsObjectStatusSet((*fp)->objectList(field), rowsPerField,
			    overlapRows, matchRadius, nSigma, &handlesNew,
			    &arrayNew, statsNew) == -1)
	{
	  // Error.  Clean up and return an error.
	  shChainCursorDel(ch, cursor);
	  taGenericChainDestroy(ch);
	  if (handlesOld != NULL) delete handlesOld;
	  if (arrayOld != NULL) shFree(arrayOld);
	  char fieldString[10];
	  sprintf(fieldString, "%d", field);
	  Tcl_AppendResult(interp, setObjectStatusDontCallCmd,
			   ": error in field ", fieldString, NULL);
	  return TCL_ERROR;
	}

      // Skip remainder if this is the first field
      if (field == field0) continue;

      // Match objects between this and the previous field.  Just search in
      // the lower overlap region +- the match radius.  Sets the DUPLICATE
      // flag for all matches and guarantees that exactly one object in each
      // matched pair has the OK_RUN flag set.
      tsSelfMerge(handlesNew, arrayNew, handlesOld, arrayOld, matchRadius,
		  field * rowsPerField - matchRadius,
		  field * rowsPerField + overlapRows + matchRadius,
		  statsNew);

      // Count the number of GOOD unmatched objects in the overlap regions for
      // the previous field (now that we've finished matching it).
      int size = handlesOld->size();
      for (int k = 0; k < size; k++)
	{
	  int status = (*handlesOld)[arrayOld[k].index]->status();
	  if (! (status & AR_OBJECT_STATUS_OK_RUN) &&
	      ! (status & AR_OBJECT_STATUS_DUPLICATE))
	    {
	      statsOld->nOverlap++;
	      if (arrayOld[k].bright) statsOld->nOverlapBright++;
	    }
	}

      // Delete info for previous field
      delete handlesOld;
      shFree(arrayOld);
    }
  shChainCursorDel(ch, cursor);
  printf("\n");
  fflush(stdout);
  
  // Count the number of GOOD unmatched objects in the overlap regions for
  // the last field (now that we've finished matching it).
  int size = handlesNew->size();
  for (int k = 0; k < size; k++)
    {
      int status = (*handlesNew)[arrayNew[k].index]->status();
      if (! (status & AR_OBJECT_STATUS_OK_RUN) &&
	  ! (status & AR_OBJECT_STATUS_DUPLICATE))
	{
	  statsNew->nOverlap++;
	  if(arrayNew[k].bright) statsNew->nOverlapBright++;
	}
    }

  // Delete info for last field
  delete handlesNew;
  shFree(arrayNew);

  // Return a handle to the matching statistics chain
  char vName[HANDLE_NAMELEN];
  if (shTclHandleNew(interp, vName, "CHAIN", (void *)ch) != TCL_OK)
    {
      // Error.  Delete the chain and return an error.
      taGenericChainDestroy(ch);
      Tcl_SetResult(interp,
		    "setObjectStatusDontCall: error allocating chain TCL handle",
		    TCL_VOLATILE);
      return TCL_ERROR;
    }
  Tcl_SetResult(interp, vName, TCL_VOLATILE);
  return TCL_OK;
}

/*============================================================================
**<AUTO EXTRACT>
**
** TCL VERB: unsetObjectStatus
**
**<HTML>
**	Undo a "setObjectStatus" operation for a single Frames Pipeline Run.
**	This clears the "arObjectStatus" flags for all objects in the
**	run.  It also removes all links for all objects in
**	the run to other objects belonging to the run.  It does not remove
**	links to objects in other Frames Pipeline Runs.
**</HTML>
**
**</AUTO>
**============================================================================
*/
static char runUnsetObjectStatusCmd[] = "unsetObjectStatus";
static int runUnsetObjectStatusFlg = FTCL_ARGV_NO_LEFTOVERS;
static ftclArgvInfo runUnsetObjectStatusTbl[] = {
   {NULL, FTCL_ARGV_HELP, NULL, NULL,
    "Reset object status to UNDETERMINED for all objects in a Frames Pipeline Run, and break links to other objects in the same run\n"},
   {"<run>", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of Frames Pipeline run"},
   {"-force", FTCL_ARGV_CONSTANT, (void *) 1, NULL,
    "Force operation, even if run hasn't been through 'setObjectStatus'"},
   {NULL, FTCL_ARGV_END, NULL, NULL, NULL}
};

static int
tclRunUnsetObjectStatus (ClientData clientData,
			 Tcl_Interp *interp,
			 int argc,
			 char **argv)
{
   ooHandle(arFramesPipeRun) *framesH = NULL;
   char *name = NULL;
   int force = 0;
   int status;

   // Parse the arguments
   runUnsetObjectStatusTbl[1].dst = &name;
   runUnsetObjectStatusTbl[2].dst = &force;
   if ((status = shTclParseArgv(interp, &argc, argv, runUnsetObjectStatusTbl,
				runUnsetObjectStatusFlg,
				runUnsetObjectStatusCmd))
       != FTCL_ARGV_SUCCESS)
     return status;
   
   // Get handle to frames pipeline run
   if (shTclAddrGetFromName (interp, name, (void **)&framesH,
			     "ooHandle(arFramesPipeRun)") != TCL_OK)
     {
       Tcl_AppendResult(interp, runUnsetObjectStatusCmd,
			": error getting handle to run", NULL);
       return TCL_ERROR;
     }
   
   // Illegal to do so for a run which has already been merged against
   // other runs.
   if ((*framesH)->merge().isValid() == oocTrue)
     {
       Tcl_AppendResult(interp, runUnsetObjectStatusCmd,
			": illegal: run has already been merged", NULL);
       return TCL_ERROR;
     }

   // Verify that this run has been through "setObjectStatus"
   if (! force)
     {
       if ((*framesH)->setObjectStatus().isValid() == oocFalse)
	 {
	   Tcl_AppendResult(interp, runUnsetObjectStatusCmd,
			    ": run has not been through 'setObjectStatus'",
			    NULL);
	   return TCL_ERROR;
	 }
     }
   
   // Fetch the range of fields
   uint16 field0 = (*framesH)->field0();
   uint16 nFields = (*framesH)->nFields();

   // For each field ...
   for (int32 field = field0; field < field0 + nFields; field++)
     {
       printf("%4d", field);
       fflush(stdout);
       // Set the status of all objects in the field to undetermined, and break
       // links to other objects in the same run.
       ooHandle(arObjectList) olist = (*framesH)->objectList(field);
       if (tsFieldUnsetObjectStatus(olist)!=SH_SUCCESS)
	 {
	   // Error.
	   char fieldString[10];
	   sprintf(fieldString, "%d", field);
	   Tcl_AppendResult(interp, runUnsetObjectStatusCmd,
			    ": error in field ", fieldString, NULL);
	   return TCL_ERROR;
	 }
     }
   printf("\n");
   fflush(stdout);

   // Remove link to "setObjectStatus" report
   ooHandle(arProduct) report = (*framesH)->setObjectStatus();
   (*framesH)->delSetObjectStatus();
   if (report.isValid() == oocTrue)
     {
       ooDelete(report);
     }

   // Return.
   return TCL_OK;
}

/*============================================================================
**<AUTO EXTRACT>
**
** TCL VERB: selfMergeQA
**
**<HTML>
**	Generate matching statistics for when "setObjectStatus" was run
**	on the specified Frames Pipeline Run.  This returns a CHAIN of
**	TSMATCHSTATS, with one set of statistics for each field in the run.
**	The returned statistics give the row and column statistics in pixels,
**	and the uncalibrated magnitude statistics (PSF for stars, petrosian
**	for galaxies) in magnitudes.
**</HTML>
**
**</AUTO>
**============================================================================
*/
static char selfMergeQACmd[] = "selfMergeQA";
static int selfMergeQAFlg = FTCL_ARGV_NO_LEFTOVERS;
static ftclArgvInfo selfMergeQATbl[] = {
   {NULL, FTCL_ARGV_HELP, NULL, NULL,
    "Generate matching statistics between adjacent fields for a Frames Pipeline run\n"},
   {"<run>", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of Frames Pipeline run"},
   {"-firstField", FTCL_ARGV_INT, NULL, NULL, "First field"},
   {"-lastField", FTCL_ARGV_INT, NULL, NULL, "Last field"},
   {"-nSigma", FTCL_ARGV_DOUBLE, NULL, NULL,
    "Compile matching statistics for nSigma detections (PSF, in r') only (defaults to 10)"},
   {"-noBlends", FTCL_ARGV_CONSTANT, (void *) 1, NULL,
    "Don't include blended objects (DFLAG_BLENDED || DFLAG_CHILD) in matching statistics and chains"},
   {NULL, FTCL_ARGV_END, NULL, NULL, NULL}
};

static int
tclSelfMergeQA (ClientData clientData,
		Tcl_Interp *interp,
		int argc,
		char **argv)
{
   ooHandle(arFramesPipeRun) *fp = NULL;
   char *name = NULL;
   int firstField = -1, lastField = -1;
   double nSigma = 10;
   int noBlends = 0;
   int status;

   // Parse the arguments
   selfMergeQATbl[1].dst = &name;
   selfMergeQATbl[2].dst = &firstField;
   selfMergeQATbl[3].dst = &lastField;
   selfMergeQATbl[4].dst = &nSigma;
   selfMergeQATbl[5].dst = &noBlends;
   if ((status = shTclParseArgv(interp, &argc, argv, selfMergeQATbl,
				selfMergeQAFlg,	selfMergeQACmd))
       != FTCL_ARGV_SUCCESS)
     return status;
   
   // Get handle to frames pipeline run
   if (shTclAddrGetFromName (interp, name, (void **)&fp,
			     "ooHandle(arFramesPipeRun)") != TCL_OK)
     {
       Tcl_AppendResult(interp, selfMergeQACmd,
			": error getting handle to run", NULL);
       return TCL_ERROR;
     }
   
   // Verify that this run has been through "setObjectStatus"
   if ((*fp)->setObjectStatus().isValid() == oocFalse)
     {
       Tcl_AppendResult(interp, selfMergeQACmd,
			": run has not been through 'setObjectStatus'", NULL);
       return TCL_ERROR;
     }
   
   // Get number of rows per field.
   int camRow = 1;  // Any photometric CCD row will do
   int rowsPerField =
     (*fp)->imagingRun()->program()->ccd(camRow,(*fp)->camCol())->nRows();

   // Get/verify field range
   int field0 = (*fp)->field0();
   int field1 = field0 + (*fp)->nFields() - 1;
   if (firstField > -1)
     {
       if (firstField < field0 || firstField > field1)
	 {
	   Tcl_AppendResult(interp, selfMergeQACmd,
			    ": 'firstField' outside range of fields for this run",
			    NULL);
	   return TCL_ERROR;
	 }
     }
   else
     {
       firstField = field0;
     }
   if (lastField > -1)
     {
       if (lastField < field0 || lastField > field1)
	 {
	   Tcl_AppendResult(interp, selfMergeQACmd,
			    ": 'lastField' outside range of fields for this run",
			    NULL);
	   return TCL_ERROR;
	 }
     }
   else
     {
       lastField = field1;
     }
   int nFields = lastField - firstField + 1;

   // Allocate a CHAIN to return the matching statistics for each field.
   CHAIN *ch = shChainNewBlock("TSMATCHSTATS", nFields);

   // For each field ...
   TSMATCHSTATS *stats;
   CURSOR_T cursor = shChainCursorNew(ch);
   for (int32 field = firstField; field <= lastField; field++)
     {
       printf("%4d", field);
       fflush(stdout);
       // Get the stats structure for this field
       if ((stats = (TSMATCHSTATS*) shChainWalk(ch, cursor, NEXT)) == NULL)
	 {
	   // Error.  Clean up and return an error.
	   shChainCursorDel(ch, cursor);
	   taGenericChainDestroy(ch);
	   char fieldString[10];
	   sprintf(fieldString, "%d", field);
	   Tcl_AppendResult(interp, selfMergeQACmd,
			    ": error getting stats in field ", fieldString,
			    NULL);
	   return TCL_ERROR;
	 }
       stats->field = field;

       // Calculate the matching statistics for this field.  Set the previous
       // field for the first field to itself --- this allows the subroutine to
       // work, but it will find no matches.
       if (tsFieldMatchStatsByPixel((*fp)->objectList(field),
				    (*fp)->objectList(field == field0 ?
						      field0 : field-1),
				    rowsPerField, nSigma, stats, noBlends)
	   != SH_SUCCESS)
	 {
	   // Error.  Delete the chain and return an error.
	   shChainCursorDel(ch, cursor);
	   taGenericChainDestroy(ch);
	   char fieldString[10];
	   sprintf(fieldString, "%d", field);
	   Tcl_AppendResult(interp, selfMergeQACmd,
			    ": error in field ", fieldString, NULL);
	   return TCL_ERROR;
	 }
     }
   shChainCursorDel(ch, cursor);
   printf("\n");
   fflush(stdout);

   // Return a handle to the matching statistics chain
   char vName[HANDLE_NAMELEN];
   if (shTclHandleNew(interp, vName, "CHAIN", (void *)ch) != TCL_OK)
     {
       // Error.  Delete the chain and return an error.
       taGenericChainDestroy(ch);
       Tcl_AppendResult(interp, selfMergeQACmd,
			": error allocating chain TCL handle", NULL);
       return TCL_ERROR;
     }
   Tcl_SetResult(interp, vName, TCL_VOLATILE);
   return TCL_OK;
}

/*============================================================================
**<AUTO EXTRACT>
**
** TCL VERB: fieldPlot
**
**<HTML>
**	Plot the objects for a range of fields in a single Frames Pipeline Run.
**	Two styles of plots are supported, the first to examine the results of
**	"setObjectStatus", and the second the result of "resolve".
**	In both plot styles, even numbered fields' objects are plotted in green,
**	odd numbered fields in red.
**	The symbol used for each point depends on the object's status.
**	In the first style of plot (the default, useful for examine the results
**	of "setObjectStatus"), 
**	triangle = matching object between fields (only one symbol is plotted
**	for the pair, plotted in color of field to which the OK_RUN object in
**	the pair belongs); dot = OK_RUN; plus = GOOD but not OK_RUN;
**	cross = BAD.  In the second style of plot, selected by specifying the
**	"-resolve" flag (useful for examining the results of "resolve"),
**	plus = unmatched primary, asterik = primary which matches one or more
**	secondaries in another run, triangle = secondary which matches a
**	primary	in another run, cross = secondary which does not match a
**	primary.  Primaries which have also been selected as spectroscopic
**	targets are also marked with a circle.
**<p>
**	In both types of plots, field boundaries are plotted in white.  For the
**	"-resolve" plots, scanline boundaries are plotted in yellow and stripe
**	boundaries in blue.
**</HTML>
**
**</AUTO>
**============================================================================
*/
static char fieldPlotCmd[] = "fieldPlot";
static int fieldPlotFlg = FTCL_ARGV_NO_LEFTOVERS;
static ftclArgvInfo fieldPlotTbl[] = {
   {NULL, FTCL_ARGV_HELP, NULL, NULL,
    "Plot objects for a single Frames Pipeline Run\n"},
   {"<run>", FTCL_ARGV_STRING, NULL, NULL, "Handle of Frames Pipeline run"},
   {"[device]", FTCL_ARGV_STRING, "/XWINDOW", NULL, "plotting device"},
   {"-firstField", FTCL_ARGV_INT, NULL, NULL,
    "First field to plot (defaults to first field in the run)"},
   {"-lastField", FTCL_ARGV_INT, NULL, NULL,
    "Last field to plot (defaults to last field in the run)"},
   {"-resolve", FTCL_ARGV_CONSTANT, (void *) 1, NULL,
    "Plot symbols for debugging resolve"},
   {NULL, FTCL_ARGV_END, NULL, NULL, NULL}
};

static int
tclFieldPlot (ClientData clientData,
	      Tcl_Interp *interp,
	      int argc,
	      char **argv)
{
   ooHandle(arFramesPipeRun) *fp = NULL;
   char *name = NULL;
   char *device = NULL;
   int firstField = -1, lastField = -1;
   int resolve = 0;
   int status;

   // Parse the arguments
   fieldPlotTbl[1].dst = &name;
   fieldPlotTbl[2].dst = &device;
   fieldPlotTbl[3].dst = &firstField;
   fieldPlotTbl[4].dst = &lastField;
   fieldPlotTbl[5].dst = &resolve;
   if ((status = shTclParseArgv(interp, &argc, argv, fieldPlotTbl,
				fieldPlotFlg,
				fieldPlotCmd))
       != FTCL_ARGV_SUCCESS)
     return status;
   
   // Get handle to frames pipeline run
   if (shTclAddrGetFromName (interp, name, (void **)&fp,
			     "ooHandle(arFramesPipeRun)") != TCL_OK)
     {
       Tcl_AppendResult(interp, fieldPlotCmd,
			": error getting handle to run", NULL);
       return TCL_ERROR;
     }
   
   // Verify that this run has been through "setObjectStatus"
   if ((*fp)->setObjectStatus().isValid() == oocFalse)
     {
       Tcl_AppendResult(interp, fieldPlotCmd,
			": run has not been through 'setObjectStatus'",
			NULL);
       return TCL_ERROR;
     }

   // If testing resolve, get astrometric calibration
   ooHandle(arAstromCalib) astrom;
   int ccd;
   double etaMin, etaMax, nuMin, nuMax, stripeNode, stripeIncl;
   if (resolve)
     {
       astrom = (*fp)->mergeAstrometricCalibration();
       if (! astrom.isValid()) astrom = (*fp)->astrometricCalibration();

       // All Frames Pipeline Runs are required when stuffed to have used r' as
       // their reference filter.  Thus, we always want the TRANS structure for
       // the r' CCD.
       int camRow = 1;
       int camCol = (*fp)->camCol();
       ccd = 10 * camRow + camCol;

       // Get the minimum and maximum survey latitude (eta) bounding this
       // stripe, and the minimum and maximum great circle latitude (nu)
       // bounding this scanline.
       int stripe = (*fp)->imagingRun()->stripe();
       double etaCenter = at_etaMin + stripe * at_stripeSeparation;
       etaMin = etaCenter - at_stripeSeparation / 2.;
       etaMax = etaCenter + at_stripeSeparation / 2.;
       int scanline = (camCol - 1) * 2 +
	 ((*fp)->imagingRun()->strip() == AR_STRIP_NORTH ? 1 : 0);
       nuMin = (scanline - 6) * at_scanSeparation;
       nuMax = nuMin + at_scanSeparation;
       stripeNode = 95.;
       stripeIncl = at_etaMin + stripe * at_stripeSeparation +
	 at_surveyCenterDec;
       int southern = (stripe > 45);
       if (southern) stripeIncl -= 180.;
     }

   // Get/verify field range
   int field0 = (*fp)->field0();
   int field1 = field0 + (*fp)->nFields() - 1;
   if (firstField > -1)
     {
       if (firstField < field0 || firstField > field1)
	 {
	   Tcl_AppendResult(interp, fieldPlotCmd,
			    ": 'firstField' outside range of fields for this run",
			    NULL);
	   return TCL_ERROR;
	 }
     }
   else
     {
       firstField = field0;
     }
   if (lastField > -1)
     {
       if (lastField < field0 || lastField > field1)
	 {
	   Tcl_AppendResult(interp, fieldPlotCmd,
			    ": 'lastField' outside range of fields for this run",
			    NULL);
	   return TCL_ERROR;
	 }
     }
   else
     {
       lastField = field1;
     }
   int nFields = lastField - firstField + 1;

  // Get number of rows and columns per field.
  int camRow = 1;  // Any photometric CCD row will do
  ooHandle(arCCDProgram) program =
    (*fp)->imagingRun()->program()->ccd(camRow, (*fp)->camCol());
  int rowsPerField = program->nRows();
  // HARDWIRED:
  int colsPerField = 2048;

  // HARDWIRED: number of overlap rows per frame.
  int overlapRows = 128;

  // Open PGPLOT window
  if (cpgbeg(0, device, 1, 1) != 1)
    {
      Tcl_AppendResult(interp, fieldPlotCmd,
		       ": couldn't open plotting device", NULL);
      return TCL_ERROR;
    }
  cpgask(0);

  // Plot the axes and label
  cpgenv(0, rowsPerField * nFields + overlapRows, -100, colsPerField+100, 1, 0);
  char label[100];
  sprintf(label, "Run %d, camCol %d, rerun %d, fields %d-%d", (*fp)->run(),
	  (*fp)->camCol(), (*fp)->rerun(), firstField, lastField);
  cpglab("x", "y", label);
  
  // Plot first field boundary
  float boundary = overlapRows/2.;
  cpgsci(1);
  cpgmove(boundary, 0);
  cpgdraw(boundary, colsPerField);

   // Plot points
   int offset = 0;
   cpgsfs(2);  // plot symbols as outlines
   for (int field = firstField; field <= lastField; field++)
     {
       // Plot field boundary
       boundary = (field-firstField+1) * rowsPerField + overlapRows/2.;
       cpgsci(1);
       cpgmove(boundary, 0);
       cpgdraw(boundary, colsPerField);

       // If plotting for resolve, plot scanline and stripe limits
       double muMinOrig;
       TRANS trans;
       if (resolve)
	 {
	   // Plot scanline limits
	   cpgsci(7); // yellow
	   trans = astrom->astrotoolsTrans(ccd, field);
	   double muMin = 0, muMax = 0, nu = 0;
	   atTransApply(&trans, 'r', 0., 0., colsPerField / 2., 0., NULL, NULL,
			&muMin, NULL, &nu, NULL);
	   atTransApply(&trans, 'r', rowsPerField + overlapRows, 0.,
			colsPerField / 2., 0., NULL, NULL, &muMax, NULL, &nu,
			NULL);
	   muMinOrig = muMin;
	   // Parked scans don't scan the stripe great circle, so they
	   // need to be converted.
	   if (fabs(astrom->node() - stripeNode) > 0.00001 ||
	       fabs(astrom->incl() - stripeIncl) > 0.00001)
	     {
	       double ra, dec;
	       atGCToEq(muMin, 0, &ra, &dec, astrom->node(), astrom->incl());
	       atEqToGC(ra, dec, &muMin, &nu, stripeNode, stripeIncl);
	       atGCToEq(muMax, 0, &ra, &dec, astrom->node(), astrom->incl());
	       atEqToGC(ra, dec, &muMax, &nu, stripeNode, stripeIncl);
	     }
	   int nSteps = 100;
	   double dMu = (muMax - muMin) / nSteps;
	   muMin -= dMu;
	   muMax += dMu;
	   dMu = (muMax - muMin) / nSteps;
	   for (int i = 0; i <= nSteps; i++)
	     {
	       double mu = muMin + i * dMu;
	       nu = nuMin;
	       // Parked scans don't scan the stripe great circle, so they
	       // need to be converted.
	       if (fabs(astrom->node() - stripeNode) > 0.00001 ||
		   fabs(astrom->incl() - stripeIncl) > 0.00001)
		 {
		   double ra, dec;
		   atGCToEq(mu, nu, &ra, &dec, stripeNode, stripeIncl);
		   atEqToGC(ra, dec, &mu, &nu, astrom->node(), astrom->incl());
		   atBound(&mu, muMinOrig, muMinOrig + 360.);
		 }
	       double row, col;
	       atTransInverseApply(&trans, 'r', mu, 0., nu, 0., NULL, NULL,
				   &row, NULL, &col, NULL);
	       if (i == 0)
		 cpgmove(row + offset, col);
	       else
		 cpgdraw(row + offset, col);
	     }
	   for (i = 0; i <= nSteps; i++)
	     {
	       double mu = muMin + i * dMu;
	       nu = nuMax;
	       // Parked scans don't scan the stripe great circle, so they
	       // need to be converted.
	       if (fabs(astrom->node() - stripeNode) > 0.00001 ||
		   fabs(astrom->incl() - stripeIncl) > 0.00001)
		 {
		   double ra, dec;
		   atGCToEq(mu, nu, &ra, &dec, stripeNode, stripeIncl);
		   atEqToGC(ra, dec, &mu, &nu, astrom->node(), astrom->incl());
		   atBound(&mu, muMinOrig, muMinOrig+360.);
		 }
	       double row, col;
	       atTransInverseApply(&trans, 'r', mu, 0., nu, 0., NULL, NULL,
				   &row, NULL, &col, NULL);
	       if (i == 0)
		 cpgmove(row + offset, col);
	       else
		 cpgdraw(row + offset, col);
	     }

	   // Plot stripe limits
	   cpgsci(4); // blue
	   double lambdaMin, lambdaMax, eta;
	   atGCToSurvey(muMin, 0, stripeNode, stripeIncl, &lambdaMin, &eta);
	   atGCToSurvey(muMax, 0, stripeNode, stripeIncl, &lambdaMax, &eta);
	   double dLambda = (lambdaMax - lambdaMin) / nSteps;
	   for (i = 0; i <= nSteps; i++)
	     {
	       double lambda = lambdaMin + i * dLambda;
	       double mu, nu;
	       atSurveyToGC(lambda, etaMin, astrom->node(), astrom->incl(),
			    &mu, &nu);
	       atBound(&mu, muMinOrig, muMinOrig+360.);
	       double row, col;
	       atTransInverseApply(&trans, 'r', mu, 0., nu, 0., NULL, NULL,
				   &row, NULL, &col, NULL);
	       if (i == 0)
		 cpgmove(row + offset, col);
	       else
		 cpgdraw(row + offset, col);
	     }
	   for (i = 0; i <= nSteps; i++)
	     {
	       double lambda = lambdaMin + i * dLambda;
	       double mu, nu;
	       atSurveyToGC(lambda, etaMax, astrom->node(), astrom->incl(),
			    &mu, &nu);
	       atBound(&mu, muMinOrig, muMinOrig+360.);
	       double row, col;
	       atTransInverseApply(&trans, 'r', mu, 0., nu, 0., NULL, NULL,
				   &row, NULL, &col, NULL);
	       if (i == 0)
		 cpgmove(row + offset, col);
	       else
		 cpgdraw(row + offset, col);
	     }
	 }

       // Get object list
       ooHandle(arObjectList) list = (*fp)->objectList(field);

       // Plot segment boundary for each segments for which this list is the
       // first or last field in the segment.
       if (resolve)
	 {
	   ooItr(arSegment) segments;
	   list->segments(segments);
	   while (segments.next())
	     {
	       if (field == segments->field0() ||
		   field == segments->field0() + segments->nFields() - 1)
		 {
		   cpgsci(7);
		   double mu;
		   if (field == segments->field0())
		     mu = segments->startMu() * at_asec2Deg;
		   else
		     mu = segments->endMu() * at_asec2Deg;
		   int nSteps = 100;
		   double dNu = (nuMax - nuMin) / nSteps;
		   for (int i = 0; i <= nSteps; i++)
		     {
		       double mu1 = mu;
		       double nu1 = nuMin + i * dNu;
		       if (fabs(astrom->node() - stripeNode) > 0.00001 ||
			   fabs(astrom->incl() - stripeIncl) > 0.00001)
			 {
			   double ra, dec;
			   atGCToEq(mu1, nu1, &ra, &dec, stripeNode,stripeIncl);
			   atEqToGC(ra, dec, &mu1, &nu1, astrom->node(),
				    astrom->incl());
			   atBound(&mu1, muMinOrig, muMinOrig+360.);
			 }
		       double row, col;
		       atTransInverseApply(&trans, 'r', mu1, 0., nu1, 0., NULL,
					   NULL, &row, NULL, &col, NULL);
		       if (i == 0)
			 cpgmove(row + offset, col);
		       else
			 cpgdraw(row + offset, col);
		     }
		 }
	     }
	 }

       // Set color of objects in this and next field.
       int thisColor, nextColor;
       if ((field / 2) * 2 == field)
	 {
	   thisColor = 2;
	   nextColor = 3;
	 }
       else
	 {
	   thisColor = 3;
	   nextColor = 2;
	 }

       // Get iterator to objects for this field
       ooItr(arObject) itr;
       if (list->objects(itr, oocReadOnly) != oocSuccess)
	 {
	   Tcl_SetResult(interp,
			 "selfMergePlot: couldn't initialize iterator over objects",
			 TCL_VOLATILE);
	   return TCL_ERROR;
	 }

       // For each object ...
       while (itr.next())
	 {
	   // Set color to default color for this field.
	   cpgsci(thisColor);

	   // Get row and col
	   float x = itr->objc_rowc().value() + offset;
	   float y = itr->objc_colc().value();

	   // Get object status
	   int status = itr->status();

	   int symbol;
	   int secondSymbol = 0;
	   if (resolve)
	     {
	       // Only plot OK_RUN objects
	       if (! (status & AR_OBJECT_STATUS_OK_RUN)) continue;

	       // Plot unresolved objects as dots
	       if (! (status & AR_OBJECT_STATUS_RESOLVED))
		 symbol = 1;
	       else
		 {
		   // Check for matching objects.
		   int matchStatus = -1;
		   ooItr(arObject) matches;
		   itr->matches(matches);
		   while (matches.next())
		     {
		       if (matches->status() & AR_OBJECT_STATUS_PRIMARY)
			 {
			   matchStatus = 1;
			   break;
			 }
		       else if (matches->status() & AR_OBJECT_STATUS_SECONDARY)
			 matchStatus = 0;
		     }

		   // Plot unmatched primaries as a plus, primaries which
		   // match a secondary as an asterik, secondaries which match
		   // a primary as a triangle, unmatched secondaries as a
		   // cross, and objects which are neither primary nor
		   // secondary as a dot.
		   if (status & AR_OBJECT_STATUS_PRIMARY)
		     {
		       if (matchStatus == 1)
			 {
			   Tcl_AppendResult(interp, fieldPlotCmd,
					    ": two matched primaries", NULL);
			   return TCL_ERROR;
			 }
		       else if (matchStatus == 0)
			 symbol = 3;
		       else
			 symbol = 2;

		       // If a spectroscopic target, mark with a circle if a
		       // QSO, else a square if a galaxy, else a diamond.
		       if (status & AR_OBJECT_STATUS_TARGET)
			 {
			   ooHandle(arTarget) target;
			   target = itr->spectroscopicTarget(target);
			   uint32 primaryType = target->primaryType();
			   if (primaryType &
			       (AR_TARGET_QSO_HIZ | AR_TARGET_QSO_CAP |
				AR_TARGET_QSO_SKIRT | AR_TARGET_QSO_FIRST_CAP |
				AR_TARGET_QSO_FIRST_SKIRT))
			     secondSymbol = 4;
			   else if (primaryType &
				    (AR_TARGET_GALAXY_RED | AR_TARGET_GALAXY |
				     AR_TARGET_GALAXY_BIG |
				     AR_TARGET_GALAXY_BRIGHT_CORE))
			     secondSymbol = 6;
			   else
			     secondSymbol = 11;
			 }
		     }
		   else if (status & AR_OBJECT_STATUS_SECONDARY)
		     {
		       if (matchStatus == 1)
			 symbol = 7;
		       else
			 symbol = 5;
		     }
		   else
		     {
		       symbol = 1;
		     }
		 }
	     }
	   else
	     {
	       if (! (status & AR_OBJECT_STATUS_GOOD))
		 symbol = 5;   // plot bad objects as crosses
	       else if (status & AR_OBJECT_STATUS_DUPLICATE)
		 {
		   // Plot matching objects just once.   Only plot when in
		   // upper overlap region, except in first field.  Plot as
		   // triangles.
		   if (itr->objc_rowc().value() < rowsPerField / 2. &&
		       field != firstField)
		     continue;
		   symbol = 7;
		   if (! (status & AR_OBJECT_STATUS_OK_RUN))
		     cpgsci(nextColor);
		 }
	       else if (status & AR_OBJECT_STATUS_OK_RUN)
		 symbol = 1;   // plot OK_RUN objects as dots
	       else
		 symbol = 2;   // plot non-matching !OK_RUN objects as pluses
	     }
	       
	   // Plot the point
	   cpgpt(1, &x, &y, symbol);
	   if (secondSymbol > 0) cpgpt(1, &x, &y, secondSymbol);
	}
      offset = offset + rowsPerField;
     }
   cpgask(1);

   // Return nothing
   Tcl_SetResult(interp, "", TCL_VOLATILE);
   return TCL_OK;
}

/*============================================================================
**<AUTO EXTRACT>
**
** TCL VERB: xRunMergeDontCall
**
**<HTML>
**	Merge the OK_RUN objects (only those classified as STAR, GALAXY, or
**	UNK) in a Frames Pipeline Run with OK_RUN stars and galaxies in all
**	other previously merged Frames Pipeline Runs.  Objects match if their
**	separation is less than or
**	equal to the specified match 'radius'.  If an object in a single
**	previously-merged Frames Pipeline Run matches more than one object in
**	the target run, then only the closest match is merged (an object in
**	the target run can match more than one object in each of the
**	previously-merged run).
**<p>
**	The specified astrometric calibration is used to derive calibrated
**	positions for matching.
**	DCR corrections are applied using cosmic colors, since
**	a final photometric calibration can not be guaranteed to exist for
**	the run at this time.
**<P>
**	This verb should NOT be called directly.  The verb "runMerge" should
**      should be used; it calls this verb.
**</HTML>
**
**</AUTO>
**============================================================================
*/
static char runMergeDontCallCmd[] = "xRunMergeDontCall";
static int runMergeDontCallFlg = FTCL_ARGV_NO_LEFTOVERS;
static ftclArgvInfo runMergeDontCallTbl[] = {
   {NULL, FTCL_ARGV_HELP, NULL, NULL,
    "Merge the objects in a Frames Pipeline Run with objects in all other previously merged Frames Pipeline Runs\n"},
   {"<run>", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of Frames Pipeline run (ooHandle(arFramesPipeRun))"},
   {"<astrom>", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of astrometric calibration to use (ooHandle(arAstromCalib))"},
   {"<matchRadius>", FTCL_ARGV_DOUBLE, NULL, NULL,
    "Match 'radius' (arcsecs)"},
   {NULL, FTCL_ARGV_END, NULL, NULL, NULL}
};

static int
tclRunMergeDontCall (ClientData clientData,
	     Tcl_Interp *interp,
	     int argc,
	     char **argv)
{
  ooHandle(arFramesPipeRun) *framesH = NULL;
  ooHandle(arAstromCalib) *astromH = NULL;
  char *framesName = NULL;
  char *astromName = NULL;
  double matchRadius;
  int status;

  // Parse the arguments
  runMergeDontCallTbl[1].dst = &framesName;
  runMergeDontCallTbl[2].dst = &astromName;
  runMergeDontCallTbl[3].dst = &matchRadius;
  if ((status = shTclParseArgv(interp, &argc, argv, runMergeDontCallTbl,
			       runMergeDontCallFlg, runMergeDontCallCmd))
      != FTCL_ARGV_SUCCESS)
    return status;
   
  // Get handle to frames pipeline run
  if (shTclAddrGetFromName (interp, framesName, (void **)&framesH,
			    "ooHandle(arFramesPipeRun)") != TCL_OK)
    {
      Tcl_AppendResult(interp, runMergeDontCallCmd,
		       ": error getting handle to frames pipeline run", NULL);
      return TCL_ERROR;
    }
   
  // Verify that this run has not previously been merged
  if ((*framesH)->merge().isValid() == oocTrue)
    {
      Tcl_AppendResult(interp, runMergeDontCallCmd, ": run '", framesName,
		       "' has previously been merged", NULL);
      return TCL_ERROR;
    }

  // Verify that this run has been through setObjectStatus.
  if ((*framesH)->setObjectStatus().isValid() == oocFalse)
    {
      Tcl_AppendResult(interp, runMergeDontCallCmd, ": run '", framesName,
		       "' has not been through setObjectStatus", NULL);
      return TCL_ERROR;
    }

  // Get handle to astrometric calibration
  if (shTclAddrGetFromName (interp, astromName, (void **)&astromH,
			    "ooHandle(arAstromCalib)") != TCL_OK)
    {
      Tcl_AppendResult(interp, runMergeDontCallCmd,
		       ": error getting handle to astrometric calibration",
		       NULL);
      return TCL_ERROR;
    }
   
  // Verify that astrometric calibration is for the target imaging run and
  // spans the full range of fields
  if ((*astromH)->run() != (*framesH)->run())
    {
      Tcl_AppendResult(interp, runMergeDontCallCmd,
		       ": astrometric calibration for wrong imaging run",
		       NULL);
      return TCL_ERROR;
    }
  int astromField0 = (*astromH)->field0();
  int framesField0 = (*framesH)->field0();
  if (astromField0 > framesField0 ||
      astromField0+(*astromH)->nFields() < framesField0+(*framesH)->nFields())
    {
      Tcl_AppendResult(interp, runMergeDontCallCmd,
		       ": astrometric calibration doesn't cover field range of the frames pipeline run",
		       NULL);
      return TCL_ERROR;
    }

  // Fetch the global TCL variable "photo_frames" for container with all
  // Frames Pipeline Runs.  Then fetch the associated handle.
  char *photo_frames = Tcl_GetVar(interp, "photo_frames", TCL_GLOBAL_ONLY);
  if (photo_frames == NULL)
    {
      Tcl_AppendResult(interp, runMergeDontCallCmd,
		       ": global variable $photo_frames doesn't exist", NULL);
      return TCL_ERROR;
    }
  ooHandle(ooContObj) *framesContH;
  if (shTclAddrGetFromName (interp, photo_frames, (void **)&framesContH,
			    "ooHandle(ooContObj)") != TCL_OK)
    {
      Tcl_AppendResult(interp, runMergeDontCallCmd,
		       ": error getting handle to $photo_frames", NULL);
      return TCL_ERROR;
    }

  // Fetch great circle parameters for target run
  double equinox = (*astromH)->equinox();
  double node = (*astromH)->node();
  double incl = (*astromH)->incl();

  // Initialize iterator over all frames pipeline runs
  ooItr(arFramesPipeRun) itr;
  if (itr.scan(*framesContH, oocRead) != oocSuccess)
    {
      Tcl_AppendResult(interp, runMergeDontCallCmd,
		       ": couldn't initalize iterator over frames pipeline runs",
		       NULL);
      return TCL_ERROR;
    }

  // Set up cosmic color for use by atTransApply().
  taTransCosmicSet(cosmicColor, cosmicScatter);

  // Declare some needed variables
  TSCALIBOBJ *array;
  ooTVArray(ooHandle(arObject)) *handles;
  int nObjects;
  double minNu, maxNu;

  // For each frames pipeline run ...
  int first = 1;
  while (itr.next())
    {
      // Skip if this is the same run as the target run
      if (itr == *framesH) continue;

      // Skip if this run has not yet been merged itself.
      ooHandle(arProduct) report;
      report = itr->merge(report);
      if (report.isValid() == oocFalse) continue;

      // TODO: we should test whether this run could possibly overlap
      // the target run, and skip it if it can't.

      // If we haven't, fetch the array of sorted (by increasing mu)
      // calibrated objects for the target run.
      if (first == 1)
	{
	  if ((nObjects = tsRunToCalibArray(*framesH, *astromH,
					    &handles, &array, &minNu, &maxNu))
	      == -1)
	    {
	      shTclInterpAppendWithErrStack(interp);
	      return TCL_ERROR;
	    }
	  if (nObjects == 0)
	    {
	      delete handles;
	      shFree(array);
	      Tcl_SetResult(interp, "", TCL_VOLATILE);
	      return TCL_OK;
	    }
	  first = 0;
	}

      // Merge the target run against this run
      if (tsRunMergeOne(handles, array, nObjects, minNu, maxNu, node, incl,
			equinox, itr, matchRadius) != SH_SUCCESS)
	{
	  delete handles;
	  shFree(array);
	  shTclInterpAppendWithErrStack(interp);
	  return TCL_ERROR;
	}
    }
  
  // Delete the array
  if (first == 0)
    {
      delete handles;
      shFree(array);
    }

  // Set the calibrations used
  (*framesH)->setMergeAstrometricCalibration(*astromH);

  // Return
  Tcl_SetResult(interp, "", TCL_VOLATILE);
  return TCL_OK;
}

/*============================================================================
**<AUTO EXTRACT>
**
** TCL VERB: runUnmerge
**
**<HTML>
**	Undo a "merge" operation of a single Frames Pipeline run.  This
**	unmerges a Frames Pipeline run from other merged Frames Pipeline runs.
**	It breaks all links between objects in the Frames Pipeline run with
**	objects in other Frames Pipeline runs.  It does not remove
**	links to objects in the same Frames Pipeline Run.
**</HTML>
**
**</AUTO>
**============================================================================
*/
static char runUnmergeCmd[] = "runUnmerge";
static int runUnmergeFlg = FTCL_ARGV_NO_LEFTOVERS;
static ftclArgvInfo runUnmergeTbl[] = {
   {NULL, FTCL_ARGV_HELP, NULL, NULL,
    "Remove links to objects in other Frames Pipeline runs\n"},
   {"<run>", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of Frames Pipeline run"},
   {"-force", FTCL_ARGV_CONSTANT, (void *) 1, NULL,
    "Force operation, even if run hasn't been through 'merge'"},
   {NULL, FTCL_ARGV_END, NULL, NULL, NULL}
};

static int
tclRunUnmerge (ClientData clientData,
	       Tcl_Interp *interp,
	       int argc,
	       char **argv)
{
   ooHandle(arFramesPipeRun) *framesH = NULL;
   char *name = NULL;
   int force = 0;
   int status;

   // Parse the arguments
   runUnmergeTbl[1].dst = &name;
   runUnmergeTbl[2].dst = &force;
   if ((status = shTclParseArgv(interp, &argc, argv, runUnmergeTbl,
				runUnmergeFlg,
				runUnmergeCmd))
       != FTCL_ARGV_SUCCESS)
     return status;
   
   // Get handle to frames pipeline run
   if (shTclAddrGetFromName (interp, name, (void **)&framesH,
			     "ooHandle(arFramesPipeRun)") != TCL_OK)
     {
       Tcl_AppendResult(interp, runUnmergeCmd,
			": error getting handle to run", NULL);
       return TCL_ERROR;
     }
   
   // Verify that this run does not belong to any resolved segments.
   ooItr(arSegment) itr;
   if ((*framesH)->segments(itr) == oocError)
     {
       Tcl_AppendResult(interp, runUnmergeCmd,
			": error initialize iterator over segments", NULL);
       return TCL_ERROR;
     }
   if (itr.next() == oocTrue)
     {
       Tcl_AppendResult(interp, runUnmergeCmd,
			": run belongs to one or more resolved segments", NULL);
       return TCL_ERROR;
     }

   // Verify that this run has been through "merge"
   ooHandle(arProduct) mergeReport = (*framesH)->merge();
   if (! force)
     {
       if (mergeReport.isValid() == oocFalse)
	 {
	   Tcl_AppendResult(interp, runUnmergeCmd,
			    ": run has not been merged", NULL);
	   return TCL_ERROR;
	 }
     }
   
   // Fetch the range of fields
   uint16 field0 = (*framesH)->field0();
   uint16 nFields = (*framesH)->nFields();

   // For each field ...
   for (int32 field = field0; field < field0 + nFields; field++)
     {
       printf("%4d", field);
       fflush(stdout);
       // Remove links to objects in other runs.
       ooHandle(arObjectList) olist = (*framesH)->objectList(field);
       if (tsFieldUnmerge(olist) != SH_SUCCESS)
	 {
	   // Error.
	   char fieldString[10];
	   sprintf(fieldString, "%d", field);
	   Tcl_AppendResult(interp, runUnmergeCmd,
			    ": error in field ", fieldString, NULL);
	   return TCL_ERROR;
	 }
     }
   printf("\n");
   fflush(stdout);

   // Remove link to "merge" report and the calibrations used for the merge,
   // and delete the merge report from the database.
   (*framesH)->delMerge();
   (*framesH)->delMergeAstrometricCalibration();
   if (mergeReport.isValid() == oocTrue)
     {
       ooDelete(mergeReport);
     }

   // Return.
   return TCL_OK;
}

/*============================================================================
**<AUTO EXTRACT>
**
** TCL VERB: runMergeQA
**
**<HTML>
**	Generate matching statistics for the matching between two Frames
**	Pipeline runs.  This returns a CHAIN of TSMATCHSTAT, with one set of
**	statistics for each field in the first run.  If the astrometric or
**	photometric calibrations to use for either run are not specified, then
**	the "current" calibrations are used.
**	Position statistics are in calibrated great circle coordinates
**	(arcsecs, where the coordinate system of "run1" is used) and magnitude
**	statistics are calibrated PSF magnitudes for stars, petrosian
**	magnitudes for galaxies.  Only GOOD objects
**	of class STAR or GALAXY are used, with separate statistics
**	calculated for stars and galaxies.
**<P>
**	If the "-obj" flag is specified, then two chains containing each
**	matched pair are returned, one for stars and one for galaxies.  In this
**	case, the procedure returns a TCL list, where the first element is
**	the stats chain, the second element is the matched stars chain, and
**	the third element is the matched galaxies chain.
**</HTML>
**
**</AUTO>
**============================================================================
*/
static char runMergeQACmd[] = "runMergeQA";
static int runMergeQAFlg = FTCL_ARGV_NO_LEFTOVERS;
static ftclArgvInfo runMergeQATbl[] = {
   {NULL, FTCL_ARGV_HELP, NULL, NULL,
    "Generate matching statistics between two Frames Pipeline runs\n"},
   {"<run1>", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of first Frames Pipeline run (ooHandle(arFramesPipeRun))"},
   {"<run2>", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of second Frames Pipeline run (ooHandle(arFramesPipeRun))"},
   {"-astrom1", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of astrometric calibration to use for run 1 (ooHandle(arAstromCalib))"},
   {"-astrom2", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of astrometric calibration to use for run 2 (ooHandle(arAstromCalib))"},
   {"-photo1", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of photometric calibration to use for run 1 (ooHandle(arFinalPhotoCalib))"},
   {"-photo2", FTCL_ARGV_STRING, NULL, NULL,
    "Handle of photometric calibration to use for run 2 (ooHandle(arFinalPhotoCalib))"},
   {"-firstField", FTCL_ARGV_INT, NULL, NULL, "First field"},
   {"-lastField", FTCL_ARGV_INT, NULL, NULL, "Last field"},
   {"-min", FTCL_ARGV_DOUBLE, "15", NULL,
    "Minimum magnitude for 'bright' objects"},
   {"-max", FTCL_ARGV_DOUBLE, "19", NULL,
    "Maximum magnitude for 'bright' objects"},
   {"-filter", FTCL_ARGV_STRING, "r", NULL,
    "Filter used for 'bright' magnitude test (ugriz)"},
   {"-obj", FTCL_ARGV_CONSTANT, (void *) 1, NULL,
    "Return chains for matched stars and galaxies"},
   {"-full", FTCL_ARGV_CONSTANT, (void *) 1, NULL,
    "If returning chains for matched stars and galaxies, return the full objects, rather than just mags and positions"},
   {"-noBlends", FTCL_ARGV_CONSTANT, (void *) 1, NULL,
    "Don't include blended objects (DFLAG_BLENDED || DFLAG_CHILD) in matching statistics and chains"},
   {NULL, FTCL_ARGV_END, NULL, NULL, NULL}
};

static int
tclRunMergeQA (ClientData clientData,
			 Tcl_Interp *interp,
			 int argc,
			 char **argv)
{
   ooHandle(arFramesPipeRun) *fp[2] = {NULL, NULL};
   ooHandle(arAstromCalib) *astrom[2] = {NULL, NULL};
   ooHandle(arAstromCalib) astromH[2];
   ooHandle(arFinalPhotoCalib) *final[2] = {NULL, NULL};
   ooHandle(arFinalPhotoCalib) finalH[2];
   char *fpName[2] = {NULL, NULL};
   char *astromName[2] = {NULL, NULL};
   char *finalName[2] = {NULL, NULL};
   int firstField = -1, lastField = -1;
   double minMag, maxMag;
   char *filter;
   int status;
   double expTime[2];
   int objChains = 0;
   int full = 0;
   int noBlends = 0;

   // Parse the arguments
   runMergeQATbl[1].dst = &fpName[0];
   runMergeQATbl[2].dst = &fpName[1];
   runMergeQATbl[3].dst = &astromName[0];
   runMergeQATbl[4].dst = &astromName[1];
   runMergeQATbl[5].dst = &finalName[0];
   runMergeQATbl[6].dst = &finalName[1];
   runMergeQATbl[7].dst = &firstField;
   runMergeQATbl[8].dst = &lastField;
   runMergeQATbl[9].dst = &minMag;
   runMergeQATbl[10].dst = &maxMag;
   runMergeQATbl[11].dst = &filter;
   runMergeQATbl[12].dst = &objChains;
   runMergeQATbl[13].dst = &full;
   runMergeQATbl[14].dst = &noBlends;
   if ((status = shTclParseArgv(interp, &argc, argv, runMergeQATbl,
				runMergeQAFlg, runMergeQACmd))
       != FTCL_ARGV_SUCCESS)
     return status;
   
   // For each input frames pipeline run ...
   for (int i = 0; i < 2; i++)
     {
       // Get handle to the run
       if (shTclAddrGetFromName (interp, fpName[i], (void **)&fp[i],
				 "ooHandle(arFramesPipeRun)") != TCL_OK)
	 {
	   Tcl_AppendResult(interp, runMergeQACmd,
			    ": error getting handle to run '", fpName[i],
			    ",", NULL);
	   return TCL_ERROR;
	 }
   
       // Verify that this run has been merged
       if ((*fp[i])->merge().isValid() == oocFalse)
	 {
	   Tcl_AppendResult(interp, runMergeQACmd, ": run '", fpName[i],
			    "' has not been merged", NULL);
	   return TCL_ERROR;
	 }
       
       // Get handle to the astrometric calibration
       char buffer[10];
       sprintf(buffer, "-astrom%d", i+1);
       if (ftcl_ArgIsPresent(argc, argv, buffer, runMergeQATbl) == 1)
	 {
	   // If one was specified, use it.
	   if (shTclAddrGetFromName (interp, astromName[i],
				     (void **)&astrom[i],
				     "ooHandle(arAstromCalib)") != TCL_OK)
	     {
	       Tcl_AppendResult(interp, runMergeQACmd,
				": error getting handle to astrom run '",
				astromName[i], "'", NULL);
	       return TCL_ERROR;
	     }
	   if ((*astrom[i])->run() != (*fp[i])->run())
	     {
	       Tcl_AppendResult(interp, runMergeQACmd,
				": astrometric calibration '",astromName[i],
				"' is for wrong run", NULL);
	       return TCL_ERROR;
	     }
	 }
       else
	 {
	   // If none was specified, use the current one.
	   (*fp[i])->currentAstrometricCalibration(astromH[i], oocRead);
	   astrom[i] = &astromH[i];
	 }

       // Get handle to the final photometric calibration
       sprintf(buffer, "-photo%d", i+1);
       if (ftcl_ArgIsPresent(argc, argv, buffer, runMergeQATbl) == 1)
	 {
	   // If one was specified, use it.
	   if (shTclAddrGetFromName (interp, finalName[i], (void **)&final[i],
				     "ooHandle(arFinalPhotoCalib)")
	       != TCL_OK)
	     {
	       Tcl_AppendResult(interp, runMergeQACmd,
				": error getting handle to final photometric calibration run '",
				finalName[i], "'", NULL);
	       return TCL_ERROR;
	     }
	   if ((*final[i])->framesPipelineRun() != (*fp[i]))
	     {
	       Tcl_AppendResult(interp, runMergeQACmd,
				": photometric calibration '",finalName[i],
				"' is for wrong run", NULL);
	       return TCL_ERROR;
	     }
	 }
       else
	 {
	   // If none was specified, use the current one.
	   (*fp[i])->currentPhotometricCalibration(finalH[i], oocRead);
	   if (finalH[i].isValid() == oocFalse)
	     {
	       Tcl_AppendResult(interp, runMergeQACmd, ": run '", fpName[i],
				"' has not been through final photometric calibration", NULL);
	       return TCL_ERROR;
	     }
	   final[i] = &finalH[i];
	 }

       // Get the exposure time (seconds)
       // HARDWIRE: All photometric chips have 2048 physical rows.
       int nRowsPhysical = 2048;
       expTime[i] = nRowsPhysical * (*fp[i])->imagingRun()->c_obs() * 1.e-6;
     }

   // Great circles must be for the same equinox.
   if (abs((*astrom[0])->equinox() - (*astrom[1])->equinox()) > 0.000001)
     {
       Tcl_AppendResult(interp, runMergeQACmd, ": mismatched equinoxes", NULL);
       return TCL_ERROR;
     }

   // Set up cosmic color for use by atTransApply().
   taTransCosmicSet(cosmicColor, cosmicScatter);

   // Get/verify field range
   int field0 = (*fp[0])->field0();
   int field1 = field0 + (*fp[0])->nFields() - 1;
   if (firstField > -1)
     {
       if (firstField < field0 || firstField > field1)
	 {
	   Tcl_AppendResult(interp, selfMergeQACmd,
			    ": 'firstField' outside range of fields for this run",
			    NULL);
	   return TCL_ERROR;
	 }
     }
   else
     {
       firstField = field0;
     }
   if (lastField > -1)
     {
       if (lastField < field0 || lastField > field1)
	 {
	   Tcl_AppendResult(interp, selfMergeQACmd,
			    ": 'lastField' outside range of fields for this run",
			    NULL);
	   return TCL_ERROR;
	 }
     }
   else
     {
       lastField = field1;
     }
   int nFields = lastField - firstField + 1;

   // Fetch the database containing the objects for the second run, and the
   // camera column.
   ooHandle(ooDBObj) db2;
   (*fp[1])->objectList((*fp[1])->field0())->objectCont().containedIn(db2);
   uint8 camCol = (*fp[1])->camCol();

   // Allocate a CHAIN to return the matching statistics for each field in the
   // first run
   CHAIN *ch = shChainNewBlock("TSMATCHSTATS", nFields);

   // If object chains are desired, allocate them.
   CHAIN *starChain = NULL;
   CHAIN *galaxyChain = NULL;
   if (objChains)
     {
       starChain = shChainNew("TSMATCHEDPAIR");
       galaxyChain = shChainNew("TSMATCHEDPAIR");
     }

   // For each field ...
   TSMATCHSTATS *stats;
   CURSOR_T cursor = shChainCursorNew(ch);
   for (int32 field = firstField; field <= lastField; field++)
     {
       printf("%4d", field);
       fflush(stdout);
       // Get the stats structure for this field
       if ((stats = (TSMATCHSTATS*) shChainWalk(ch, cursor, NEXT)) == NULL)
	 {
	   // Error.  Delete the chain and return an error.
	   taGenericChainDestroy(ch);
	   if (objChains)
	     {
	       taGenericChainDestroy(starChain);
	       taGenericChainDestroy(galaxyChain);
	     }
	   char fieldString[10];
	   sprintf(fieldString, "%d", field);
	   Tcl_AppendResult(interp, runMergeQACmd,
			    ": can't get stats for field ", fieldString, NULL);
	   return TCL_ERROR;
	 }
       stats->field = field;

       // Calculate the matching statistics for this field.
       if (tsFieldMatchStatsByArcsec((*fp[0])->objectList(field), *astrom[0],
				     *final[0], expTime[0], db2, camCol,
				     *astrom[1], *final[1], expTime[1],
				     minMag, maxMag, filter[0], stats,
				     starChain, galaxyChain, full, noBlends)
	   != SH_SUCCESS)
	 {
	   // Error.  Delete the chain and return an error.
	   taGenericChainDestroy(ch);
	   if (objChains)
	     {
	       taGenericChainDestroy(starChain);
	       taGenericChainDestroy(galaxyChain);
	     }
	   char fieldString[10];
	   sprintf(fieldString, "%d", field);
	   Tcl_AppendResult(interp, runMergeQACmd,
			    ": error in field ", fieldString, NULL);
	   return TCL_ERROR;
	 }
     }
   printf("\n");
   fflush(stdout);

   // Return a handle to the matching statistics chain
   char vName[HANDLE_NAMELEN];
   if (shTclHandleNew(interp, vName, "CHAIN", (void *)ch) != TCL_OK)
     {
       // Error.  Delete the chain and return an error.
       taGenericChainDestroy(ch);
       if (objChains)
	 {
	   taGenericChainDestroy(starChain);
	   taGenericChainDestroy(galaxyChain);
	 }
       Tcl_AppendResult(interp, runMergeQACmd,
			": error allocating chain TCL handle", NULL);
       return TCL_ERROR;
     }
   Tcl_AppendElement(interp, vName);
   if (objChains)
     {
       if (shTclHandleNew(interp, vName, "CHAIN", (void *)starChain) != TCL_OK)
	 {
	   // Error.  Delete the chain and return an error.
	   taGenericChainDestroy(ch);
	   taGenericChainDestroy(starChain);
	   taGenericChainDestroy(galaxyChain);
	   Tcl_AppendResult(interp, runMergeQACmd,
			    ": error allocating chain TCL handle", NULL);
	   return TCL_ERROR;
	 }
       Tcl_AppendElement(interp, vName);
       if (shTclHandleNew(interp, vName, "CHAIN", (void *)galaxyChain) != TCL_OK)
	 {
	   // Error.  Delete the chain and return an error.
	   taGenericChainDestroy(ch);
	   taGenericChainDestroy(starChain);
	   taGenericChainDestroy(galaxyChain);
	   Tcl_AppendResult(interp, runMergeQACmd,
			    ": error allocating chain TCL handle", NULL);
	   return TCL_ERROR;
	 }
       Tcl_AppendElement(interp, vName);
     }
   return TCL_OK;
}

/*============================================================================
**<AUTO EXTRACT>
**
** TCL VERB: matchedPairChainDel
**
**<HTML>
**	Delete a TSMATCHPAIR chain and all its entries.
**</HTML>
**
**</AUTO>
**============================================================================
*/
static char matchedPairChainDelCmd[] = "matchedPairChainDel";
static int matchedPairChainDelFlg = FTCL_ARGV_NO_LEFTOVERS;
static ftclArgvInfo matchedPairChainDelTbl[] = {
   {NULL, FTCL_ARGV_HELP, NULL, NULL,
    "Delete a TSMATCHPAIR chain and all its entries\n"},
   {"<chain>", FTCL_ARGV_STRING, NULL, NULL,
    "Handle to the chain"},
   {NULL, FTCL_ARGV_END, NULL, NULL, NULL}
};

static int
tclMatchedPairChainDel (ClientData clientData,
			Tcl_Interp *interp,
			int argc,
			char **argv)
{
  CHAIN *chain = NULL;
  char *chainName = NULL;
  CURSOR_T crsr;
  TSMATCHEDPAIR *elem = NULL;
  int status;

  // Parse the arguments
  matchedPairChainDelTbl[1].dst = &chainName;
  if ((status = shTclParseArgv(interp, &argc, argv, matchedPairChainDelTbl,
			       matchedPairChainDelFlg, matchedPairChainDelCmd))
      != FTCL_ARGV_SUCCESS)
    return status;

   
  // Translate the handle to the CHAIN
  if (shTclAddrGetFromName (interp, chainName, (void **)&chain,
			    "CHAIN") != TCL_OK)
    {
      Tcl_SetResult(interp,"arg does not appear to be a CHAIN!",TCL_STATIC);
      return(TCL_ERROR);
    }


  // Check that we've got the right CHAIN type
  if (shChainTypeGet(chain) != shTypeGetFromName ("TSMATCHEDPAIR"))
    {
      Tcl_SetResult(interp,"wrong chain type for chain!",TCL_STATIC);
      return(TCL_ERROR);
    }

  // Done with all the handles, etc. Now the action begins
  crsr = shChainCursorNew(chain);
  while ((elem = (TSMATCHEDPAIR*) shChainWalk(chain, crsr, NEXT)) != NULL)
    {
      if (elem->obj1 != NULL) shFree(elem->obj1);
      if (elem->obj2 != NULL) shFree(elem->obj2);
      shFree(elem);
    }
  shChainCursorDel(chain,crsr);
  shChainDel(chain);
  p_shTclHandleDel(interp, chainName); 

  Tcl_SetResult(interp,"",TCL_STATIC);
  return TCL_OK;
}

/*============================================================================
**<AUTO EXTRACT>
**
** TCL VERB: mergeQABin
**
**<HTML>
**	Bin a TSMATCHSTATS chain by the specified number of entries.
**</HTML>
**
**</AUTO>
**============================================================================
*/
static char mergeQABinCmd[] = "mergeQABin";
static int mergeQABinFlg = FTCL_ARGV_NO_LEFTOVERS;
static ftclArgvInfo mergeQABinTbl[] = {
   {NULL, FTCL_ARGV_HELP, NULL, NULL,
    "Bin a TSMATCHSTATS chain by the specified number of entries\n"},
   {"<chain>", FTCL_ARGV_STRING, NULL, NULL,
    "Handle to the TSMATCHSTATS chain"},
   {"<binSize>", FTCL_ARGV_INT, NULL, NULL,
    "number of entries per bin"},
   {NULL, FTCL_ARGV_END, NULL, NULL, NULL}
};

static int
tclMergeQABin (ClientData clientData,
			Tcl_Interp *interp,
			int argc,
			char **argv)
{
  CHAIN *inChain = NULL;
  char *inChainName = NULL;
  int binSize = 0;
  int status;

  // Parse the arguments
  mergeQABinTbl[1].dst = &inChainName;
  mergeQABinTbl[2].dst = &binSize;
  if ((status = shTclParseArgv(interp, &argc, argv, mergeQABinTbl,
			       mergeQABinFlg, mergeQABinCmd))
      != FTCL_ARGV_SUCCESS)
    return status;

  // Translate the handle to the CHAIN
  if (shTclAddrGetFromName (interp, inChainName, (void **)&inChain,
			    "CHAIN") != TCL_OK)
    {
      Tcl_AppendResult(interp, runMergeQACmd,
		       ": arg 'chain' not a CHAIN", NULL);
      return(TCL_ERROR);
    }

  // Check that we've got the right CHAIN type and a valid bin size
  if (shChainTypeGet(inChain) != shTypeGetFromName ("TSMATCHSTATS"))
    {
      Tcl_AppendResult(interp, runMergeQACmd,
		       ": arg 'chain' wrong CHAIN type", NULL);
      return(TCL_ERROR);
    }
  if (binSize < 2)
    {
      Tcl_AppendResult(interp, runMergeQACmd,
		       ": arg 'binSize' must be > 1", NULL);
      return(TCL_ERROR);
    }

  // Allocate output chain
  CHAIN *outChain = shChainNew("TSMATCHSTATS");

  // Loop through the input chain, binning as we go.
  int j = 0;
  int nElems = shChainSize(inChain);
  TSMATCHSTATS *bin = NULL;
  TSMATCHSTATS *elem = NULL;
  CURSOR_T crsr = shChainCursorNew(inChain);
  while ((elem = (TSMATCHSTATS*) shChainWalk(inChain, crsr, NEXT)) != NULL)
    {
      if (j % binSize == 0)
	{
	  // If new bin, allocate new structure and initialize all to zero.
	  // If not enough elements left-over for a new bin, then break out.
	  if (nElems - j < binSize) break;
	  bin = (TSMATCHSTATS*) shMalloc(sizeof(TSMATCHSTATS));
	  bin->field = 0;
	  tsMatchStatsInit(bin);
	}

      // Update counts.
      bin->field += elem->field;
      bin->nOverlap += elem->nOverlap;
      bin->nOverlapBright += elem->nOverlapBright;
      bin->nMismatches += elem->nMismatches;
      bin->nMismatchesBright += elem->nMismatchesBright;
      bin->nUnknown += elem->nMismatches;
      bin->nUnknownBright += elem->nMismatchesBright;
      bin->stars.nObjects += elem->stars.nObjects;
      bin->stars.nObjectsBright += elem->stars.nObjectsBright;
      bin->stars.nMatches += elem->stars.nMatches;
      bin->stars.nMatchesBright += elem->stars.nMatchesBright;
      bin->stars.row.nPairs += elem->stars.row.nPairs;
      bin->stars.row.av += elem->stars.row.av * elem->stars.row.nPairs;
      bin->stars.row.rms += (elem->stars.row.rms * elem->stars.row.rms +
			     elem->stars.row.av * elem->stars.row.av) *
	elem->stars.row.nPairs;
      bin->stars.col.nPairs += elem->stars.col.nPairs;
      bin->stars.col.av += elem->stars.col.av * elem->stars.col.nPairs;
      bin->stars.col.rms += (elem->stars.col.rms * elem->stars.col.rms +
			     elem->stars.col.av * elem->stars.col.av) *
	elem->stars.col.nPairs;
      for (int i = 0; i < nFilters; i++)
	{
	  bin->stars.mag[i].nPairs += elem->stars.mag[i].nPairs;
	  bin->stars.mag[i].av +=
	    elem->stars.mag[i].av * elem->stars.mag[i].nPairs;
	  bin->stars.mag[i].rms +=
	    (elem->stars.mag[i].rms * elem->stars.mag[i].rms +
	     elem->stars.mag[i].av * elem->stars.mag[i].av) *
	    elem->stars.mag[i].nPairs;
	}
      for (i = 0; i < nColors; i++)
	{
	  bin->stars.color[i].nPairs += elem->stars.color[i].nPairs;
	  bin->stars.color[i].av +=
	    elem->stars.color[i].av * elem->stars.color[i].nPairs;
	  bin->stars.color[i].rms +=
	    (elem->stars.color[i].rms * elem->stars.color[i].rms +
	     elem->stars.color[i].av * elem->stars.color[i].av) *
	    elem->stars.color[i].nPairs;
	}
      bin->galaxies.nObjects += elem->galaxies.nObjects;
      bin->galaxies.nObjectsBright += elem->galaxies.nObjectsBright;
      bin->galaxies.nMatches += elem->galaxies.nMatches;
      bin->galaxies.nMatchesBright += elem->galaxies.nMatchesBright;
      bin->galaxies.row.nPairs += elem->galaxies.row.nPairs;
      bin->galaxies.row.av +=
	elem->galaxies.row.av * elem->galaxies.row.nPairs;
      bin->galaxies.row.rms +=
	(elem->galaxies.row.rms * elem->galaxies.row.rms +
	 elem->galaxies.row.av * elem->galaxies.row.av) *
	elem->galaxies.row.nPairs;
      bin->galaxies.col.nPairs += elem->galaxies.col.nPairs;
      bin->galaxies.col.av +=
	elem->galaxies.col.av * elem->galaxies.col.nPairs;
      bin->galaxies.col.rms +=
	(elem->galaxies.col.rms * elem->galaxies.col.rms +
	 elem->galaxies.col.av * elem->galaxies.col.av) *
	elem->galaxies.col.nPairs;
      for (i = 0; i < nFilters; i++)
	{
	  bin->galaxies.mag[i].nPairs += elem->galaxies.mag[i].nPairs;
	  bin->galaxies.mag[i].av +=
	    elem->galaxies.mag[i].av * elem->galaxies.mag[i].nPairs;
	  bin->galaxies.mag[i].rms +=
	    (elem->galaxies.mag[i].rms * elem->galaxies.mag[i].rms +
	     elem->galaxies.mag[i].av * elem->galaxies.mag[i].av) *
	    elem->galaxies.mag[i].nPairs;
	}
      for (i = 0; i < nColors; i++)
	{
	  bin->galaxies.color[i].nPairs += elem->galaxies.color[i].nPairs;
	  bin->galaxies.color[i].av +=
	    elem->galaxies.color[i].av * elem->galaxies.color[i].nPairs;
	  bin->galaxies.color[i].rms +=
	    (elem->galaxies.color[i].rms * elem->galaxies.color[i].rms +
	     elem->galaxies.color[i].av * elem->galaxies.color[i].av) *
	    elem->galaxies.color[i].nPairs;
	}

      // If bin is finished, finish off stats
      j++;
      if (j % binSize == 0)
	{
	  bin->field /= binSize;

	  // Final star statistics
	  if (bin->stars.row.nPairs > 0)
	    {
	      bin->stars.row.av /= bin->stars.row.nPairs;
	      double meanSq = bin->stars.row.rms / bin->stars.row.nPairs -
		bin->stars.row.av * bin->stars.row.av;
	      bin->stars.row.rms = (meanSq > 0. ? sqrt(meanSq) : 0.);
	    }
	  if (bin->stars.col.nPairs > 0)
	    {
	      bin->stars.col.av /= bin->stars.col.nPairs;
	      double meanSq = bin->stars.col.rms / bin->stars.col.nPairs -
		bin->stars.col.av * bin->stars.col.av;
	      bin->stars.col.rms = (meanSq > 0. ? sqrt(meanSq) : 0.);
	    }
	  for (i = 0; i < nFilters; i++)
	    if (bin->stars.mag[i].nPairs > 0)
	      {
		bin->stars.mag[i].av /= bin->stars.mag[i].nPairs;
		double meanSq =
		  bin->stars.mag[i].rms / bin->stars.mag[i].nPairs -
		  bin->stars.mag[i].av * bin->stars.mag[i].av;
		bin->stars.mag[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.);
	      }
	  for (i = 0; i < nColors; i++)
	    if (bin->stars.color[i].nPairs > 0)
	      {
		bin->stars.color[i].av /= bin->stars.color[i].nPairs;
		double meanSq =
		  bin->stars.color[i].rms / bin->stars.color[i].nPairs -
		  bin->stars.color[i].av * bin->stars.color[i].av;
		bin->stars.color[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.);
	      }

	  // Final galaxy statistics
	  if (bin->galaxies.row.nPairs > 0)
	    {
	      bin->galaxies.row.av /= bin->galaxies.row.nPairs;
	      double meanSq =
		bin->galaxies.row.rms / bin->galaxies.row.nPairs -
		bin->galaxies.row.av * bin->galaxies.row.av;
	      bin->galaxies.row.rms = (meanSq > 0. ? sqrt(meanSq) : 0.);
	    }
	  if (bin->galaxies.col.nPairs > 0)
	    {
	      bin->galaxies.col.av /= bin->galaxies.col.nPairs;
	      double meanSq =
		bin->galaxies.col.rms / bin->galaxies.col.nPairs -
		bin->galaxies.col.av * bin->galaxies.col.av;
	      bin->galaxies.col.rms = (meanSq > 0. ? sqrt(meanSq) : 0.);
	    }
	  for (i = 0; i < nFilters; i++)
	    if (bin->galaxies.mag[i].nPairs > 0)
	      {
		bin->galaxies.mag[i].av /= bin->galaxies.mag[i].nPairs;
		double meanSq =
		  bin->galaxies.mag[i].rms / bin->galaxies.mag[i].nPairs -
		  bin->galaxies.mag[i].av * bin->galaxies.mag[i].av;
		bin->galaxies.mag[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.);
	      }
	  for (i = 0; i < nColors; i++)
	    if (bin->galaxies.color[i].nPairs > 0)
	      {
		bin->galaxies.color[i].av /= bin->galaxies.color[i].nPairs;
		double meanSq =
		  bin->galaxies.color[i].rms / bin->galaxies.color[i].nPairs -
		  bin->galaxies.color[i].av * bin->galaxies.color[i].av;
		bin->galaxies.color[i].rms = (meanSq > 0. ? sqrt(meanSq) : 0.);
	      }

	  // Add to chain
	  if (shChainElementAddByPos(outChain, bin, "TSMATCHSTATS", TAIL,
				     AFTER) != SH_SUCCESS)
	    {
	      shFree(bin);
	      taGenericChainDestroy(outChain);
	      Tcl_AppendResult(interp, runMergeQACmd,
			       ": couldn't append to output chain", NULL);
	      return TCL_ERROR;
	    }
	}
    }
  shChainCursorDel(inChain, crsr);

  // Return a handle to the binned statistics chain
  char vName[HANDLE_NAMELEN];
  if (shTclHandleNew(interp, vName, "CHAIN", (void *)outChain) != TCL_OK)
    {
      // Error.  Delete the chain and return an error.
      taGenericChainDestroy(outChain);
      Tcl_AppendResult(interp, runMergeQACmd,
		       ": error allocating output chain TCL handle", NULL);
      return TCL_ERROR;
    }
  Tcl_SetResult(interp, vName, TCL_VOLATILE);
  return TCL_OK;
}

//****************************************************************************/
//
// Declare my new tcl verbs to tcl
//
void
tsTclMergeDeclare(Tcl_Interp *interp) 
{
  char *tclHelpFacil = "ts";
  shTclDeclare(interp, setObjectStatusDontCallCmd,
	       (Tcl_CmdProc *)tclSetObjectStatusDontCall,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, setObjectStatusDontCallTbl,
			       setObjectStatusDontCallFlg,
			       setObjectStatusDontCallCmd),
	       shTclGetUsage(interp, setObjectStatusDontCallTbl,
			     setObjectStatusDontCallFlg,
			     setObjectStatusDontCallCmd));
  shTclDeclare(interp, runUnsetObjectStatusCmd,
	       (Tcl_CmdProc *)tclRunUnsetObjectStatus,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, runUnsetObjectStatusTbl,
			       runUnsetObjectStatusFlg,
			       runUnsetObjectStatusCmd),
	       shTclGetUsage(interp, runUnsetObjectStatusTbl,
			     runUnsetObjectStatusFlg,
			     runUnsetObjectStatusCmd));
  shTclDeclare(interp, selfMergeQACmd,
	       (Tcl_CmdProc *)tclSelfMergeQA,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, selfMergeQATbl, selfMergeQAFlg,
			       selfMergeQACmd),
	       shTclGetUsage(interp, selfMergeQATbl, selfMergeQAFlg,
			     selfMergeQACmd));
  shTclDeclare(interp, fieldPlotCmd,
	       (Tcl_CmdProc *)tclFieldPlot,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, fieldPlotTbl, fieldPlotFlg,
			       fieldPlotCmd),
	       shTclGetUsage(interp, fieldPlotTbl, fieldPlotFlg,
			     fieldPlotCmd));
  shTclDeclare(interp, runMergeDontCallCmd,
	       (Tcl_CmdProc *)tclRunMergeDontCall,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, runMergeDontCallTbl,
			       runMergeDontCallFlg, runMergeDontCallCmd),
	       shTclGetUsage(interp, runMergeDontCallTbl, runMergeDontCallFlg,
			     runMergeDontCallCmd));
  shTclDeclare(interp, runUnmergeCmd,
	       (Tcl_CmdProc *)tclRunUnmerge,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, runUnmergeTbl, runUnmergeFlg,
			       runUnmergeCmd),
	       shTclGetUsage(interp, runUnmergeTbl, runUnmergeFlg,
			     runUnmergeCmd));
  shTclDeclare(interp, runMergeQACmd,
	       (Tcl_CmdProc *)tclRunMergeQA,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, runMergeQATbl, runMergeQAFlg,
			       runMergeQACmd),
	       shTclGetUsage(interp, runMergeQATbl, runMergeQAFlg,
			     runMergeQACmd));
  shTclDeclare(interp, matchedPairChainDelCmd,
	       (Tcl_CmdProc *)tclMatchedPairChainDel,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, matchedPairChainDelTbl,
			       matchedPairChainDelFlg,
			       matchedPairChainDelCmd),
	       shTclGetUsage(interp, matchedPairChainDelTbl,
			     matchedPairChainDelFlg,
			     matchedPairChainDelCmd));
  shTclDeclare(interp, mergeQABinCmd,
	       (Tcl_CmdProc *)tclMergeQABin,
	       (ClientData) 0,	(Tcl_CmdDeleteProc *)NULL, tclHelpFacil,
	       shTclGetArgInfo(interp, mergeQABinTbl, mergeQABinFlg,
			       mergeQABinCmd),
	       shTclGetUsage(interp, mergeQABinTbl, mergeQABinFlg,
			     mergeQABinCmd));
}