Index: trunk/Ohana/src/kapa2/include/prototypes.h
===================================================================
--- trunk/Ohana/src/kapa2/include/prototypes.h	(revision 32687)
+++ trunk/Ohana/src/kapa2/include/prototypes.h	(revision 32695)
@@ -126,4 +126,5 @@
 void 	      PSPixmap24 	  PROTO((Graphic *graphic, KapaImageWidget *image, FILE *f));
 void 	      PSPixmap32 	  PROTO((Graphic *graphic, KapaImageWidget *image, FILE *f));
+void          PSPixmap_3byte      PROTO((Graphic *graphic, KapaImageWidget *image, FILE *f));
 
 /* kapa bDraw Functions */
Index: trunk/Ohana/src/kapa2/include/structures.h
===================================================================
--- trunk/Ohana/src/kapa2/include/structures.h	(revision 32687)
+++ trunk/Ohana/src/kapa2/include/structures.h	(revision 32695)
@@ -26,6 +26,8 @@
   XFontStruct   *font;
   Cursor         cursor;
-  int            x,  y;
-  unsigned int   dx, dy;
+  int            x,  y;	      // corner coord in X world coords
+  unsigned int   dx, dy;      // size of window in X coords
+  int            xwin, ywin;  // corner coord of display subregion (eg, png, ps plot)
+  int            dxwin, dywin; // corner coord of display subregion (eg, png, ps plot)
   CCNode        *cube;
   XColor        *cmap;
Index: trunk/Ohana/src/kapa2/src/PSPixmap.c
===================================================================
--- trunk/Ohana/src/kapa2/src/PSPixmap.c	(revision 32687)
+++ trunk/Ohana/src/kapa2/src/PSPixmap.c	(revision 32695)
@@ -1,5 +1,3 @@
 # include "Ximage.h"
-
-// XXX this stuff has been broken by the conversion to the pixmap stuff
 
 void PSPixmap8 (Graphic *graphic, KapaImageWidget *image, FILE *f) {
@@ -146,40 +144,127 @@
 }
 
-# if (0)
-// XXX needs work!
-void PSPixmap32_RGB (Graphic *graphic, KapaImageWidget *image, FILE *f) {
-
-  int i, k, m, val;
-  double Nchar, Npix, start, slope, frac;
-  unsigned int *buff;
-  unsigned long back;
-  unsigned char 
-
-  ALLOCATE (pixelR, unsigned char, graphic[0].Npixels);
-  ALLOCATE (pixelG, unsigned char, graphic[0].Npixels);
-  ALLOCATE (pixelB, unsigned char, graphic[0].Npixels);
+# define WHITE_R 255
+# define WHITE_G 255
+# define WHITE_B 255
+
+void PSPixmap_3byte (Graphic *graphic, KapaImageWidget *image, FILE *f) {
+
+  int i, j, ii, jj;
+  int i_start, i_end, j_start, j_end;
+  int I_start, J_start;
+  int dropback;  /* this is a bit of a kludge... */
+  int dx, dy, DX, DY, inDX, inDY, Xs, Ys;
+  int expand_in, expand_out;
+  double expand, Ix, Iy;
+  unsigned short *in_pix, *in_pix_ref;
+  unsigned char *pixel1, *pixel2, *pixel3;
+
+  if (image == NULL) return;
+
+  ALLOCATE (pixel1, unsigned char, graphic[0].Npixels);
+  ALLOCATE (pixel2, unsigned char, graphic[0].Npixels);
+  ALLOCATE (pixel3, unsigned char, graphic[0].Npixels);
 
   /** cmap[i].pixel must be defined even if X is not used **/
   for (i = 0; i < graphic[0].Npixels; i++) { /* set up pixel array */
-    pixelR[i] = graphic[0].cmap[i].red >> 8;
-    pixelG[i] = graphic[0].cmap[i].green >> 8;
-    pixelB[i] = graphic[0].cmap[i].blue >> 8;
-  }
-
-  for (i = 0; i < image[0].picture.dy; i++) {
-    for (k = 0; k < image[0].picture.dx; k++, buff++) {
-      if (*buff == back) 
-	val = Nchar;
-      else {
-	for (m = 0; (graphic[0].cmap[m].pixel != *buff) && (m < Npix); m++);
-	val = Nchar - frac * MIN (MAX (start + m * slope, 0), Npix);
-      }
-      fprintf (f, "%02x", val);
-      if (!((k+1) % 40)) fprintf (f, "\n"); 
-    }
-    fprintf (f, "\n");
-    buff -= 2*image[0].picture.dx;
-  }
-  return;
-}
-# endif
+    pixel1[i] = graphic[0].cmap[i].red >> 8;
+    pixel2[i] = graphic[0].cmap[i].green >> 8;
+    pixel3[i] = graphic[0].cmap[i].blue >> 8;
+  }
+
+  assert ((image[0].picture.expand >= 1) || (image[0].picture.expand <= -2));
+  expand = expand_in = expand_out = 1.0;
+  if (image[0].picture.expand > 0) {
+    expand = 1 / (1.0*image[0].picture.expand);
+    expand_out = image[0].picture.expand;
+    expand_in  = 1;
+  }
+  if (image[0].picture.expand < 0) {
+    expand = fabs((double)image[0].picture.expand);
+    expand_out = 1;
+    expand_in  = -image[0].picture.expand;
+  }
+
+  Xs = image[0].picture.x;
+  Ys = image[0].picture.y;
+  dx = image[0].picture.dx;
+  dy = image[0].picture.dy;
+  DX = image[0].image[0].matrix.Naxis[0];
+  DY = image[0].image[0].matrix.Naxis[1];
+
+  // i_start, j_start are the closest lit screen pixel to 0,0
+  // I_start, J_start are the image pixel corresponding to i_start, j_start
+  Picture_Lower (&i_start, &j_start, &I_start, &J_start, &image[0].image[0].matrix, &image[0].picture);
+
+  // i_end, j_end are the closest lit screen pixel to dx, dy
+  // I_end, J_end are the image pixel corresponding to i_end, j_end
+  Picture_Upper (&i_end, &j_end, i_start, j_start, &image[0].image[0].matrix, &image[0].picture);
+
+  assert (i_start <= i_end);
+  assert (j_start <= j_end);
+
+  Ix = image[0].picture.flipx ? I_start - 1 : I_start;
+  Iy = image[0].picture.flipy ? J_start - 1 : J_start;
+
+  inDX = image[0].picture.flipx ? -1 : +1;
+  inDY = image[0].picture.flipy ? -1 : +1;
+
+  dropback = expand_out - (i_end - i_start) % expand_out;
+  if ((i_end - i_start) % expand_out == 0) dropback = 0;
+
+  in_pix_ref  = &image[0].pixmap[DX*(int)MAX(Iy,0) + (int)MAX(Ix,0)];
+
+  /********** below we do the mapping from buffer pixels (in) to picture pixels (out) **********/
+
+  // add in occasional return chars
+
+  /**** fill in bottom area ****/
+  for (j = 0; j < j_start; j++) {
+    for (i = 0; i < dx; i++) {
+      fprintf (f, "%02x%02x%02x", WHITE_R, WHITE_G, WHITE_B);
+    }
+  }
+  
+  // probably could do this all smarter with scale operations in PS...
+
+  /*** fill in the image data region ***/
+  for (j = j_start; j < j_end; j+= expand_out, in_pix_ref += inDY*expand_in*DX) {
+    
+    // repeat the section below 'expand_out' times
+    for (jj = 0; jj < expand_out; jj++) {
+
+      /* create one output image line */
+      in_pix = in_pix_ref;
+
+      /**** fill in area to the left of the picture ****/
+      for (i = 0; i < i_start; i++) {
+	fprintf (f, "%02x%02x%02x", WHITE_R, WHITE_G, WHITE_B);
+      }
+    
+      /*** fill in the picture region ***/
+      for (i = i_start; i < i_end; i+=expand_out, in_pix += inDX*expand_in) {
+	for (ii = 0; ii < expand_out; ii++) {
+	  fprintf (f, "%02x%02x%02x", pixel1[*in_pix], pixel2[*in_pix], pixel3[*in_pix]);
+	}
+      }
+    
+      /**** fill in area to the right of the picture ****/
+      for (i = i_end; i < dx; i++) {
+	fprintf (f, "%02x%02x%02x", WHITE_R, WHITE_G, WHITE_B);
+      }
+    }
+  }
+
+  /**** fill in top area ****/
+  for (j = j_end; j < dy; j++) {
+    for (i = 0; i < dx; i++) { 
+      fprintf (f, "%02x%02x%02x", WHITE_R, WHITE_G, WHITE_B);
+    }
+  }
+
+  free (pixel1);
+  free (pixel2);
+  free (pixel3);
+
+  return;
+}
Index: trunk/Ohana/src/kapa2/src/PSimage.c
===================================================================
--- trunk/Ohana/src/kapa2/src/PSimage.c	(revision 32687)
+++ trunk/Ohana/src/kapa2/src/PSimage.c	(revision 32695)
@@ -13,4 +13,11 @@
   graphic = GetGraphic();
 
+  // Update this to generate color PS images: I need to test if the image if color or BW (based
+  // on the colormap).  If it is a color image, I need to generate a 24 bit image. the header
+  // should look the same (Nx Ny 8 [1 0 0 1 0 0]) but then it should finish with "false 3
+  // colorimage" instead of just "image" (see ../doc/color.image.s for an example).  the hex
+  // string needs to have RGB represented as chars (2 hex chars per R,G,B). 
+  // the values come from the following map: image[0].pixmap[i] -> cmap[pixel].red -> R (0-255)
+
   fprintf (f, " newpath %d %d moveto %d %d lineto %d %d lineto %d %d lineto closepath\n\n", 
 	   (int) image[0].picture.x,                       graphic->dy - (int) image[0].picture.y, 
@@ -21,23 +28,18 @@
   fprintf (f, "%d %d translate\n", (int) image[0].picture.x, graphic->dy - (int) image[0].picture.y - image[0].picture.dy);
   fprintf (f, "%d %d 8\n", image[0].picture.dx, image[0].picture.dy);
-  fprintf (f, "[1 0 0 1 0 0]\n");
+  fprintf (f, "[1 0 0 -1 0 %d]\n", image[0].picture.dy);
+  // write out the image in normal order, but flip in PS
+
+# if (0)
+
   fprintf (f, "{currentfile %d string readhexstring pop} image\n\n", image[0].picture.dx);
+  PSPixmap_1byte (graphic, image, f);
 
-  /******** First we draw the picture itself ********/
-  /* in !USE_XWINDOW, we'll have to change this to use the JPEG function */
-  switch (graphic[0].Nbits) {
-  case 8:
-    PSPixmap8 (graphic, image, f);
-    break;
-  case 16:
-    PSPixmap16 (graphic, image, f);
-    break;
-  case 24:
-    PSPixmap24 (graphic, image, f);
-    break;
-  case 32:
-    PSPixmap32 (graphic, image, f);
-    break;
-  }
+# else
+
+  fprintf (f, "{currentfile %d string readhexstring pop} false 3 colorimage\n\n", 3*image[0].picture.dx);
+  PSPixmap_3byte (graphic, image, f);
+
+# endif
 
   fprintf (f, "grestore %% end of image\n");
Index: trunk/Ohana/src/kapa2/src/Sections.c
===================================================================
--- trunk/Ohana/src/kapa2/src/Sections.c	(revision 32687)
+++ trunk/Ohana/src/kapa2/src/Sections.c	(revision 32695)
@@ -237,2 +237,35 @@
   return (TRUE);
 }
+
+int SectionMinBoundary (Graphic *graphic) {
+
+    int i;
+    int Xs = graphic->dx;
+    int Ys = graphic->dy;
+    int Xe = 0;
+    int Ye = 0;
+
+    // the boundary for a single section should probably be adjusted depending on the 
+    // image status (should not include the imtool portion)
+    for (i = 0; i < Nsections; i++) {
+	Xs = MIN (Xs, sections[i][0].x);
+	Ys = MIN (Ys, sections[i][0].y);
+	Xe = MAX (Xe, sections[i][0].x + sections[i][0].dx);
+	Ye = MAX (Ye, sections[i][0].y + sections[i][0].dy);
+    }
+
+    if ((Xs >= Xe) || (Ys >- Ye)) {
+	// default values for the region window
+	graphic->xwin  = 0;
+	graphic->ywin  = 0;
+	graphic->dxwin = graphic->dx;
+	graphic->dywin = graphic->dy; 
+    }	
+    
+    // set min/max boundary (min window containing max range of sections)
+    graphic->xwin = Xs;
+    graphic->ywin = Ys;
+    graphic->dxwin = Xe - Xs;
+    graphic->dywin = Ye - Ys;
+    return TRUE;
+}
Index: trunk/Ohana/src/kapa2/src/SetUpGraphic.c
===================================================================
--- trunk/Ohana/src/kapa2/src/SetUpGraphic.c	(revision 32687)
+++ trunk/Ohana/src/kapa2/src/SetUpGraphic.c	(revision 32695)
@@ -16,4 +16,10 @@
   graphic->dx = 512;
   graphic->dy = 512; 
+
+  // default values for the region window
+  graphic->xwin  = 0;
+  graphic->ywin  = 0;
+  graphic->dxwin = graphic->dx;
+  graphic->dywin = graphic->dy; 
 
   if (!USE_XWINDOW) {
Index: trunk/Ohana/src/kapa2/src/bDrawImage.c
===================================================================
--- trunk/Ohana/src/kapa2/src/bDrawImage.c	(revision 32687)
+++ trunk/Ohana/src/kapa2/src/bDrawImage.c	(revision 32695)
@@ -47,6 +47,7 @@
   }
 
-  Xs = image[0].picture.x;
-  Ys = image[0].picture.y;
+  // Xs,Ys are in full-frame coords.  can we trim here?
+  Xs = image[0].picture.x - graphic[0].xwin;
+  Ys = image[0].picture.y - graphic[0].ywin;
   dx = image[0].picture.dx;
   dy = image[0].picture.dy;
@@ -54,4 +55,5 @@
   DY = image[0].image[0].matrix.Naxis[1];
 
+  // the created buffer is supposed to contain the output windows
   if (buffer[0].Nx < Xs + dx) {
     fprintf (stderr, "invalid condition\n");
Index: trunk/Ohana/src/kapa2/src/bDrawIt.c
===================================================================
--- trunk/Ohana/src/kapa2/src/bDrawIt.c	(revision 32687)
+++ trunk/Ohana/src/kapa2/src/bDrawIt.c	(revision 32695)
@@ -13,5 +13,9 @@
   black = KapaColorByName ("black");
 
-  buffer = bDrawBufferCreate (graphic->dx, graphic->dy, Nbyte, palette, Npalette);
+  // SectionMinBoundary (graphic);
+
+  // if we want to trim, we'll need to carry about the start in graphic coords and
+  // the dx,dy size.  
+  buffer = bDrawBufferCreate (graphic->dxwin, graphic->dywin, Nbyte, palette, Npalette);
   bDrawSetStyle (buffer, black, 0, 0);
   
Index: trunk/Ohana/src/libdvo/src/dbExtractMeasures.c
===================================================================
--- trunk/Ohana/src/libdvo/src/dbExtractMeasures.c	(revision 32687)
+++ trunk/Ohana/src/libdvo/src/dbExtractMeasures.c	(revision 32695)
@@ -245,8 +245,8 @@
       break;
     case MEAS_PAR: /* OK */
-      value.Flt = average[0].R;
+      value.Flt = average[0].P;
       break;
     case MEAS_PAR_ERR: /* OK */
-      value.Flt = average[0].D;
+      value.Flt = average[0].dP;
       break;
     case MEAS_CHISQ_POS: /* OK */
Index: trunk/Ohana/src/libkapa/src/PSRotFont.c
===================================================================
--- trunk/Ohana/src/libkapa/src/PSRotFont.c	(revision 32687)
+++ trunk/Ohana/src/libkapa/src/PSRotFont.c	(revision 32695)
@@ -32,4 +32,24 @@
   dY = currentfont[65].ascent;
   
+  /***
+
+      update this to use the following PS code:
+
+      /ceshow { % (string) fontsize fontname x y 
+      gsave 
+      moveto findfont exch scalefont setfont % s 
+      gsave 
+      dup false charpath flattenpath pathbbox % s x0 y0 x1 y1 
+      grestore 
+      3 -1 roll sub % s x0 x1 dy 
+      3 1 roll sub % s dy -dx 
+      2 div exch % s -dx/2 dy 
+      -2 div % s -dx/2 -dy/2 
+      rmoveto show 
+      grestore 
+      } bind def 
+
+   ***/
+
   /* apply appropriate offset */
   Xoff = Yoff = 0;
Index: trunk/Ohana/src/relastro/include/relastro.h
===================================================================
--- trunk/Ohana/src/relastro/include/relastro.h	(revision 32687)
+++ trunk/Ohana/src/relastro/include/relastro.h	(revision 32695)
@@ -317,7 +317,7 @@
 int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon);
 int ParFactor (double *pR, double *pD, double R, double D, double T, double Tmean);
-int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts);
+int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int XVERB);
 int FitPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *pR, double *pD, int Npts);
-int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts);
+int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int XVERB);
 
 Mosaic *getMosaicForImage (off_t N);
Index: trunk/Ohana/src/relastro/src/FitPM.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPM.c	(revision 32687)
+++ trunk/Ohana/src/relastro/src/FitPM.c	(revision 32695)
@@ -2,5 +2,5 @@
 
 /* do we want an init function which does the alloc and a clear function to free? */
-int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts) {
+int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int XVERB) {
 
   int i;
@@ -77,4 +77,5 @@
     chisq += SQ(X[i] - Xf) / SQ(dX[i]);
     chisq += SQ(Y[i] - Yf) / SQ(dY[i]);
+    if (XVERB) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq);
   }
   fit[0].Nfit = Npts;
Index: trunk/Ohana/src/relastro/src/FitPMandPar.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPMandPar.c	(revision 32687)
+++ trunk/Ohana/src/relastro/src/FitPMandPar.c	(revision 32695)
@@ -2,5 +2,5 @@
 
 /* do we want an init function which does the alloc and a clear function to free? */
-int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts) {
+int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int XVERB) {
 
   int i;
@@ -99,8 +99,10 @@
   chisq = 0.0;
   for (i = 0; i < Npts; i++) {
-    Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].dp*pR[i];
-    Yf = fit[0].Do + fit[0].uD*T[i] + fit[0].dp*pD[i];
+    Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].p*pR[i];
+    Yf = fit[0].Do + fit[0].uD*T[i] + fit[0].p*pD[i];
     chisq += SQ(X[i] - Xf) / SQ(dX[i]);
     chisq += SQ(Y[i] - Yf) / SQ(dY[i]);
+    if (XVERB) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq);
+
   }
   fit[0].Nfit = Npts;
Index: trunk/Ohana/src/relastro/src/ImageOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/ImageOps.c	(revision 32687)
+++ trunk/Ohana/src/relastro/src/ImageOps.c	(revision 32695)
@@ -913,4 +913,5 @@
   if (!finite(measure[0].dR) || !finite(measure[0].dD)) return FALSE;
   if (!finite(measure[0].M)) return FALSE; //XXX is this necessary for all relastro tasks?
+  if (!finite(measure[0].dM)) return FALSE; //XXX is this necessary for all relastro tasks?
   
   /* select measurements by photcode, or equiv photcode, if specified */
Index: trunk/Ohana/src/relastro/src/UpdateObjects.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 32687)
+++ trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 32695)
@@ -84,4 +84,6 @@
       /* calculate the average value of R,D for a single star */
 
+      XVERB = FALSE;
+
       // skip objects which are known to be problematic
       // XXX include this code or not?
@@ -196,5 +198,8 @@
       }
       
-      XVERB = (catalog[i].measure[m].dM < 0.01) && (N == 6) && (mode == FIT_PM_ONLY);
+      // XVERB |= (catalog[i].averge[j].objID == 0xc90) && (catalog[i].average[j].catID == 0x2a1e);
+      XVERB |= (catalog[i].average[j].objID == OBJ_ID_SRC) && (catalog[i].average[j].catID == CAT_ID_SRC);
+      XVERB |= (catalog[i].average[j].objID == OBJ_ID_DST) && (catalog[i].average[j].catID == CAT_ID_DST);
+      // XVERB = (catalog[i].measure[m].dM < 0.01) && (N == 6) && (mode == FIT_PM_ONLY);
 
       // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
@@ -211,7 +216,7 @@
 	}	  
 
-	FitPM (&fitPM, X, dX, Y, dY, T, N);
-
-	if (XVERB) fprintf (stderr, "fitted:  %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq);
+	FitPM (&fitPM, X, dX, Y, dY, T, N, XVERB);
+
+	if (XVERB) fprintf (stderr, "fitted PM:  %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq);
 
 	// project Ro, Do back to RA,DEC
@@ -230,5 +235,5 @@
 	  ParFactor (&pX[k], &pY[k], R[k], D[k], T[k], Tmean);
 	}
-	FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N);
+	FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB);
 	XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords);
 	catalog[i].average[j].flags |= ID_STAR_FIT_PAR;
@@ -275,5 +280,5 @@
 	break;
       }
-      if (XVERB) fprintf (stderr, "A%f %f -> %f %f (%f,%f) pm=(%f %f)\n",
+      if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f)\n",
 			  catalog[i].average[j].R, 
 			  catalog[i].average[j].D, 
Index: trunk/Ohana/src/relastro/src/args.c
===================================================================
--- trunk/Ohana/src/relastro/src/args.c	(revision 32687)
+++ trunk/Ohana/src/relastro/src/args.c	(revision 32695)
@@ -14,4 +14,7 @@
   FIT_MODE = FIT_AVERAGE;
 
+  OBJ_ID_SRC = OBJ_ID_DST = 0;
+  CAT_ID_SRC = CAT_ID_DST = 0;
+
   if ((N = get_argument (argc, argv, "-merge-source"))) {
     if (N > argc - 6) usage_merge_source();
@@ -40,4 +43,26 @@
     remove_argument (N, &argc, argv);
     FIT_TARGET = TARGET_OBJECTS;
+  }
+
+  if ((N = get_argument (argc, argv, "-testobj1"))) {
+    if (N > argc - 3) usage ();
+    remove_argument (N, &argc, argv);
+    OBJ_ID_SRC = strtol(argv[N], &endptr, 0);
+    if (*endptr) usage ();
+    remove_argument (N, &argc, argv);
+    CAT_ID_SRC = strtol(argv[N], &endptr, 0);
+    if (*endptr) usage (); 
+    remove_argument (N, &argc, argv);
+  }
+
+  if ((N = get_argument (argc, argv, "-testobj2"))) {
+    if (N > argc - 3) usage ();
+    remove_argument (N, &argc, argv);
+    OBJ_ID_DST = strtol(argv[N], &endptr, 0);
+    if (*endptr) usage ();
+    remove_argument (N, &argc, argv);
+    CAT_ID_DST = strtol(argv[N], &endptr, 0);
+    if (*endptr) usage (); 
+    remove_argument (N, &argc, argv);
   }
 
Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 32687)
+++ trunk/psModules/src/imcombine/pmSubtraction.c	(revision 32695)
@@ -1487,4 +1487,5 @@
             }
         }
+	psFree(out1->covariance);
         out1->covariance = psImageCovarianceAverage(covars);
         psFree(covars);
@@ -1527,4 +1528,5 @@
             }
         }
+	psFree(out2->covariance);
         out2->covariance = psImageCovarianceAverage(covars);
         psFree(covars);
Index: trunk/psModules/src/objects/pmSource.c
===================================================================
--- trunk/psModules/src/objects/pmSource.c	(revision 32687)
+++ trunk/psModules/src/objects/pmSource.c	(revision 32695)
@@ -121,4 +121,5 @@
     source->modelFits = NULL;
     source->extFitPars = NULL;
+
     source->type = PM_SOURCE_TYPE_UNKNOWN;
     source->mode = PM_SOURCE_MODE_DEFAULT;
@@ -138,4 +139,5 @@
     source->apFluxErr 	     = NAN; 
 
+    source->windowRadius     = NAN;
     source->skyRadius  	     = NAN;
     source->skyFlux    	     = NAN;
@@ -203,8 +205,10 @@
 
     source->imageID          = in->imageID;
+
     source->type     	     = in->type;
     source->mode     	     = in->mode;
     source->mode2    	     = in->mode2;
     source->tmpFlags 	     = in->tmpFlags;
+
     source->psfMag     	     = in->psfMag;
     source->psfMagErr 	     = in->psfMagErr;
@@ -217,6 +221,13 @@
     source->apFlux    	     = in->apFlux;
     source->apFluxErr 	     = in->apFluxErr;
+
+    source->windowRadius     = in->windowRadius;
+    source->skyRadius  	     = in->skyRadius;  	
+    source->skyFlux    	     = in->skyFlux;    	
+    source->skySlope   	     = in->skySlope;   	
+
     source->pixWeightNotBad  = in->pixWeightNotBad;
     source->pixWeightNotPoor = in->pixWeightNotPoor;
+
     source->psfChisq         = in->psfChisq;
     source->crNsigma         = in->crNsigma;
@@ -259,4 +270,5 @@
     }
     mySource->region   = srcRegion;
+    mySource->windowRadius = Radius;
 
     return true;
@@ -328,4 +340,5 @@
         mySource->psfImage = NULL;
     }
+    mySource->windowRadius = Radius;
     return extend;
 }
Index: trunk/psModules/src/objects/pmSource.h
===================================================================
--- trunk/psModules/src/objects/pmSource.h	(revision 32687)
+++ trunk/psModules/src/objects/pmSource.h	(revision 32695)
@@ -97,4 +97,5 @@
     float apFluxErr;                    ///< apFluxErr corresponding to psfMag or extMag (depending on type)
 
+    float windowRadius;			///< size of box used for full analysis
     float skyRadius;			///< radius at which profile hits local sky (or goes flat)
     float skyFlux;			///< mean flux per pixel in aperture at which profile hits local sky (or goes flat)
Index: trunk/psphot/src/psphot.h
===================================================================
--- trunk/psphot/src/psphot.h	(revision 32687)
+++ trunk/psphot/src/psphot.h	(revision 32695)
@@ -108,6 +108,6 @@
 bool            psphotBlendFit_Threaded (psThreadJob *job);
 
-bool            psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule);
-bool            psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe);
+bool            psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState);
+bool            psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool ignoreState);
 
 bool            psphotAddNoise (pmConfig *config, const pmFPAview *view, const char *filerule);
@@ -176,5 +176,7 @@
 
 // in psphotReplaceUnfit.c:
-bool            psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe);
+bool            psphotRemoveAllSourcesByArray (const psArray *sources, const psMetadata *recipe);
+bool            psphotRemoveAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState);
+bool            psphotRemoveAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState);
 bool            psphotReplaceUnfitSources (psArray *sources, psImageMaskType maskVal);
 
@@ -453,4 +455,10 @@
 psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc);
 
+bool psphotSourceParents (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc);
+bool psphotSourceParentsReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index);
+
+bool psphotCopyPeaks (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc);
+bool psphotCopyPeaksReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index);
+
 bool psphotSersicModelClassGuessPCM (pmPCMdata *pcm, pmSource *source);
 void psphotSersicModelClassInit ();
Index: trunk/psphot/src/psphotBlendFit.c
===================================================================
--- trunk/psphot/src/psphotBlendFit.c	(revision 32687)
+++ trunk/psphot/src/psphotBlendFit.c	(revision 32695)
@@ -236,13 +236,15 @@
         pmSource *source = sources->data[i];
 
-# define TEST_X -420.0
-# define TEST_Y 300.0
+# if (0)
+# define TEST_X 34
+# define TEST_Y 28
    
 	if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
-	    fprintf (stderr, "test galaxy\n");
+	    fprintf (stderr, "test object\n");
 	}
 
 # undef TEST_X
 # undef TEST_Y
+# endif
 
         // skip non-astronomical objects (very likely defects)
Index: trunk/psphot/src/psphotEfficiency.c
===================================================================
--- trunk/psphot/src/psphotEfficiency.c	(revision 32687)
+++ trunk/psphot/src/psphotEfficiency.c	(revision 32695)
@@ -255,5 +255,5 @@
 
     // remove all sources, adding noise for subtracted sources
-    psphotRemoveAllSources(realSources, recipe);
+    psphotRemoveAllSourcesByArray(realSources, recipe);
 
 #if TESTING
Index: trunk/psphot/src/psphotFake.c
===================================================================
--- trunk/psphot/src/psphotFake.c	(revision 32687)
+++ trunk/psphot/src/psphotFake.c	(revision 32695)
@@ -166,5 +166,5 @@
 
     // remove all sources, adding noise for subtracted sources
-    psphotRemoveAllSources(realSources, recipe);
+    psphotRemoveAllSourcesByArray(realSources, recipe);
     psphotAddNoise(readout, realSources, recipe);
 
Index: trunk/psphot/src/psphotFitSourcesLinear.c
===================================================================
--- trunk/psphot/src/psphotFitSourcesLinear.c	(revision 32687)
+++ trunk/psphot/src/psphotFitSourcesLinear.c	(revision 32695)
@@ -342,9 +342,8 @@
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "solve matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
 
-    // XXXX **** philosophical question:
-    // we measure bright objects in three passes: 1) linear fit; 2) non-linear fit; 3) linear fit:
-    // should retain the chisq and errors from the intermediate non-linear fit?
-    // the non-linear fit provides better values for the position errors, and for
-    // extended sources, the shape errors
+    // Philosophical question: we measure bright objects in three passes: 1) linear fit; 2)
+    // non-linear fit; 3) linear fit: should we retain the chisq and errors from the
+    // intermediate non-linear fit?  the non-linear fit provides better values for the position
+    // errors, and for extended sources, the shape errors
 
     // adjust I0 for fitSources and subtract
Index: trunk/psphot/src/psphotMergeSources.c
===================================================================
--- trunk/psphot/src/psphotMergeSources.c	(revision 32687)
+++ trunk/psphot/src/psphotMergeSources.c	(revision 32695)
@@ -520,4 +520,176 @@
 }
 
+// copy the newPeaks from the detections of one pmFPAfile to another
+bool psphotCopyPeaks (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc)
+{
+    bool status = true;
+
+    int num = psphotFileruleCount(config, ruleSrc);
+
+    // skip the chisq image because it is a duplicate of the detection version
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+	if (i == chisqNum) continue; // skip chisq image
+        if (!psphotCopyPeaksReadout (config, view, ruleOut, ruleSrc, i)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to copy peaks from %s to %s entry %d", ruleSrc, ruleOut, i);
+            return false;
+        }
+    }
+    return true;
+}
+
+// add newly detected peaks to the existing list of sources
+bool psphotCopyPeaksReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index) {
+
+    bool status;
+
+    // find the currently selected readout
+    pmFPAfile *fileSrc = pmFPAfileSelectSingle(config->files, ruleSrc, index); // File of interest
+    psAssert (fileSrc, "missing file?");
+
+    pmReadout *readoutSrc = pmFPAviewThisReadout(view, fileSrc->fpa);
+    psAssert (readoutSrc, "missing readout?");
+
+    pmDetections *detections = psMetadataLookupPtr (&status, readoutSrc->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    // find the currently selected readout
+    pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, ruleOut, index); // File of interest
+    psAssert (fileOut, "missing file?");
+
+    pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa);
+    psAssert (readoutOut, "missing readout?");
+
+    // generate a new detection structure for the output filerule
+    pmDetections *detectionsOut = psMetadataLookupPtr (&status, readoutOut->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detectionsOut, "missing PSPHOT.DETECTIONS?");
+
+    psAssert (detectionsOut->peaks, "programming error");
+    psAssert (!detectionsOut->oldPeaks, "programming error");
+    psAssert (detections->peaks, "programming error");
+
+    // save the OUT existing peaks on oldPeaks
+    detectionsOut->oldPeaks = detectionsOut->peaks;
+    detectionsOut->peaks = psArrayAllocEmpty(detections->peaks->n);
+
+    for (int i = 0; i < detections->peaks->n; i++) {
+	psAssert (detections->peaks->data[i], "programming error");
+	pmPeak *peak = pmPeakAlloc (0, 0, 0.0, PM_PEAK_LONE);
+	pmPeakCopy(peak, detections->peaks->data[i]);
+	psArrayAdd (detectionsOut->peaks, 100, peak);
+	psFree (peak);
+    }
+    return true;
+}
+
+// create source parents children from ruleSrc for ruleOut for orphans
+bool psphotSourceParents (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc)
+{
+    bool status = true;
+
+    int num = psphotFileruleCount(config, ruleSrc);
+
+    // skip the chisq image because it is a duplicate of the detection version
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+	if (i == chisqNum) continue; // skip chisq image
+        if (!psphotSourceParentsReadout (config, view, ruleOut, ruleSrc, i)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to copy sources from %s to %s entry %d", ruleSrc, ruleOut, i);
+            return false;
+        }
+    }
+    return true;
+}
+
+// create source parents from ruleSrc for ruleOut for orphaned children for this readout.  
+bool psphotSourceParentsReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index) {
+
+    bool status;
+    int nParents = 0;
+    int nNonOrphans = 0;
+
+    // find the currently selected readout
+    pmFPAfile *fileSrc = pmFPAfileSelectSingle(config->files, ruleSrc, index); // File of interest
+    psAssert (fileSrc, "missing file?");
+
+    pmReadout *readoutSrc = pmFPAviewThisReadout(view, fileSrc->fpa);
+    psAssert (readoutSrc, "missing readout?");
+
+    pmDetections *detectionsSrc = psMetadataLookupPtr (&status, readoutSrc->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detectionsSrc, "missing detections?");
+
+    psArray *sourcesSrc = detectionsSrc->allSources;
+    psAssert (sourcesSrc, "missing sources?");
+
+    // find the currently selected readout
+    pmFPAfile *fileOut = pmFPAfileSelectSingle(config->files, ruleOut, index); // File of interest
+    psAssert (fileOut, "missing file?");
+
+    pmReadout *readoutOut = pmFPAviewThisReadout(view, fileOut->fpa);
+    psAssert (readoutOut, "missing readout?");
+
+    // generate a new detection structure for the output filerule
+    pmDetections *detectionsOut = psMetadataLookupPtr (&status, readoutOut->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detectionsOut, "missing PSPHOT.DETECTIONS?");
+
+    // loop over the sources, redefine their pixels to point at the new filerule image,
+    // copy the source data, and add a reference back to the original source
+    
+    // copy the sources from sourceSrcs to the new detection structure
+    for (int i = 0; i < sourcesSrc->n; i++) {
+      pmSource *sourceSrc = sourcesSrc->data[i];
+      if (sourceSrc->parent) {
+	  nNonOrphans ++;
+	  continue; // Not an orphan
+      }
+
+      pmSource *sourceOut = pmSourceCopy(sourceSrc);
+      sourceOut->parent = sourceSrc;
+      
+      // keep the original source flags
+      sourceOut->seq      = sourceSrc->seq;
+      sourceOut->type     = sourceSrc->type;
+      sourceOut->mode     = sourceSrc->mode;
+      sourceOut->mode2    = sourceSrc->mode2;
+      sourceOut->tmpFlags = sourceSrc->tmpFlags;
+
+      // does this copy all model data? (NO)
+      sourceOut->modelPSF = pmModelCopy(sourceSrc->modelPSF);
+      sourceOut->modelEXT = pmModelCopy(sourceSrc->modelEXT);
+
+      if (sourceSrc->modelFits) {
+	  sourceOut->modelFits = psArrayAlloc(sourceSrc->modelFits->n);
+	  for (int j = 0; j < sourceSrc->modelFits->n; j++) {
+	      sourceOut->modelFits->data[j] = pmModelCopy(sourceSrc->modelFits->data[j]);
+	  }
+      }
+
+      // drop the references to the original image pixels:
+      pmSourceFreePixels (sourceOut);
+
+      // allocate image, weight, mask for the new image for each peak
+      if (sourceOut->modelPSF) {
+	pmSourceRedefinePixels (sourceOut, readoutOut, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius);
+      }
+
+      // child sources have not been subtracted in this image, but this flag may be raised if
+      // they were subtracted in the parent's image
+      sourceOut->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+
+      nParents ++;
+      psArrayAdd (detectionsOut->allSources, 100, sourceOut);
+      psFree (sourceOut);
+    }
+    psLogMsg ("psphot", 3, "%d parents created, %d unorphaned children, %ld input vs %ld output", nParents, nNonOrphans, sourcesSrc->n, detectionsOut->allSources->n);
+
+    return true;
+}
+
 // create source children from ruleSrc for ruleOut
 bool psphotSourceChildren (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc)
@@ -542,7 +714,7 @@
 }
 
-// create source children from ruleSrc for ruleOut for this entry.  XXX currently, this is only
+// Create source children from ruleSrc for ruleOut for this entry.  Currently, this is only
 // used by psphotStackReadout (sources go on allSources so that psphotChoosePSF can be called
-// repeatedly)
+// repeatedly).
 bool psphotSourceChildrenReadout (pmConfig *config, const pmFPAview *view, const char *ruleOut, const char *ruleSrc, int index) {
 
@@ -569,15 +741,10 @@
     psAssert (readoutOut, "missing readout?");
 
-    // generate a new detection structure for the output filerule
-    pmDetections *detectionsOut = psMetadataLookupPtr (&status, readoutOut->analysis, "PSPHOT.DETECTIONS");
-    if (!detectionsOut) {
-	detectionsOut = pmDetectionsAlloc();
-	detectionsOut->allSources = psArrayAllocEmpty (100);
-	// save detections on the readout->analysis
-	if (!psMetadataAddPtr (readoutOut->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detectionsOut)) {
-	    psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout");
-	    return false;
-	}
-	psFree(detectionsOut); // a copy remains on the analysis metadata
+    // replace any existing DETECTION container on readoutOut->analysis with the new one
+    pmDetections *detectionsOut = pmDetectionsAlloc();
+    detectionsOut->allSources = psArrayAllocEmpty (100);
+    if (!psMetadataAddPtr (readoutOut->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detectionsOut)) {
+	psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout");
+	return false;
     }
 
@@ -593,4 +760,5 @@
       
       // keep the original source flags
+      sourceOut->seq      = sourceSrc->seq;
       sourceOut->type     = sourceSrc->type;
       sourceOut->mode     = sourceSrc->mode;
@@ -612,8 +780,8 @@
       pmSourceFreePixels (sourceOut);
 
+      // XXX do we need to skip the Chisq image sources?
+
       // allocate image, weight, mask for the new image for each peak
-      if (sourceOut->modelPSF) {
-	pmSourceRedefinePixels (sourceOut, readoutOut, sourceOut->peak->x, sourceOut->peak->y, sourceOut->modelPSF->fitRadius);
-      }
+      pmSourceRedefinePixels (sourceOut, readoutOut, sourceOut->peak->x, sourceOut->peak->y, sourceOut->windowRadius);
 
       // child sources have not been subtracted in this image, but this flag may be raised if
@@ -624,13 +792,14 @@
       psFree (sourceOut);
     }
-    psLogMsg ("psphot", 3, "%ld known sources supplied", detectionsOut->allSources->n);
-
-
-    return true;
-}
-
-// create source children associated with 'filerule' from the objectsSrc.  XXX currently, this
-// is only used by psphotStackReadout (sources go on allSources so that psphotChoosePSF can be
-// called repeatedly)
+    psLogMsg ("psphot", 3, "created %ld children", detectionsOut->allSources->n);
+
+    psFree(detectionsOut); // a copy remains on the analysis metadata
+
+    return true;
+}
+
+// create source children associated with 'filerule' from the objectsSrc.  returns a new object
+// array containing the child sources.  XXX currently, this is only used by psphotStackReadout
+// (sources go on allSources so that psphotChoosePSF can be called repeatedly)
 psArray *psphotSourceChildrenByObject (pmConfig *config, const pmFPAview *view, const char *filerule, psArray *objectsSrc) {
 
@@ -652,4 +821,5 @@
 	psAssert (readout, "missing readout?");
 
+	// create DETECTIONS containers for each image, in case one lacks it
 	pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
 	if (!detections) {
@@ -665,4 +835,6 @@
 	    psAssert (detections, "missing detections?");
 	}
+
+	// we need to save the new sources on the detection arrays of the appropriate image
 	detArrays->data[i] = psMemIncrRefCounter(detections);
 	readouts->data[i] = psMemIncrRefCounter(readout);
@@ -695,7 +867,8 @@
 
 	    pmSource *sourceOut = pmSourceCopy(sourceSrc);
+	    sourceOut->parent = sourceSrc;
+
+	    // save on the output object array at the same location
 	    objectOut->sources->data[i] = sourceOut;
-
-	    sourceOut->parent = sourceSrc;
 
 	    // keep the original source flags and sequence ID (if set)
@@ -720,5 +893,5 @@
 	    pmSourceFreePixels (sourceOut);
 
-	    // set the output readotu
+	    // set the output readout
 	    int index = sourceOut->imageID;
 	    if (index >= readouts->n) continue; // skip the sources generated by the chisq image
Index: trunk/psphot/src/psphotOutput.c
===================================================================
--- trunk/psphot/src/psphotOutput.c	(revision 32687)
+++ trunk/psphot/src/psphotOutput.c	(revision 32695)
@@ -407,4 +407,71 @@
     return true;
 }
+
+// for now, let's store the detections on the readout->analysis for each readout
+bool psphotDumpTest (pmConfig *config, const pmFPAview *view, const char *filerule)
+{
+    static int npass = 0;
+    char filename[64];
+
+    // XXX dump tests are disabled unless this is commented out:
+    return true;
+
+    bool status = true;
+
+    int num = psphotFileruleCount(config, filerule);
+
+    snprintf (filename, 64, "testdump.%02d.dat", npass);
+    FILE *f = fopen (filename, "w");
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+
+        // find the currently selected readout
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
+        psAssert (file, "missing file?");
+
+        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+        psAssert (readout, "missing readout?");
+
+        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+        psAssert (detections, "missing detections?");
+
+        psArray *sources = detections->newSources ? detections->newSources : detections->allSources;
+        psAssert (sources, "missing sources?");
+
+	if (detections->newSources) {
+	    fprintf (f, "## --- from new sources ---\n");
+	} else {
+	    fprintf (f, "## --- from all sources ---\n");
+	}
+
+	for (int i = 0; i < sources->n; i++) {
+	    pmSource *source = sources->data[i];
+	    if (!source) continue;
+
+	    pmPeak *peak = source->peak;
+	    if (!peak) continue;
+
+	    // XXX only dump a given region
+	    // if (peak->xf < 20) continue;
+	    // if (peak->yf < 20) continue;
+	    // if (peak->xf > 40) continue;
+	    // if (peak->yf > 70) continue;
+
+	    float Msum = source->moments ? source->moments->Sum : NAN;
+	    float Mx   = source->moments ? source->moments->Mx : NAN;
+	    float My   = source->moments ? source->moments->My : NAN;
+	    // float Npix = source->moments ? source->moments->nPixels : NAN;
+	    float Io = source->modelPSF ? source->modelPSF->params->data.F32[PM_PAR_I0] : NAN;
+	    fprintf (f, "%d %f %f  : %f %f : %f %f\n", source->imageID, peak->xf, peak->yf, Mx, My, Msum, Io);
+	}
+    }
+    fclose (f);
+    npass ++;
+
+    return true;
+}
+
+# if (0)
 bool psphotDumpTest (pmConfig *config, const pmFPAview *view, const char *filerule, char *filename) {
 
@@ -449,3 +516,3 @@
     return true;
 }
-
+# endif
Index: trunk/psphot/src/psphotReadout.c
===================================================================
--- trunk/psphot/src/psphotReadout.c	(revision 32687)
+++ trunk/psphot/src/psphotReadout.c	(revision 32695)
@@ -9,55 +9,4 @@
 }
 
-// for now, let's store the detections on the readout->analysis for each readout
-bool psphotDumpChisqs (pmConfig *config, const pmFPAview *view, const char *filerule)
-{
-    static int npass = 0;
-    char filename[64];
-
-    return true;
-
-    bool status = true;
-
-    int num = psphotFileruleCount(config, filerule);
-
-    snprintf (filename, 64, "chisq.%02d.dat", npass);
-    FILE *f = fopen (filename, "w");
-
-    // loop over the available readouts
-    for (int i = 0; i < num; i++) {
-
-        // find the currently selected readout
-        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
-        psAssert (file, "missing file?");
-
-        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
-        psAssert (readout, "missing readout?");
-
-        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
-        psAssert (detections, "missing detections?");
-
-        psArray *sources = detections->allSources;
-        psAssert (sources, "missing sources?");
-
-	for (int i = 0; i < sources->n; i++) {
-	    pmSource *source = sources->data[i];
-	    if (!source) continue;
-
-	    pmModel *model = pmSourceGetModel (NULL, source);
-	    if (!model) continue;
-	
-	    if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
-		fprintf (f, "%f %f %f %d %d %f  1 NONLINEAR\n", model->mag, model->params->data.F32[1], model->chisq, model->nDOF, model->nPix, model->chisqNorm);
-	    } else {
-		fprintf (f, "%f %f %f %d %d %f  0 LINEAR\n", model->mag, model->params->data.F32[1], model->chisq, model->nDOF, model->nPix, model->chisqNorm);
-	    }
-	}
-    }
-    fclose (f);
-    npass ++;
-
-    return true;
-}
-
 bool psphotReadout(pmConfig *config, const pmFPAview *view, const char *filerule) {
 
@@ -188,5 +137,4 @@
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
     psphotFitSourcesLinear (config, view, filerule, false); // pass 1 (detections->allSources)
-    psphotDumpChisqs (config, view, filerule);
 
     // measure the radial profiles to the sky
@@ -208,13 +156,11 @@
     // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
     psphotBlendFit (config, view, filerule); // pass 1 (detections->allSources)
-    psphotDumpChisqs (config, view, filerule);
 
     // replace all sources
-    psphotReplaceAllSources (config, view, filerule); // pass 1 (detections->allSources)
+    psphotReplaceAllSources (config, view, filerule, false); // pass 1 (detections->allSources)
 
     // linear fit to include all sources (subtract again)
     // NOTE : apply to ALL sources (extended + psf)
     psphotFitSourcesLinear (config, view, filerule, true); // pass 2 (detections->allSources)
-    psphotDumpChisqs (config, view, filerule);
 
     // if we only do one pass, skip to extended source analysis
@@ -253,5 +199,5 @@
 	// replace all sources so fit below applies to all at once
 	// NOTE: apply only to OLD sources (which have been subtracted)
-	psphotReplaceAllSources (config, view, filerule); // pass 2
+	psphotReplaceAllSources (config, view, filerule, false); // pass 2
 
 	// merge the newly selected sources into the existing list
@@ -262,5 +208,4 @@
 	// NOTE: apply to ALL sources
 	psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources)
-	psphotDumpChisqs (config, view, filerule);
     }
 
@@ -295,5 +240,5 @@
 	// replace all sources so fit below applies to all at once
 	// NOTE: apply only to OLD sources (which have been subtracted)
-	psphotReplaceAllSources (config, view, filerule); // pass 2
+	psphotReplaceAllSources (config, view, filerule, false); // pass 2
 
 	// merge the newly selected sources into the existing list
Index: trunk/psphot/src/psphotReplaceUnfit.c
===================================================================
--- trunk/psphot/src/psphotReplaceUnfit.c	(revision 32687)
+++ trunk/psphot/src/psphotReplaceUnfit.c	(revision 32695)
@@ -9,9 +9,9 @@
 
     for (int i = 0; i < sources->n; i++) {
-      source = sources->data[i];
-
-      // replace other sources?
-      if (source->mode & PM_SOURCE_MODE_FAIL) goto replace;
-      continue;
+	source = sources->data[i];
+
+	// replace other sources?
+	if (source->mode & PM_SOURCE_MODE_FAIL) goto replace;
+	continue;
 
     replace:
@@ -23,5 +23,5 @@
 
 // for now, let's store the detections on the readout->analysis for each readout
-bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule)
+bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState)
 {
     bool status = true;
@@ -35,5 +35,5 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-	if (!psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe)) {
+	if (!psphotReplaceAllSourcesReadout (config, view, filerule, i, recipe, ignoreState)) {
 	    psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i);
 	    return false;
@@ -43,8 +43,7 @@
 }
 
-bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
-
-    bool status;
-    pmSource *source;
+bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState) {
+
+    bool status;
 
     psTimerStart ("psphot.replace");
@@ -75,10 +74,24 @@
 
     for (int i = 0; i < sources->n; i++) {
-      source = sources->data[i];
-
-      // replace other sources?
-      if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
-
-      pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+	pmSource *source = sources->data[i];
+
+	if (ignoreState) {
+	    // rely on the type of source to decide if we subtract it or not
+
+	    // skip non-astronomical objects (very likely defects)
+	    if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+	    if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+      
+	    // do not include CRs in the full ensemble fit
+	    if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
+	
+	    // do not include MOMENTS_FAILURES in the fit
+	    if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+	} else {
+	    // if we respect the state, do not replace unsubtracted sources
+	    if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
+	}
+
+	pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     }
 
@@ -88,5 +101,84 @@
 }
 
-bool psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe) {
+// for now, let's store the detections on the readout->analysis for each readout
+bool psphotRemoveAllSources (pmConfig *config, const pmFPAview *view, const char *filerule, bool ignoreState)
+{
+    bool status = true;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int num = psphotFileruleCount(config, filerule);
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+	if (!psphotRemoveAllSourcesReadout (config, view, filerule, i, recipe, ignoreState)) {
+	    psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for %s entry %d", filerule, i);
+	    return false;
+	}
+    }
+    return true;
+}
+
+bool psphotRemoveAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool ignoreState) {
+
+    bool status;
+
+    psTimerStart ("psphot.replace");
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
+    psAssert (file, "missing file?");
+
+    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *sources = detections->allSources;
+    psAssert (sources, "missing sources?");
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    psAssert (maskVal, "missing mask value?");
+
+    // bit-mask to mark pixels not used in analysis
+    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
+    assert (markVal);
+
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
+
+    for (int i = 0; i < sources->n; i++) {
+	pmSource *source = sources->data[i];
+
+	if (ignoreState) {
+	    // rely on the type of source to decide if we subtract it or not
+
+	    // skip non-astronomical objects (very likely defects)
+	    if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+	    if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+      
+	    // do not include CRs in the full ensemble fit
+	    if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
+	
+	    // do not include MOMENTS_FAILURES in the fit
+	    if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+	} else {
+	    // if we respect the state, only remove unsubtracted sources
+	    if ((source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) continue;
+	}
+
+	pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+    }
+
+    psphotVisualShowImage(readout);
+    psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace"));
+    return true;
+}
+
+bool psphotRemoveAllSourcesByArray (const psArray *sources, const psMetadata *recipe) {
 
     bool status;
@@ -100,10 +192,10 @@
 
     for (int i = 0; i < sources->n; i++) {
-      source = sources->data[i];
-
-      // replace other sources?
-      if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
-
-      pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+	source = sources->data[i];
+
+	// replace other sources?
+	if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
+
+	pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     }
     psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace"));
@@ -158,17 +250,17 @@
 
     for (int i = 0; i < sources->n; i++) {
-      source = sources->data[i];
-
-      // sources have not yet been subtracted in this image (but this flag may be raised)
-      source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
-      if (!source->modelPSF) continue;
-
-      float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
-      float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
-      float radius = source->modelPSF->fitRadius;
-
-      // force a redefine to this image
-      pmSourceFreePixels(source);
-      pmSourceRedefinePixels (source, readout, Xo, Yo, radius);
+	source = sources->data[i];
+
+	// sources have not yet been subtracted in this image (but this flag may be raised)
+	source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+	if (!source->modelPSF) continue;
+
+	float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+	float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+	float radius = source->modelPSF->fitRadius;
+
+	// force a redefine to this image
+	pmSourceFreePixels(source);
+	pmSourceRedefinePixels (source, readout, Xo, Yo, radius);
     }
     return true;
@@ -229,62 +321,62 @@
 
     for (int i = 0; i < sources->n; i++) {
-      source = sources->data[i];
-
-      // *** we need to cache the 'best' model, and we have 3 cases:
-      // 1) model is the psf model --> generate from the new psf
-      // 2) model is an unconvolved extended model --> just cache the copy (not perfect)
-      // 3) model is a convolved extended model --> re-generate
-
-      // use the 'best' model to cache the model (PSF or EXT : EXT may point at one of modelFits
-      bool isPSF = false;
-      pmModel *model = pmSourceGetModel(&isPSF, source);
-      if (!model) continue;
-
-      float radius = model->fitRadius; // save for future use below
-
-      // regenerate the PSF if the model is a PSF, or if we need the PSF for a PCM
-      if (isPSF || model->isPCM) {
-	  // the guess central intensity comes from the peak:
-	  float Io = source->peak->rawFlux;
-	  float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
-	  float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
-
-	  // generate a model for this object with Io = 1.0
-	  pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io);
-	  if (modelPSF == NULL) {
-	      psWarning ("Failed to determine PSF model for source at (%f,%f), skipping", Xo, Yo);
-	      continue;
-	  }
-
-	  // set the source PSF model
-	  psFree (source->modelPSF);
-	  source->modelPSF = modelPSF;
-	  source->modelPSF->fitRadius = radius;
-      }
-
-      if (model->isPCM) {
-	  psAssert(false, "this section is not complete");
-
-	  pmSourceCachePSF (source, maskVal);
-
-	  psKernel *psfKernel = pmPCMkernelFromPSF(source, psfSize);
-	  if (!psfKernel) { 
-	      psWarning ("no psf kernel");
-	  }
-
-	  // generate an image of the right size
-	  psImage *rawModelFlux = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
-	  psImageInit (rawModelFlux, 0.0);
+	source = sources->data[i];
+
+	// *** we need to cache the 'best' model, and we have 3 cases:
+	// 1) model is the psf model --> generate from the new psf
+	// 2) model is an unconvolved extended model --> just cache the copy (not perfect)
+	// 3) model is a convolved extended model --> re-generate
+
+	// use the 'best' model to cache the model (PSF or EXT : EXT may point at one of modelFits
+	bool isPSF = false;
+	pmModel *model = pmSourceGetModel(&isPSF, source);
+	if (!model) continue;
+
+	float radius = model->fitRadius; // save for future use below
+
+	// regenerate the PSF if the model is a PSF, or if we need the PSF for a PCM
+	if (isPSF || model->isPCM) {
+	    // the guess central intensity comes from the peak:
+	    float Io = source->peak->rawFlux;
+	    float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+	    float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+
+	    // generate a model for this object with Io = 1.0
+	    pmModel *modelPSF = pmModelFromPSFforXY(psf, Xo, Yo, Io);
+	    if (modelPSF == NULL) {
+		psWarning ("Failed to determine PSF model for source at (%f,%f), skipping", Xo, Yo);
+		continue;
+	    }
+
+	    // set the source PSF model
+	    psFree (source->modelPSF);
+	    source->modelPSF = modelPSF;
+	    source->modelPSF->fitRadius = radius;
+	}
+
+	if (model->isPCM) {
+	    psAssert(false, "this section is not complete");
+
+	    pmSourceCachePSF (source, maskVal);
+
+	    psKernel *psfKernel = pmPCMkernelFromPSF(source, psfSize);
+	    if (!psfKernel) { 
+		psWarning ("no psf kernel");
+	    }
+
+	    // generate an image of the right size
+	    psImage *rawModelFlux = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
+	    psImageInit (rawModelFlux, 0.0);
 	  
-	  // insert the model image normalized to 1.0
-	  pmModelAdd (rawModelFlux, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
-
-	  psImageConvolveFFT (source->modelFlux, rawModelFlux, NULL, 0, psfKernel);
+	    // insert the model image normalized to 1.0
+	    pmModelAdd (rawModelFlux, NULL, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
+
+	    psImageConvolveFFT (source->modelFlux, rawModelFlux, NULL, 0, psfKernel);
 	  
-	  psFree (psfKernel);
-	  psFree (rawModelFlux);
-      }
-
-      pmSourceCacheModel (source, maskVal);  // ALLOC x14 (!)
+	    psFree (psfKernel);
+	    psFree (rawModelFlux);
+	}
+
+	pmSourceCacheModel (source, maskVal);  // ALLOC x14 (!)
     }
 
Index: trunk/psphot/src/psphotSourceFits.c
===================================================================
--- trunk/psphot/src/psphotSourceFits.c	(revision 32687)
+++ trunk/psphot/src/psphotSourceFits.c	(revision 32695)
@@ -257,4 +257,5 @@
 	okDBL  = psphotEvalDBL (tmpSrc, DBL->data[0]);
 	okDBL &= psphotEvalDBL (tmpSrc, DBL->data[1]);
+	okDBL = false; // XXX this is failing badly...
 	// XXX should I keep / save the flags set in the eval functions?
 
Index: trunk/psphot/src/psphotStackChisqImage.c
===================================================================
--- trunk/psphot/src/psphotStackChisqImage.c	(revision 32687)
+++ trunk/psphot/src/psphotStackChisqImage.c	(revision 32695)
@@ -13,4 +13,6 @@
     psAssert (chisqFile, "missing chisq image FPA?");
 
+    // the readout containing the chisq image is generated in the first pass of this loop and
+    // used by the successive passes
     pmReadout *chiReadout = NULL;
 
Index: trunk/psphot/src/psphotStackReadout.c
===================================================================
--- trunk/psphot/src/psphotStackReadout.c	(revision 32687)
+++ trunk/psphot/src/psphotStackReadout.c	(revision 32695)
@@ -82,8 +82,4 @@
     }
 
-    // XXX I think this is not defined correctly for an array of images.
-    // XXX I probably need to subtract the model (same model?) for both RAW and OUT.
-    // XXX But, probably not a problem in practice since the stacks are constructed with 0.0 mean level.
-
     // generate a background model (median, smoothed image)
     if (!psphotModelBackground (config, view, STACK_DET)) {
@@ -91,4 +87,7 @@
     }
     if (!psphotSubtractBackground (config, view, STACK_DET)) {
+	return psphotReadoutCleanup (config, view, STACK_SRC);
+    }
+    if (!psphotSubtractBackground (config, view, STACK_SRC)) {
 	return psphotReadoutCleanup (config, view, STACK_SRC);
     }
@@ -114,5 +113,7 @@
     }
 
-    // if DET and SRC are different images, copy the detections from DET to SRC 
+    // If DET and SRC are different images, copy the detections from DET to SRC.  This 'copy'
+    // is just a copy of the container pointer; the sources on both DET and SRC are the same
+    // memory objects
     if (strcmp(STACK_SRC, STACK_DET)) {
 	if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
@@ -130,4 +131,5 @@
 	return psphotReadoutCleanup (config, view, STACK_SRC);
     }
+    // psphotDumpTest (config, view, STACK_SRC);
     psMemDump("sourcestats");
 
@@ -166,6 +168,4 @@
 
     // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
-    // XXX why do this as a stack operation?
-    // psphotFitSourcesLinearStack (config, objects, false);
     psphotFitSourcesLinear (config, view, STACK_SRC, false);
     psphotStackVisualFilerule(config, view, STACK_SRC);
@@ -186,6 +186,6 @@
     psphotBlendFit (config, view, STACK_SRC); // pass 1 (detections->allSources)
 
-    // replace all sources
-    psphotReplaceAllSources (config, view, STACK_SRC); // pass 1 (detections->allSources)
+    // replace all sources (do NOT ignore subtraction state)
+    psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 1 (detections->allSources)
 
     // if we only do one pass, skip to extended source analysis
@@ -194,4 +194,5 @@
     // linear fit to include all sources (subtract again)
     // NOTE : apply to ALL sources (extended + psf)
+    // NOTE 2 : this function subtracts the models from the given filerule (SRC), not DET
     psphotFitSourcesLinear (config, view, STACK_SRC, true); // pass 2 (detections->allSources)
 
@@ -200,4 +201,16 @@
     // NOTE: this block performs the 2nd pass low-significance PSF detection stage
     { 
+	// if DET and SRC are different images, generate children sources for all sources in
+	// the SRC image.  This operation replaces the existing DETECTION container on DET
+	// which is currently a view to the one on SRC).  children sources go to
+	// det->allSources
+	if (strcmp(STACK_SRC, STACK_DET)) {
+	    psphotSourceChildren (config, view, STACK_DET, STACK_SRC); 
+
+	    //  subtract all sources from DET (this will subtract using the psf model for SRC, which
+	    //  will somewhat oversubtract the sources -- this is OK
+	    psphotRemoveAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources
+	}
+
 	// add noise for subtracted objects
 	psphotAddNoise (config, view, STACK_DET); // pass 1 (detections->allSources)
@@ -212,7 +225,11 @@
 
 	// if DET and SRC are different images, copy the detections from DET to SRC 
+	// (this operation just ensures the metadata container has a view on SRC as well
 	if (strcmp(STACK_SRC, STACK_DET)) {
-	    // XXX how does this handle 1st vs 2nd pass sources?
-	    if (!psphotCopySources (config, view, STACK_SRC, STACK_DET)) {
+	    // replace all sources in DET
+	    psphotReplaceAllSources (config, view, STACK_DET, false); // ignore subtraction state for sources
+
+	    // copy the newly detected peaks from DET to SRC so SourceStats below can operate on them
+	    if (!psphotCopyPeaks (config, view, STACK_SRC, STACK_DET)) {
 		psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
 		return psphotReadoutCleanup (config, view, STACK_SRC);
@@ -237,5 +254,5 @@
 	// replace all sources so fit below applies to all at once
 	// NOTE: apply only to OLD sources (which have been subtracted)
-	psphotReplaceAllSources (config, view, STACK_SRC); // pass 2
+	psphotReplaceAllSources (config, view, STACK_SRC, false); // pass 2
 
 	// merge the newly selected sources into the existing list
@@ -312,5 +329,5 @@
 
 	// replace the flux in the image so it is returned to its original state
-	psphotReplaceAllSources (config, view, STACK_OUT);
+	psphotReplaceAllSources (config, view, STACK_OUT, false);
 
 	// smooth to the next FWHM, or set 'smoothAgain' to false if no more