Index: /trunk/psphot/src/Makefile.am
===================================================================
--- /trunk/psphot/src/Makefile.am	(revision 12664)
+++ /trunk/psphot/src/Makefile.am	(revision 12665)
@@ -49,9 +49,10 @@
 	psphotMergeSources.c	 \
 	psphotReadoutCleanup.c	 \
+	psphotSourcePlots.c	 \
+	psphotRadialPlot.c	 \
+	psphotMosaicSubimage.c	 \
 	psphotAddNoise.c	 \
 	psphotCleanup.c	   	
 
-#	psphotSourcePlots.c	 
-#	psphotRadialPlot.c	 
 
 include_HEADERS = \
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 12664)
+++ /trunk/psphot/src/psphot.h	(revision 12665)
@@ -104,4 +104,9 @@
 bool psphotSetHeaderNstars (psMetadata *recipe, psArray *sources);
 bool psphotAddNoise (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add);
-bool psphotRadialPlot (const char *filename, pmSource *source);
+bool psphotRadialPlot (int *kapa, char *filename, pmSource *source);
 bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe);
+bool psphotMosaicSubimage (psImage *outImage, pmSource *source, int Xo, int Yo, int DX, int DY);
+
+bool psphotAddWithTest (pmSource *source, bool useState);
+bool psphotSubWithTest (pmSource *source, bool useState);
+bool psphotSetState (pmSource *source, bool curState);
Index: /trunk/psphot/src/psphotChoosePSF.c
===================================================================
--- /trunk/psphot/src/psphotChoosePSF.c	(revision 12664)
+++ /trunk/psphot/src/psphotChoosePSF.c	(revision 12665)
@@ -74,7 +74,14 @@
 
     // select the candidate PSF stars (pointers to original sources)
-    for (int i = 0; (i < sources->n) && (stars->n < NSTARS); i++) {
+    for (int i = 0; i < sources->n; i++) {
         pmSource *source = sources->data[i];
-        if (source->mode & PM_SOURCE_MODE_PSFSTAR) psArrayAdd (stars, 200, source);
+        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
+	    // keep NSTARS PSF stars, unmark the rest
+	    if (stars->n < NSTARS) {
+		psArrayAdd (stars, 200, source);
+	    } else {
+		source->mode &= ~PM_SOURCE_MODE_PSFSTAR;
+	    }
+	} 
     }
     psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "selected candidate %ld PSF objects\n", stars->n);
Index: /trunk/psphot/src/psphotMosaicSubimage.c
===================================================================
--- /trunk/psphot/src/psphotMosaicSubimage.c	(revision 12665)
+++ /trunk/psphot/src/psphotMosaicSubimage.c	(revision 12665)
@@ -0,0 +1,45 @@
+# include "psphot.h"
+
+// insert the source image into the outimage at Xo, Yo
+bool psphotMosaicSubimage (psImage *outImage, pmSource *source, int Xo, int Yo, int DX, int DY) {
+
+    psRegion inRegion, outRegion;
+    psImage *inImage = source->pixels;
+
+    // identify the region in the output image
+    outRegion = psRegionSet (Xo, Xo + DX, Yo, Yo + DX);
+    outRegion = psRegionForImage (outImage, outRegion);
+    int DXo = outRegion.x1 - outRegion.x0;
+    int DYo = outRegion.y1 - outRegion.y0;
+    if (DXo <= 0) return false;
+    if (DYo <= 0) return false;
+    
+    // center the input source in the output box
+    int dX = (DXo - 1) / 2;
+    int dY = (DYo - 1) / 2;
+
+    // int xo = inImage->col0 + inImage->numCols / 2;
+    // int yo = inImage->row0 + inImage->numRows / 2;
+
+    int xo = source->peak->xf;
+    int yo = source->peak->yf;
+
+    // adjust region to overlay input image pixels
+    inRegion = psRegionSet (xo - dX, xo + dX + 1, yo - dY, yo + dY + 1);
+    inRegion = psRegionForImage (inImage, inRegion);
+
+    float peak = source->peak->flux;
+
+    psImage *subImage = psImageSubset (inImage, inRegion);
+    psImage *newImage = psImageAlloc (subImage->numCols, subImage->numRows, PS_TYPE_F32);
+    for (int iy = 0; iy < newImage->numRows; iy++) {
+	for (int ix = 0; ix < newImage->numCols; ix++) {
+	    newImage->data.F32[iy][ix] = subImage->data.F32[iy][ix] / peak;
+	}
+    }
+
+    psImageOverlaySection (outImage, newImage, Xo, Yo, "=");
+
+    psFree (newImage);
+    return true;
+}
Index: /trunk/psphot/src/psphotRadialPlot.c
===================================================================
--- /trunk/psphot/src/psphotRadialPlot.c	(revision 12664)
+++ /trunk/psphot/src/psphotRadialPlot.c	(revision 12665)
@@ -6,47 +6,66 @@
 # include <kapa.h>
 
-bool psphotRadialPlot (const char *filename, pmSource *source) {
+static int nCount = 0;
+
+bool psphotRadialPlot (int *kapa, char *filename, pmSource *source) {
+
+    Graphdata graphdata;
+
+    // only plot 50 stars for now...
+    if (nCount > 00) {
+	if (*kapa != 0) {
+	    KiiClose (*kapa);
+	    *kapa = 0;
+	}
+	return true;
+    }
 
     // XXX get the 'showWindow' option from the recipes somewhere
-    int kapa = pmKapaOpen (false);
-    if (kapa == -1) {
+    // XXX 'showWindow = false' is broken
+    if (*kapa == 0) {
+	*kapa = pmKapaOpen (true);
+	KapaResize (*kapa, 500, 500);
+	unlink (filename);
+    }
+    if (*kapa == -1) {
 	psError(PSPHOT_ERR_UNKNOWN, true, "failure to open kapa");
 	return false;
     }
 
-    KapaResize (kapa, 500, 500);
     KapaInitGraph (&graphdata);
+    KapaClear(*kapa, true);
 
     // examine sources to set data range
-    graphdata.xmin = -0.05;
+    graphdata.xmin =  -0.05;
+    graphdata.xmax = +30.05;
     graphdata.ymin = -0.05;
-    graphdata.xmax = +2.05;
-    graphdata.ymax = +2.05;
-    KapaSetLimits (kapa, &graphdata);
+    graphdata.ymax = +5.05;
+    KapaSetLimits (*kapa, &graphdata);
   
-    KapaSetFont (kapa, "helvetica", 14);
-    KapaBox (kapa, &graphdata);
-    KapaSendLabel (kapa, "&ss&h_x| (pixels)", KAPA_LABEL_XM);
-    KapaSendLabel (kapa, "&ss&h_y| (pixels)", KAPA_LABEL_YM);
+    KapaSetFont (*kapa, "helvetica", 14);
+    KapaBox (*kapa, &graphdata);
+    KapaSendLabel (*kapa, "radius (pixels)", KAPA_LABEL_XM);
+    KapaSendLabel (*kapa, "log flux (counts)", KAPA_LABEL_YM);
 	       
     int nPts = source->pixels->numRows * source->pixels->numCols;
-    psVector *radius = psVectorAllocEmpty (nPts, PS_TYPE_F32);
-    psVector *lflux  = psVectorAllocEmpty (nPts, PS_TYPE_F32);
-    psVector *rmask  = psVectorAllocEmpty (nPts, PS_TYPE_F32);
-    psVector *lfmask = psVectorAllocEmpty (nPts, PS_TYPE_F32);
+    psVector *rg = psVectorAllocEmpty (nPts, PS_TYPE_F32);
+    psVector *fg = psVectorAllocEmpty (nPts, PS_TYPE_F32);
+    psVector *rb = psVectorAllocEmpty (nPts, PS_TYPE_F32);
+    psVector *fb = psVectorAllocEmpty (nPts, PS_TYPE_F32);
 
-    int n = 0;
-    int nm = 0;
+    int ng = 0;
+    int nb = 0;
+    float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
+    float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
     for (int iy = 0; iy < source->pixels->numRows; iy++) {
 	for (int ix = 0; ix < source->pixels->numCols; ix++) {
-	    
 	    if (source->mask->data.U8[iy][ix]) {
-		rmask->data.F32[nm] = RADIUS;
-		lfmask->data.F32[nm] = log10(source->pixels->data.F32[iy][ix]);
-		nm++;
+		rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ;
+		fb->data.F32[nb] = log10(source->pixels->data.F32[iy][ix]);
+		nb++;
 	    } else {
-		radius->data.F32[n] = RADIUS;
-		lflux->data.F32[n] = log10(source->pixels->data.F32[iy][ix]);
-		n++;
+		rg->data.F32[ng] = hypot (ix - Xo, iy - Yo) ;
+		fg->data.F32[ng] = log10(source->pixels->data.F32[iy][ix]);
+		ng++;
 	    }
 	}
@@ -56,26 +75,44 @@
 
     graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 2;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    KapaPrepPlot (*kapa, ng, &graphdata);
+    KapaPlotVector (*kapa, ng, rg->data.F32);
+    KapaPlotVector (*kapa, ng, fg->data.F32);
+
+    graphdata.color = KapaColorByName ("red");
     graphdata.ptype = 0;
     graphdata.size = 0.3;
     graphdata.style = 2;
-    KapaPrepPlot (kapa, nm, &graphdata);
-    KapaPlotVector (kapa, nm, rmask->data.F32);
-    KapaPlotVector (kapa, nm, lfmask->data.F32);
+    KapaPrepPlot (*kapa, nb, &graphdata);
+    KapaPlotVector (*kapa, nb, rb->data.F32);
+    KapaPlotVector (*kapa, nb, fb->data.F32);
   
-    graphdata.color = KapaColorByName ("red");
-    graphdata.ptype = 2;
-    graphdata.size = 0.5;
-    graphdata.style = 2;
-    KapaPrepPlot (kapa, n, &graphdata);
-    KapaPlotVector (kapa, n, radius->data.F32);
-    KapaPlotVector (kapa, n, lflux->data.F32);
+    psLogMsg ("psphot", 3, "saving plot to %s", filename);
 
-    psLogMsg ("psphot", 3, "saving plot to %s", file->filename);
-    KapaPS (kapa, false, KAPA_PS_NEWPAGE, filename, "none");
+    char pagename[16];
+    sprintf (pagename, "%02d", nCount);
+    if (nCount == 0) {
+	KiiPS (*kapa, false, KAPA_PS_NEWPLOT, filename, pagename);
+    } else {
+	KiiPS (*kapa, false, KAPA_PS_NEWPAGE, filename, pagename);
+    }
 
-    psFree (radius);
-    psFree (lflux);
-    psFree (rmask);
-    psFree (lfmask);
+    psFree (rg);
+    psFree (fg);
+    psFree (rb);
+    psFree (fb);
 
+    nCount ++;
+    return true;
 }
+
+# else
+
+bool psphotRadialPlot (int *kapa, const char *filename, pmSource *source) {
+    psLogMsg ("psphot", 3, "skipping source radial plots");
+    return true;
+}
+
+# endif
Index: /trunk/psphot/src/psphotReadout.c
===================================================================
--- /trunk/psphot/src/psphotReadout.c	(revision 12664)
+++ /trunk/psphot/src/psphotReadout.c	(revision 12665)
@@ -108,4 +108,7 @@
     }
 
+    // plot positive + negative sources (replace, then remove)
+    psphotSourcePlots (readout, sources, recipe);
+
     // non-linear PSF and EXT fit to brighter sources
     psphotBlendFit (readout, sources, recipe, psf);
@@ -113,4 +116,7 @@
     // replace all sources
     psphotReplaceAll (sources);
+
+    // plot positive sources 
+    // psphotSourcePlots (readout, sources, recipe);
 
     // linear PSF fit to remaining peaks
Index: /trunk/psphot/src/psphotReplaceUnfit.c
===================================================================
--- /trunk/psphot/src/psphotReplaceUnfit.c	(revision 12664)
+++ /trunk/psphot/src/psphotReplaceUnfit.c	(revision 12665)
@@ -53,2 +53,67 @@
     return true;
 }
+
+// add or sub source replace or if the source has 
+bool psphotAddWithTest (pmSource *source, bool useState) { 
+
+    // what is current state? (true : add; false : sub)
+    bool state = !(source->mode & PM_SOURCE_MODE_SUBTRACTED);
+    if (state && useState) return true;
+
+    // select appropriate model
+    pmModel *model = pmSourceGetModel (NULL, source);
+    if (model == NULL) return false;  // model must be defined
+    
+    psTrace ("psphot", 3, "replacing object at %f,%f\n", 
+	     model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+    
+    // replace the model if 1) state says it is missing or 2) useState is false (just do it)
+    if (!state || !useState) {
+	pmModelAdd (source->pixels, source->mask, model, false, false);
+    }
+    return true;
+}
+
+// add or sub source replace or if the source has 
+bool psphotSubWithTest (pmSource *source, bool useState) { 
+
+    // what is current state? (true : sub; false : add)
+    bool state = (source->mode & PM_SOURCE_MODE_SUBTRACTED);
+    if (state && useState) return true;
+
+    // select appropriate model
+    pmModel *model = pmSourceGetModel (NULL, source);
+    if (model == NULL) return false;  // model must be defined
+    
+    psTrace ("psphot", 3, "replacing object at %f,%f\n", 
+	     model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+    
+    // replace the model if 1) state says it is missing or 2) useState is false (just do it)
+    if (!state || !useState) {
+	pmModelSub (source->pixels, source->mask, model, false, false);
+    }
+    return true;
+}
+
+// add or sub source replace or if the source has 
+bool psphotSetState (pmSource *source, bool curState) { 
+
+    // what is desired state? (true : add; false : sub)
+    bool newState = !(source->mode & PM_SOURCE_MODE_SUBTRACTED);
+    if (curState == newState) return true;
+
+    // select appropriate model
+    pmModel *model = pmSourceGetModel (NULL, source);
+    if (model == NULL) return false;  // model must be defined
+    
+    psTrace ("psphot", 3, "replacing object at %f,%f\n", 
+	     model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+    
+    if (curState && !newState) {
+	pmModelSub (source->pixels, source->mask, model, false, false);
+    }
+    if (newState && !curState) {
+	pmModelAdd (source->pixels, source->mask, model, false, false);
+    }
+    return true;
+}
Index: /trunk/psphot/src/psphotSourcePlots.c
===================================================================
--- /trunk/psphot/src/psphotSourcePlots.c	(revision 12664)
+++ /trunk/psphot/src/psphotSourcePlots.c	(revision 12665)
@@ -3,10 +3,13 @@
 bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe) {
 
-    bool status = false;
+    // bool status = false;
     psTimerStart ("psphot");
 
     // the source images are written to an image 10x the size of a PSF object
-    float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
-    PS_ASSERT (status);
+    // float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
+    // PS_ASSERT (status, false);
+
+    int DX = 21;
+    int DY = 21;
 
     // examine PSF sources in S/N order (brightest first)
@@ -16,9 +19,8 @@
     int dX = 0;				// starting corner of next box
     int dY = 0;				// height of row so far
-    int NX = 10*(2*OUTER+1);		// full width of output image
+    int NX = 20*DX;			// full width of output image
     int NY = 0;				// total height of output image
 
     // first, examine the PSF and SAT stars:
-    // - generate radial plots (PS plots)
     // - determine bounding boxes for summary image
     for (int i = 0; i < sources->n; i++) {
@@ -32,9 +34,12 @@
 
 	// how does this subimage get placed into the output image?
-	if (dX + source->pixels->numCols >= NX) {
+	// DX = source->pixels->numCols
+	// DY = source->pixels->numRows
+
+	if (dX + DX > NX) {
 	    // too wide for the rest of this row
 	    if (dX == 0) {
 		// alone on this row
-		NY += source->pixels->numRows;
+		NY += DY;
 		dX = 0;
 		dY = 0;
@@ -42,13 +47,12 @@
 		// start the next row
 		NY += dY;
-		dX = source->pixels->numCols;
-		dY = source->pixels->numRows;
+		dX = DX;
+		dY = DY;
 	    }
 	} else {
 	    // extend this row
-	    dX += source->pixels->numCols;
-	    dY = PS_MAX (dY, source->pixels->numRows);
+	    dX += DX;
+	    dY = PS_MAX (dY, DY);
 	}
-	psphotRadialPlot ("radial.plots.ps", source);
     }
 
@@ -59,9 +63,13 @@
     int Xo = 0;				// starting corner of next box
     int Yo = 0;				// starting corner of next box
-    int dY = 0;				// height of row so far
+    dY = 0;				// height of row so far
+
+    int nPSF = 0;
+    int nSAT = 0;
+    int kapa = 0;			// file descriptor for plotting routine
 
     // first, examine the PSF and SAT stars:
     // - generate radial plots (PS plots)
-    // - determine bounding boxes for summary image
+    // - create output image array
     for (int i = 0; i < sources->n; i++) {
 
@@ -69,15 +77,31 @@
 
 	bool keep = false;
-        keep |= (source->mode & PM_SOURCE_MODE_PSFSTAR);
-        keep |= (source->mode & PM_SOURCE_MODE_SATSTAR);
+        if (source->mode & PM_SOURCE_MODE_PSFSTAR) {
+	    nPSF ++;
+	    keep = true;
+	}
+        if (source->mode & PM_SOURCE_MODE_SATSTAR) {
+	    nSAT ++;
+	    keep = true;
+	}	    
 	if (!keep) continue;
 
 	// how does this subimage get placed into the output image?
-	if (Xo + source->pixels->numCols >= NX) {
+	// DX = source->pixels->numCols
+	// DY = source->pixels->numRows
+
+	if (Xo + DX > NX) {
 	    // too wide for the rest of this row
 	    if (Xo == 0) {
-		// alone on this row
-		PlaceImage (source, Xo, Yo);
-		Yo += source->pixels->numRows;
+		// place source alone on this row
+		psphotAddWithTest (source, true); // replace source if subtracted
+		psphotRadialPlot (&kapa, "radial.plots.ps", source);
+		psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY);
+
+		psphotSubWithTest (source, false); // remove source (force)
+		psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY);
+
+		psphotSetState (source, false); // replace source (has been subtracted)
+		Yo += DY;
 		Xo = 0;
 		dY = 0;
@@ -86,17 +110,36 @@
 		Yo += dY;
 		Xo = 0;
-		PlaceImage (source, Xo, Yo);
-		Xo = source->pixels->numCols;
-		dY = source->pixels->numRows;
+		psphotAddWithTest (source, true); // replace source if subtracted
+		psphotRadialPlot (&kapa, "radial.plots.ps", source);
+		psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY);
+
+		psphotSubWithTest (source, false); // remove source (force)
+		psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY);
+		psphotSetState (source, false); // replace source (has been subtracted)
+
+		Xo = DX;
+		dY = DY;
 	    }
 	} else {
 	    // extend this row
-	    PlaceImage (source, Xo, Yo);
-	    Xo += source->pixels->numCols;
-	    dY = PS_MAX (dY, source->pixels->numRows);
+	    psphotAddWithTest (source, true); // replace source if subtracted
+	    psphotRadialPlot (&kapa, "radial.plots.ps", source);
+	    psphotMosaicSubimage (outpos, source, Xo, Yo, DX, DY);
+
+	    psphotSubWithTest (source, false); // remove source (force)
+	    psphotMosaicSubimage (outsub, source, Xo, Yo, DX, DY);
+	    psphotSetState (source, false); // replace source (has been subtracted)
+
+	    Xo += DX;
+	    dY = PS_MAX (dY, DY);
 	}
     }
 
-    psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d failed: %f sec\n", sources->n, Nmoments, Nfail, psTimerMark ("psphot"));
+    psphotSaveImage (NULL, outpos, "outpos.fits");
+    psphotSaveImage (NULL, outsub, "outsub.fits");
+    psLogMsg ("psphot", PS_LOG_INFO, "plotted %d sources (%d psf, %d sat): %f sec\n", nPSF + nSAT, nPSF, nSAT, psTimerMark ("psphot"));
+
+    psFree (outpos);
+    psFree (outsub);
     return (sources);
 }