Index: trunk/psphot/src/psphotBlendFit.c
===================================================================
--- trunk/psphot/src/psphotBlendFit.c	(revision 34258)
+++ trunk/psphot/src/psphotBlendFit.c	(revision 34404)
@@ -1,3 +1,4 @@
 # include "psphotInternal.h"
+bool psphotBlendFitSetSource(pmSource *source, psImage *IDimage, float sigSigma);
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -92,4 +93,17 @@
     float skySig = psMetadataLookupF32(&status, recipe, "SKY_SIG");
     assert (status && isfinite(skySig) && skySig > 0);
+
+# if (0)
+    // **** test block : generate an ID image where pixels are set based the source models (flux > 0.1 peak && flux > 2.0 skySig)
+    psImage *IDimage = psImageAlloc (readout->image->numCols, readout->image->numRows, PS_TYPE_S32);
+    psImageInit (IDimage, 0);
+
+    // start with the currently known moments (Mxx, Myy, Mxy) and generate a window image
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+	psphotBlendFitSetSource (source, IDimage, skySig);
+    }
+    psphotSaveImage (NULL, IDimage, "idimage.fits");
+# endif
 
     // Define source fitting parameters for extended source fits
@@ -241,12 +255,13 @@
         pmSource *source = sources->data[i];
 
-# if (0)
-# define TEST_X 34
-# define TEST_Y 28
-   
+	int TEST_ON = false;
+# if (1)
+# define TEST_X 653
+# define TEST_Y 466
 	if ((fabs(source->peak->xf - TEST_X) < 5) && (fabs(source->peak->yf - TEST_Y) < 5)) {
 	    fprintf (stderr, "test object\n");
-	}
-
+	    //  psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
+	    TEST_ON = true;
+	}
 # undef TEST_X
 # undef TEST_Y
@@ -258,4 +273,7 @@
         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+
+	// skip saturated stars modeled with a radial profile 
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
 
         // skip DBL second sources (ie, added by psphotFitBlob
@@ -306,4 +324,7 @@
         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
             if (psphotFitBlob (readout, source, newSources, psf, fitOptions, maskVal, markVal)) {
+		if (TEST_ON) {
+		    psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
+		}
                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
                 Next ++;
@@ -312,4 +333,7 @@
         } else {
             if (psphotFitBlend (readout, source, psf, fitOptions, maskVal, markVal)) {
+		if (TEST_ON) {
+		    psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
+		}
                 source->type = PM_SOURCE_TYPE_STAR;
                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
@@ -349,2 +373,82 @@
 }
 
+bool psphotBlendFitSetSource(pmSource *source, psImage *IDimage, float skySigma) {
+
+    if (!source) return false;
+    if (!source->peak) return false; // XXX how can we have a peak-less source?
+    if (!source->moments) return false;
+    if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
+    if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
+    if (source->mode2 == PM_SOURCE_MODE2_MATCHED) return false;
+    psAssert(IDimage, "need a window");
+
+    if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) {
+	// find the radius at which we hit something like 0.1*SATURATION
+	return false;
+    }
+
+    if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false;
+
+    // we have a source with moments Mx, My, Mxx, Mxy, Myy.  we just need to define a Gaussian that has 
+    // these values (and peak of 1.0)
+
+    int Nx = IDimage->numCols;
+    int Ny = IDimage->numRows;
+
+    float Xo = source->moments->Mx;
+    float Yo = source->moments->My;
+
+    psEllipseMoments moments;
+    moments.x2 = source->moments->Mxx;
+    moments.y2 = source->moments->Myy;
+    moments.xy = source->moments->Mxy;
+
+    psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0);
+    if (! (isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta)) ) {
+        fprintf (stderr, "invalid axes found id: %4d major: %f minor: %f theta: %f\n", source->id, axes.major, axes.minor, axes.theta);
+        return false;
+    }
+
+    // Use 1st radial moment as size, not sigma.  Why this factor of 0.5 ?
+    float scale = 0.5 * source->moments->Mrf / axes.major;
+    axes.major *= scale;
+    axes.minor *= scale;
+
+    psEllipseShape shape = psEllipseAxesToShape(axes);
+    if (! (isfinite(shape.sx) && isfinite(shape.sy) && isfinite(shape.sxy)) ) {
+        fprintf (stderr, "invalid shape found id: %d sx: %f sy %f sxy: %f\n", source->id, shape.sx, shape.sy, shape.sxy);
+        return false;
+    }
+
+    float Smajor = axes.major;
+
+    // we want to go to 2.15 sigma (0.1 Io)
+    int minX = PS_MIN(PS_MAX(Xo - 2.15*Smajor, 0), Nx - 1);
+    int maxX = PS_MIN(PS_MAX(Xo + 2.15*Smajor, 0), Nx - 1);
+    int minY = PS_MIN(PS_MAX(Yo - 2.15*Smajor, 0), Ny - 1);
+    int maxY = PS_MIN(PS_MAX(Yo + 2.15*Smajor, 0), Ny - 1);
+
+    float rMxx = 0.5 / PS_SQR(shape.sx);
+    float rMyy = 0.5 / PS_SQR(shape.sy);
+    float Sxy = -1. * shape.sxy;    // factor of -1 is included to match the previous window function
+                                    // implementation. XXX: Is this correct?
+
+    int ID = source->id;
+    float Io = source->peak->rawFlux;
+    for (int iy = minY; iy < maxY; iy++) {
+	for (int ix = minX; ix < maxX; ix++) {
+
+	    float dX = (ix + 0.5 - Xo);
+	    float dY = (iy + 0.5 - Yo);
+
+	    float z = rMxx * PS_SQR(dX) + rMyy * PS_SQR(dY) + Sxy*dX*dY;
+
+	    float f = Io*exp(-z);
+	    
+	    if (f < 2.0*skySigma) continue;
+            IDimage->data.S32[iy][ix] = ID;
+	}
+    }
+
+    return true;
+}