Index: trunk/Ohana/src/relastro/Makefile
===================================================================
--- trunk/Ohana/src/relastro/Makefile	(revision 40378)
+++ trunk/Ohana/src/relastro/Makefile	(revision 40382)
@@ -52,4 +52,5 @@
 $(SRC)/UpdateChips.$(ARCH).o         \
 $(SRC)/UpdateStacks.$(ARCH).o         \
+$(SRC)/UpdateStacksWithFit.$(ARCH).o \
 $(SRC)/UpdateMosaic.$(ARCH).o        \
 $(SRC)/UpdateObjects.$(ARCH).o       \
Index: trunk/Ohana/src/relastro/include/relastro.h
===================================================================
--- trunk/Ohana/src/relastro/include/relastro.h	(revision 40378)
+++ trunk/Ohana/src/relastro/include/relastro.h	(revision 40382)
@@ -367,4 +367,5 @@
 int    USE_ICRF_POLE;
 
+int    FIT_STACKS;
 int    REPAIR_STACKS;
 int    USE_IMAGE_COORDS_FOR_REPAIR;
@@ -623,4 +624,5 @@
 int UpdateChips (Catalog *catalog, int Ncatalog, int Nloop);
 int UpdateStacks (Catalog *catalog, int Ncatalog);
+int UpdateStacksWithFit (Catalog *catalog, int Ncatalog);
 int UpdateMosaic (Catalog *catalog, int Ncatalog);
 int UpdateMeasures (Catalog *catalog, int Ncatalog);
Index: trunk/Ohana/src/relastro/src/UpdateStacksWithFit.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateStacksWithFit.c	(revision 40382)
+++ trunk/Ohana/src/relastro/src/UpdateStacksWithFit.c	(revision 40382)
@@ -0,0 +1,86 @@
+# include "relastro.h"
+
+// This stack analysis allows us to fit the chip model to each stack
+int UpdateStacksWithFit (Catalog *catalog, int Ncatalog) {
+
+  if (SKIP_PS1_STACK) return TRUE;
+
+  off_t Nimage;
+  Image *image = getimages (&Nimage, NULL);
+
+  off_t Nskip   = 0;
+  off_t NnewFit = 0;
+  off_t NoldFit = 0;
+
+  // measure the scatter for each stack
+  for (off_t i = 0; i < Nimage; i++) {
+
+    /* skip all except stack images */
+    if (!isGPC1stack(image[i].photcode)) continue;
+
+    /* convert measure coordinates to raw entries */
+    off_t Nraw;
+    StarData *raw = getImageRaw (catalog, Ncatalog, i, &Nraw, MODE_SIMPLE);
+    if (!raw) {
+      Nskip ++;
+      continue;
+    }
+    if (Nraw <= IMFIT_TOO_FEW) {
+      Nskip ++;
+      free (raw);
+      continue;
+    }
+
+    /* convert average coordinates to ref entries */
+    off_t Nref;
+    StarData *ref = getImageRef (catalog, Ncatalog, i, &Nref, MODE_SIMPLE);
+    if (!ref) {
+      Nskip ++;
+      free (raw);
+      continue;
+    }
+
+    // note that Nraw & Nref must be equal: if not, we made a programming error in one of these two functions.
+    assert (Nraw == Nref);
+
+    // save these in case of failure
+    Coords oldCoords;
+    SaveCoords (&oldCoords, &image[i].coords);
+
+    float dXpixSys = image[i].dXpixSys;
+    float dYpixSys = image[i].dYpixSys;
+    int   nFitAstr = image[i].nFitAstrom;
+
+    // FitChip does iterative, clipped fitting
+    if (!FitChip (raw, ref, Nraw, &image[i])) {
+      if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT", Nused %d of %d\n", image[i].name, i, Nraw, image[i].nFitAstrom, image[i].nstar);
+
+      // restore status quo ante (replace truMap with tmpMap)
+      RestoreCoords (&image[i].coords, &oldCoords, &image[i]);
+      image[i].dXpixSys = dXpixSys;
+      image[i].dYpixSys = dYpixSys; 
+      image[i].nFitAstrom = nFitAstr;
+
+      NoldFit ++;
+      free (raw);
+      free (ref);
+      continue;
+    }
+
+    AstromOffsetMapFree (oldCoords.offsetMap);
+
+    // Apply the modified coords back to the measure.R,D.  Note that raw.R,D, ref.L,M, etc
+    // are all automatically updated in this block because they are re-generated from
+    // image.coords on each pass.
+    setImageRaw (catalog, Ncatalog, i, raw, Nraw, MODE_MOSAIC);
+
+    NnewFit ++;
+    free (raw);
+    free (ref);
+  }
+
+  fprintf (stderr, "UpdateStacksWithFit: %d fitted, %d keep old, %d skipped\n", (int) NnewFit, (int) NoldFit, (int) Nskip);
+  return (TRUE);
+}
+
+// XXX 2 hardwired hacks in this file: 1) photcode hardwired for GPC1 stacks, 2) platescale
Index: trunk/Ohana/src/relastro/src/args.c
===================================================================
--- trunk/Ohana/src/relastro/src/args.c	(revision 40378)
+++ trunk/Ohana/src/relastro/src/args.c	(revision 40382)
@@ -861,4 +861,11 @@
     remove_argument (N, &argc, argv);
     RELASTRO_OP = OP_REPAIR_WARPS;
+  }
+
+  // the default is to only measure the scatter for the stacks, not to perform a new fit
+  FIT_STACKS = FALSE;
+  if ((N = get_argument (argc, argv, "-fit-stacks"))) {
+    remove_argument (N, &argc, argv);
+    FIT_STACKS = TRUE;
   }
 
Index: trunk/Ohana/src/relastro/src/relastro_images.c
===================================================================
--- trunk/Ohana/src/relastro/src/relastro_images.c	(revision 40378)
+++ trunk/Ohana/src/relastro/src/relastro_images.c	(revision 40382)
@@ -87,6 +87,11 @@
       }
 
-      // measure scatter for stacks
-      UpdateStacks (catalog, Ncatalog);
+      if (FIT_STACKS) {
+	// model astrometry for stacks
+	UpdateStacksWithFit (catalog, Ncatalog);
+      } else {
+	// measure scatter for stacks
+	UpdateStacks (catalog, Ncatalog);
+      }
       MARKTIME("UpdateStacks: %f sec\n", dtime);
 
Index: trunk/Ohana/src/relastro/src/relastro_parallel_images.c
===================================================================
--- trunk/Ohana/src/relastro/src/relastro_parallel_images.c	(revision 40378)
+++ trunk/Ohana/src/relastro/src/relastro_parallel_images.c	(revision 40382)
@@ -150,6 +150,11 @@
       }
 
-      // measure scatter for stacks
-      UpdateStacks (catalog, Ncatalog);
+      if (FIT_STACKS) {
+	// model astrometry for stacks
+	UpdateStacksWithFit (catalog, Ncatalog);
+      } else {
+	// measure scatter for stacks
+	UpdateStacks (catalog, Ncatalog);
+      }
       LOGRTIME("UpdateStacks on %s, host %d: %f sec\n", myHostName, REGION_HOST_ID, dtime);
 
