Index: trunk/psModules/src/objects/Makefile.am
===================================================================
--- trunk/psModules/src/objects/Makefile.am	(revision 27530)
+++ trunk/psModules/src/objects/Makefile.am	(revision 27531)
@@ -22,4 +22,5 @@
 	pmSourceMasks.c \
 	pmSourceMoments.c \
+	pmSourceDiffStats.c \
 	pmSourceExtendedPars.c \
 	pmSourceUtils.c \
@@ -40,4 +41,5 @@
 	pmSourceIO_CMF_PS1_V1.c \
 	pmSourceIO_CMF_PS1_V2.c \
+	pmSourceIO_CMF_PS1_DV1.c \
 	pmSourceIO_MatchedRefs.c \
 	pmSourcePlots.c \
@@ -79,4 +81,5 @@
 	pmSource.h \
 	pmSourceMasks.h \
+	pmSourceDiffStats.h \
 	pmSourceExtendedPars.h \
 	pmSourceUtils.h \
Index: trunk/psModules/src/objects/models/pmModel_PGAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_PGAUSS.c	(revision 27530)
+++ trunk/psModules/src/objects/models/pmModel_PGAUSS.c	(revision 27531)
@@ -236,5 +236,5 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psF64 z, f;
+    psF64 z;
     int Nstep = 0;
     psEllipseShape shape;
@@ -271,19 +271,22 @@
     z0 = 0.0;
 
-    // perform a type of bisection to find the value
-    float f0 = 1.0 / (1 + z0 + z0*z0/2.0 + z0*z0*z0/6.0);
-    float f1 = 1.0 / (1 + z1 + z1*z1/2.0 + z1*z1*z1/6.0);
-    while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
-        z = 0.5*(z0 + z1);
-        f = 1.0 / (1 + z + z*z/2.0 + z*z*z/6.0);
-        if (f > limit) {
-            z0 = z;
-            f0 = f;
-        } else {
-            z1 = z;
-            f1 = f;
-        }
-        Nstep ++;
-    }
+    // starting guess:
+    z = 0.5*(z0 + z1);
+    float dz = 1.0;
+
+    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
+	// use Newton-Raphson to minimize f(z) - limit = 0
+	float dqdz = (1.0 + z + z*z/2.0);
+	float q = (dqdz + z*z*z/6.0);
+
+	float f = 1.0 / q;
+	float dfdz = -dqdz * f / q;
+
+	dz = (f - limit) / dfdz;
+
+	// fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
+	z -= dz;
+    }
+
     psF64 radius = sigma * sqrt (2.0 * z);
 
Index: trunk/psModules/src/objects/models/pmModel_PS1_V1.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_PS1_V1.c	(revision 27530)
+++ trunk/psModules/src/objects/models/pmModel_PS1_V1.c	(revision 27531)
@@ -267,5 +267,5 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psF64 z, f;
+    psF64 z;
     int Nstep = 0;
     psEllipseShape shape;
@@ -273,16 +273,8 @@
     psF32 *PAR = params->data.F32;
 
-    if (flux <= 0) {
-        return 1.0;
-    }
-    if (PAR[PM_PAR_I0] <= 0) {
-        return 1.0;
-    }
-    if (flux >= PAR[PM_PAR_I0]) {
-        return 1.0;
-    }
-    if (PAR[PM_PAR_7] == 0.0) {
-        return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
-    }
+    if (flux <= 0) return 1.0;
+    if (PAR[PM_PAR_I0] <= 0) return 1.0;
+    if (flux >= PAR[PM_PAR_I0]) return 1.0;
+    if (PAR[PM_PAR_7] == 0.0) return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
 
     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
@@ -305,18 +297,20 @@
     if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0;
 
-    // perform a type of bisection to find the value
-    float f0 = 1.0 / (1 + PAR[PM_PAR_7]*z0 + pow(z0, ALPHA));
-    float f1 = 1.0 / (1 + PAR[PM_PAR_7]*z1 + pow(z1, ALPHA));
-    while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
-        z = 0.5*(z0 + z1);
-        f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
-        if (f > limit) {
-            z0 = z;
-            f0 = f;
-        } else {
-            z1 = z;
-            f1 = f;
-        }
-        Nstep ++;
+    // starting guess:
+    z = 0.5*(z0 + z1);
+    float dz = 1.0;
+
+    // use Newton-Raphson to minimize f(z) - limit = 0
+    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
+	float q = (1.0 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
+	float dqdz = (PAR[PM_PAR_7] + ALPHA*pow(z, ALPHA - 1.0));
+
+	float f = 1.0 / q;
+	float dfdz = -dqdz * f / q;
+
+	dz = (f - limit) / dfdz;
+
+	// fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
+	z -= dz;
     }
     psF64 radius = sigma * sqrt (2.0 * z);
@@ -350,5 +344,4 @@
     // convert to shape terms (SXX,SYY,SXY)
     if (!pmPSF_FitToModel (out, 0.1)) {
-        // psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         psTrace("psModules.objects", 5, "Failed to fit object at (r,c) = (%.1f,%.1f)", in[PM_PAR_YPOS], in[PM_PAR_XPOS]);
         return false;
@@ -475,5 +468,4 @@
 }
 
-
 # undef PM_MODEL_FUNC
 # undef PM_MODEL_FLUX
Index: trunk/psModules/src/objects/models/pmModel_QGAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 27530)
+++ trunk/psModules/src/objects/models/pmModel_QGAUSS.c	(revision 27531)
@@ -39,4 +39,7 @@
 # define PM_MODEL_SET_LIMITS      pmModelSetLimits_QGAUSS
 
+# define ALPHA   2.250
+# define ALPHA_M 1.250
+
 // the model is a function of the pixel coordinate (pixcoord[0,1] = x,y)
 // 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords
@@ -44,19 +47,21 @@
 
 // Lax parameter limits
-static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.1 };
+static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -1.0 };
 static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
 
 // Moderate parameter limits
-static float *paramsMinModerate = paramsMinLax;
-static float *paramsMaxModerate = paramsMaxLax;
+// Tolerate a small divot (k < 0)
+static float paramsMinModerate[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, -0.05 };
+static float paramsMaxModerate[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
 
 // Strict parameter limits
-static float *paramsMinStrict = paramsMinLax;
-static float *paramsMaxStrict = paramsMaxLax;
+// k = PAR_7 < 0 is very undesirable (big divot in the middle)
+static float paramsMinStrict[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 0.0 };
+static float paramsMaxStrict[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 20.0 };
 
 // Parameter limits to use
 static float *paramsMinUse = paramsMinLax;
 static float *paramsMaxUse = paramsMaxLax;
-static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5 };
+static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5, 2.0 };
 
 static bool limitsApply = true;         // Apply limits?
@@ -81,5 +86,5 @@
     assert (z >= 0);
 
-    psF32 zp = pow(z,1.25);
+    psF32 zp = pow(z,ALPHA_M);
     psF32 r  = 1.0 / (1 + PAR[PM_PAR_7]*z + z*zp);
 
@@ -92,5 +97,5 @@
         // note difference from a pure gaussian: q = params->data.F32[PM_PAR_I0]*r
         psF32 t = r1*r;
-        psF32 q = t*(PAR[PM_PAR_7] + 2.25*zp);
+        psF32 q = t*(PAR[PM_PAR_7] + ALPHA*zp);
 
         dPAR[PM_PAR_SKY]  = +1.0;
@@ -246,7 +251,7 @@
     float f1, f2;
     for (z = DZ; z < 50; z += DZ) {
-        f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
+        f1 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
         z += DZ;
-        f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
+        f2 = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
         norm += f0 + 4*f1 + f2;
         f0 = f2;
@@ -263,5 +268,5 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psF64 z, f;
+    psF64 z;
     int Nstep = 0;
     psEllipseShape shape;
@@ -269,10 +274,8 @@
     psF32 *PAR = params->data.F32;
 
-    if (flux <= 0)
-        return (1.0);
-    if (PAR[PM_PAR_I0] <= 0)
-        return (1.0);
-    if (flux >= PAR[PM_PAR_I0])
-        return (1.0);
+    if (flux <= 0) return 1.0;
+    if (PAR[PM_PAR_I0] <= 0) return 1.0;
+    if (flux >= PAR[PM_PAR_I0]) return 1.0;
+    if (PAR[PM_PAR_7] == 0.0) return powf(PAR[PM_PAR_I0] / flux - 1.0, 1.0 / ALPHA);
 
     shape.sx  = PAR[PM_PAR_SXX] / M_SQRT2;
@@ -290,25 +293,25 @@
 
     // choose a z value guaranteed to be beyond our limit
-    float z0 = pow((1.0 / limit), (1.0 / 2.25));
-    psAssert (isfinite(z0), "fix this code: z0 should not be nan for %f", PAR[PM_PAR_7]);
-    float z1 = (1.0 / limit) / PAR[PM_PAR_7];
+    float z0 = 0.0;
+    float z1 = pow((1.0 / limit), (1.0 / ALPHA));
     psAssert (isfinite(z1), "fix this code: z1 should not be nan for %f", PAR[PM_PAR_7]);
-    z1 = PS_MAX (z0, z1);
-    z0 = 0.0;
-
-    // perform a type of bisection to find the value
-    float f0 = 1.0 / (1 + PAR[PM_PAR_7]*z0 + pow(z0, 2.25));
-    float f1 = 1.0 / (1 + PAR[PM_PAR_7]*z1 + pow(z1, 2.25));
-    while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
-        z = 0.5*(z0 + z1);
-        f = 1.0 / (1 + PAR[PM_PAR_7]*z + pow(z, 2.25));
-        if (f > limit) {
-            z0 = z;
-            f0 = f;
-        } else {
-            z1 = z;
-            f1 = f;
-        }
-        Nstep ++;
+    if (PAR[PM_PAR_7] < 0.0) z1 *= 2.0;
+
+    // starting guess:
+    z = 0.5*(z0 + z1);
+    float dz = 1.0;
+
+    // use Newton-Raphson to minimize f(z) - limit = 0
+    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
+	float q = (1.0 + PAR[PM_PAR_7]*z + pow(z, ALPHA));
+	float dqdz = (PAR[PM_PAR_7] + ALPHA*pow(z, ALPHA - 1.0));
+
+	float f = 1.0 / q;
+	float dfdz = -dqdz * f / q;
+
+	dz = (f - limit) / dfdz;
+
+	// fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
+	z -= dz;
     }
     psF64 radius = sigma * sqrt (2.0 * z);
@@ -475,2 +478,4 @@
 # undef PM_MODEL_FIT_STATUS
 # undef PM_MODEL_SET_LIMITS
+# undef ALPHA
+# undef ALPHA_M
Index: trunk/psModules/src/objects/models/pmModel_RGAUSS.c
===================================================================
--- trunk/psModules/src/objects/models/pmModel_RGAUSS.c	(revision 27530)
+++ trunk/psModules/src/objects/models/pmModel_RGAUSS.c	(revision 27531)
@@ -257,5 +257,5 @@
 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux)
 {
-    psF64 z, f;
+    psF64 z;
     int Nstep = 0;
     psEllipseShape shape;
@@ -291,19 +291,22 @@
     z0 = 0.0;
 
-    // perform a type of bisection to find the value
-    float f0 = 1.0 / (1 + z0 + pow(z0, PAR[PM_PAR_7]));
-    float f1 = 1.0 / (1 + z1 + pow(z1, PAR[PM_PAR_7]));
-    while ((Nstep < 10) && (fabs(z1 - z0) > 0.5)) {
-        z = 0.5*(z0 + z1);
-        f = 1.0 / (1 + z + pow(z, PAR[PM_PAR_7]));
-        if (f > limit) {
-            z0 = z;
-            f0 = f;
-        } else {
-            z1 = z;
-            f1 = f;
-        }
-        Nstep ++;
-    }
+    // starting guess:
+    z = 0.5*(z0 + z1);
+    float dz = 1.0;
+
+    for (int i = 0; (i < 10) && (fabs(dz) > 0.0001); i++) {
+	// use Newton-Raphson to minimize f(z) - limit = 0
+	float q = (1 + z + pow(z,PAR[PM_PAR_7]));
+	float dqdz = (1.0 + PAR[PM_PAR_7]*pow(z,PAR[PM_PAR_7] - 1.0));
+
+	float f = 1.0 / q;
+	float dfdz = -dqdz * f / q;
+
+	dz = (f - limit) / dfdz;
+
+	// fprintf (stderr, "%f %f %f : %f %f\n", f, z, dz, dfdz, q);
+	z -= dz;
+    }
+
     psF64 radius = sigma * sqrt (2.0 * z);
 
Index: trunk/psModules/src/objects/pmGrowthCurveGenerate.c
===================================================================
--- trunk/psModules/src/objects/pmGrowthCurveGenerate.c	(revision 27530)
+++ trunk/psModules/src/objects/pmGrowthCurveGenerate.c	(revision 27531)
@@ -157,5 +157,5 @@
 
     // measure the fitMag for this model
-    pmSourcePhotometryModel (&fitMag, model);
+    pmSourcePhotometryModel (&fitMag, NULL, model);
     growth->fitMag = fitMag;
 
Index: trunk/psModules/src/objects/pmPSF_IO.c
===================================================================
--- trunk/psModules/src/objects/pmPSF_IO.c	(revision 27530)
+++ trunk/psModules/src/objects/pmPSF_IO.c	(revision 27531)
@@ -491,30 +491,31 @@
         psMetadataAddS32 (header, PS_LIST_TAIL, "YCENTER", 0, "", psf->residuals->yCenter);
 
-        // write the residuals as three planes of the image
-        // this call creates an extension with NAXIS3 = 3
+        // write the residuals as planes of the image
+	psArray *images = psArrayAllocEmpty (1);
+	psArrayAdd (images, 1, psf->residuals->Ro);  // z = 0 is Ro
+
         if (psf->residuals->Rx) {
-            // this call creates an extension with NAXIS3 = 3
-            psArray *images = psArrayAllocEmpty (3);
-            psArrayAdd (images, 1, psf->residuals->Ro);
             psArrayAdd (images, 1, psf->residuals->Rx);
             psArrayAdd (images, 1, psf->residuals->Ry);
-
-            if (!psFitsWriteImageCube (file->fits, header, images, residName)) {
-                psError(psErrorCodeLast(), false, "Unable to write PSF residuals.");
-                psFree(images);
-                psFree(residName);
-                psFree(header);
-                return false;
-            }
-            psFree (images);
-        } else {
-            // this call creates an extension with NAXIS3 = 1
-            if (!psFitsWriteImage(file->fits, header, psf->residuals->Ro, 0, residName)) {
-                psError(psErrorCodeLast(), false, "Unable to write PSF residuals.");
-                psFree(residName);
-                psFree(header);
-                return false;
-            }
-        }
+	}
+
+	// note that all N plane are implicitly of the same type, so we convert the mask
+	if (psf->residuals->mask) {
+	    psImage *mask = psImageCopy (NULL, psf->residuals->mask, psf->residuals->Ro->type.type);
+	    psArrayAdd (images, 1, mask);
+	    psFree (mask);
+	}
+
+	// psFitsWriteImageCube (file->fits, header, images, residName);
+	// psFree (images);
+
+	if (!psFitsWriteImageCube (file->fits, header, images, residName)) {
+            psError(psErrorCodeLast(), false, "Unable to write PSF residuals.");
+            psFree(images);
+            psFree(residName);
+            psFree(header);
+            return false;
+        }
+        psFree (images);
         psFree (residName);
         psFree (header);
@@ -1017,6 +1018,21 @@
             return false;
         }
-        if (Nz > 1) {
-            assert (Nz == 3);
+
+	// note that all N plane are implicitly of the same type, so we convert the mask
+	psImage *mask = psImageCopy(NULL, psf->residuals->mask, psf->residuals->Ro->type.type);
+	psImageInit (psf->residuals->mask, 0);
+	psImageInit (psf->residuals->Rx, 0.0);
+	psImageInit (psf->residuals->Ry, 0.0);
+	switch (Nz) {
+	  case 1: // Ro only
+	    break;
+	  case 2: // Ro and mask
+            if (!psFitsReadImageBuffer(mask, file->fits, fullImage, 1)) {
+                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
+                return false;
+            }
+	    psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK);
+	    break;
+	  case 3: // Ro, Rx and Ry, no mask
             if (!psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1)) {
                 psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
@@ -1027,6 +1043,22 @@
                 return false;
             }
-        }
-        // XXX notice that we are not saving the resid->mask
+	    break;
+	  case 4: // Ro, Rx, Ry, and mask:
+            if (!psFitsReadImageBuffer(psf->residuals->Rx, file->fits, fullImage, 1)) {
+                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
+                return false;
+            }
+            if (!psFitsReadImageBuffer(psf->residuals->Ry, file->fits, fullImage, 2)) {
+                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
+                return false;
+            }
+            if (!psFitsReadImageBuffer(mask, file->fits, fullImage, 3)) {
+                psError(psErrorCodeLast(), false, "Unable to read PSF residual image.");
+                return false;
+            }
+	    psImageCopy (psf->residuals->mask, mask, PM_TYPE_RESID_MASK);
+	    break;
+        }
+	psFree (mask);
     }
 
Index: trunk/psModules/src/objects/pmSource.c
===================================================================
--- trunk/psModules/src/objects/pmSource.c	(revision 27530)
+++ trunk/psModules/src/objects/pmSource.c	(revision 27531)
@@ -46,5 +46,5 @@
     psFree(tmp->maskView);
     psFree(tmp->modelFlux);
-    psFree(tmp->psfFlux);
+    psFree(tmp->psfImage);
     psFree(tmp->moments);
     psFree(tmp->modelPSF);
@@ -52,4 +52,6 @@
     psFree(tmp->modelFits);
     psFree(tmp->extpars);
+    psFree(tmp->moments);
+    psFree(tmp->diffStats);
     psFree(tmp->blends);
     psTrace("psModules.objects", 10, "---- end ----\n");
@@ -68,5 +70,5 @@
     psFree (source->maskView);
     psFree (source->modelFlux);
-    psFree (source->psfFlux);
+    psFree (source->psfImage);
 
     source->pixels = NULL;
@@ -75,5 +77,5 @@
     source->maskView = NULL;
     source->modelFlux = NULL;
-    source->psfFlux = NULL;
+    source->psfImage = NULL;
     return;
 }
@@ -103,5 +105,5 @@
     source->maskView = NULL;
     source->modelFlux = NULL;
-    source->psfFlux = NULL;
+    source->psfImage = NULL;
     source->moments = NULL;
     source->blends = NULL;
@@ -113,9 +115,13 @@
     source->tmpFlags = 0;
     source->extpars = NULL;
+    source->diffStats = NULL;
+
     source->region = psRegionSet(NAN, NAN, NAN, NAN);
     psMemSetDeallocator(source, (psFreeFunc) sourceFree);
 
     // default values are NAN
-    source->psfMag = NAN;
+    source->psfMag     = NAN;
+    source->psfFlux    = NAN;
+    source->psfFluxErr = NAN;
     source->extMag = NAN;
     source->errMag = NAN;
@@ -259,7 +265,7 @@
         mySource->modelFlux = NULL;
 
-        // drop the old psfFlux pixels and force the user to re-create
-        psFree (mySource->psfFlux);
-        mySource->psfFlux = NULL;
+        // drop the old psfImage pixels and force the user to re-create
+        psFree (mySource->psfImage);
+        mySource->psfImage = NULL;
     }
     return extend;
@@ -873,11 +879,11 @@
 
     // if we already have a cached image, re-use that memory
-    source->psfFlux = psImageCopy (source->psfFlux, source->pixels, PS_TYPE_F32);
-    psImageInit (source->psfFlux, 0.0);
+    source->psfImage = psImageCopy (source->psfImage, source->pixels, PS_TYPE_F32);
+    psImageInit (source->psfImage, 0.0);
 
     // in some places (psphotEnsemble), we need a normalized version
     // in others, we just want the model.  which is more commonly used?
-    // psfFlux always has unity normalization (I0 = 1.0)
-    pmModelAdd (source->psfFlux, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
+    // psfImage always has unity normalization (I0 = 1.0)
+    pmModelAdd (source->psfImage, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM, maskVal);
     return true;
 }
Index: trunk/psModules/src/objects/pmSource.h
===================================================================
--- trunk/psModules/src/objects/pmSource.h	(revision 27530)
+++ trunk/psModules/src/objects/pmSource.h	(revision 27531)
@@ -16,4 +16,5 @@
 #include "pmMoments.h"
 #include "pmSourceExtendedPars.h"
+#include "pmSourceDiffStats.h"
 
 /// @addtogroup Objects Object Detection / Analysis Functions
@@ -38,7 +39,8 @@
 
 typedef enum {
-    PM_SOURCE_TMPF_MODEL_GUESS   = 0x0001,
-    PM_SOURCE_TMPF_SUBTRACTED    = 0x0002,
-    PM_SOURCE_TMPF_SIZE_MEASURED = 0x0004,
+    PM_SOURCE_TMPF_MODEL_GUESS       = 0x0001,
+    PM_SOURCE_TMPF_SUBTRACTED        = 0x0002,
+    PM_SOURCE_TMPF_SIZE_MEASURED     = 0x0004,
+    PM_SOURCE_TMPF_SIZE_CR_CANDIDATE = 0x0008,
 } pmSourceTmpF;
 
@@ -64,5 +66,5 @@
     psImage *maskView;                  ///< view into global image mask for this object region
     psImage *modelFlux;                 ///< cached copy of the best model for this source
-    psImage *psfFlux;                   ///< cached copy of the psf model for this source
+    psImage *psfImage;                   ///< cached copy of the psf model for this source
     pmMoments *moments;                 ///< Basic moments measured for the object.
     pmModel *modelPSF;                  ///< PSF Model fit (parameters and type)
@@ -74,4 +76,6 @@
     psArray *blends;                    ///< collection of sources thought to be confused with object
     float psfMag;                       ///< calculated from flux in modelPSF
+    float psfFlux;                      ///< calculated from flux in modelPSF
+    float psfFluxErr;                   ///< calculated from flux in modelPSF
     float extMag;                       ///< calculated from flux in modelEXT
     float errMag;                       ///< error in psfMag OR extMag (depending on type)
@@ -85,4 +89,5 @@
     psRegion region;                    ///< area on image covered by selected pixels
     pmSourceExtendedPars *extpars;      ///< extended source parameters
+    pmSourceDiffStats *diffStats;       ///< extra parameters for difference detections
 };
 
Index: trunk/psModules/src/objects/pmSourceDiffStats.c
===================================================================
--- trunk/psModules/src/objects/pmSourceDiffStats.c	(revision 27531)
+++ trunk/psModules/src/objects/pmSourceDiffStats.c	(revision 27531)
@@ -0,0 +1,42 @@
+/** @file  pmSourceDiffStats.c
+ *
+ *  Functions defining the pmSourceDiffStats structure and associated measurements
+
+ *  @author EAM, IfA
+ *
+ *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-10-03 20:59:16 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <pslib.h>
+#include "pmSourceDiffStats.h"
+
+/** initialization function for a pmSourceDiffStats object
+ **/
+void pmSourceDiffStatsInit(pmSourceDiffStats *diffStats)
+{
+    diffStats->fRatio = NAN;
+    diffStats->nRatioBad = NAN;
+    diffStats->nRatioMask = NAN;
+    diffStats->nRatioAll = NAN;
+    diffStats->nGood = 0;
+}
+
+/******************************************************************************
+pmSourceDiffStatsAlloc(): Allocate the pmSourceDiffStats structure and initialize the members
+to zero.
+*****************************************************************************/
+pmSourceDiffStats *pmSourceDiffStatsAlloc(void)
+{
+    pmSourceDiffStats *tmp = (pmSourceDiffStats *) psAlloc(sizeof(pmSourceDiffStats));
+    pmSourceDiffStatsInit(tmp);
+    return(tmp);
+}
Index: trunk/psModules/src/objects/pmSourceDiffStats.h
===================================================================
--- trunk/psModules/src/objects/pmSourceDiffStats.h	(revision 27531)
+++ trunk/psModules/src/objects/pmSourceDiffStats.h	(revision 27531)
@@ -0,0 +1,44 @@
+/* @file  pmSourceDiffStats.h
+ *
+ * @author EAM, IfA
+ *
+ * @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2009-02-16 22:30:50 $
+ * Copyright 2004 Institute for Astronomy, University of Hawaii
+ */
+
+# ifndef PM_SOURCE_DIFF_STATS_H
+# define PM_SOURCE_DIFF_STATS_H
+
+/// @addtogroup Objects Object Detection / Analysis Functions
+/// @{
+
+/* The elements of this structure are used to characterize diff image detections.
+ * The values in the structure are derived from the following measurements:
+ * fGood = \sum(flux) for pixels with S/N > +SN_LIMIT
+ * fBad  = \sum(flux) for pixels with S/N < -SN_LIMIT
+ * nGood = \sum(pixels) with S/N > +SN_LIMIT
+ * nBad  = \sum(pixels) with S/N < -SN_LIMIT
+ * nMask = \sum(pixels) masked
+ */
+
+typedef struct {
+    float fRatio;			// = fGood / (fGood + fBad)
+    float nRatioBad;			// = nGood / (nGood + nBad)
+    float nRatioMask;			// = nGood / (nGood + nMask)
+    float nRatioAll;			// = nGood / (nGood + nMask + nBad)
+    int   nGood;			// nGood as defined above
+} pmSourceDiffStats;
+
+
+/** pmSourceDiffStatsAlloc()
+ */
+pmSourceDiffStats *pmSourceDiffStatsAlloc(void);
+
+/** pmSourceDiffStatsInit()
+ * function to initialize the values in the structure
+ */
+void pmSourceDiffStatsInit(pmSourceDiffStats *diffStats);
+
+/// @}
+# endif
Index: trunk/psModules/src/objects/pmSourceIO.c
===================================================================
--- trunk/psModules/src/objects/pmSourceIO.c	(revision 27530)
+++ trunk/psModules/src/objects/pmSourceIO.c	(revision 27531)
@@ -539,4 +539,7 @@
                 status &= pmSourcesWrite_CMF_PS1_V2 (file->fits, readout, sources, file->header, outhead, dataname);
             }
+            if (!strcmp (exttype, "PS1_DV1")) {
+                status &= pmSourcesWrite_CMF_PS1_DV1 (file->fits, readout, sources, file->header, outhead, dataname);
+            }
 
             if (xsrcname) {
@@ -553,4 +556,7 @@
 		    status &= pmSourcesWrite_CMF_PS1_V2_XSRC (file->fits, sources, xsrcname, recipe);
 		}
+		if (!strcmp (exttype, "PS1_DV1")) {
+		    status &= pmSourcesWrite_CMF_PS1_DV1_XSRC (file->fits, sources, xsrcname, recipe);
+		}
             }
             if (xfitname) {
@@ -567,4 +573,7 @@
 		    status &= pmSourcesWrite_CMF_PS1_V2_XFIT (file->fits, sources, xfitname);
 		}
+		if (!strcmp (exttype, "PS1_DV1")) {
+		    status &= pmSourcesWrite_CMF_PS1_DV1_XFIT (file->fits, sources, xfitname);
+		}
             }
 	    psFree (outhead);
@@ -1015,4 +1024,7 @@
                 sources = pmSourcesRead_CMF_PS1_V2 (file->fits, hdu->header);
             }
+            if (!strcmp (exttype, "PS1_DV1")) {
+                sources = pmSourcesRead_CMF_PS1_DV1 (file->fits, hdu->header);
+            }
 
             if (!pmReadoutReadDetEff(file->fits, readout, deteffname)) {
Index: trunk/psModules/src/objects/pmSourceIO.h
===================================================================
--- trunk/psModules/src/objects/pmSourceIO.h	(revision 27530)
+++ trunk/psModules/src/objects/pmSourceIO.h	(revision 27531)
@@ -43,4 +43,8 @@
 bool pmSourcesWrite_CMF_PS1_V2_XFIT (psFits *fits, psArray *sources, char *extname);
 
+bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources, psMetadata *imageHeader, psMetadata *tableHeader, char *extname);
+bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe);
+bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, psArray *sources, char *extname);
+
 bool pmSource_CMF_WritePHU (const pmFPAview *view, pmFPAfile *file, pmConfig *config);
 
@@ -53,4 +57,5 @@
 psArray *pmSourcesRead_CMF_PS1_V1 (psFits *fits, psMetadata *header);
 psArray *pmSourcesRead_CMF_PS1_V2 (psFits *fits, psMetadata *header);
+psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header);
 
 bool pmSourcesWritePSFs (psArray *sources, char *filename);
Index: trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c
===================================================================
--- trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c	(revision 27531)
+++ trunk/psModules/src/objects/pmSourceIO_CMF_PS1_DV1.c	(revision 27531)
@@ -0,0 +1,660 @@
+/** @file  pmSourceIO.c
+ *
+ *  @author EAM, IfA
+ *
+ *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2009-02-18 02:44:19 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <pslib.h>
+
+#include "pmConfig.h"
+#include "pmDetrendDB.h"
+
+#include "pmHDU.h"
+#include "pmFPA.h"
+#include "pmFPALevel.h"
+#include "pmFPAview.h"
+#include "pmFPAfile.h"
+
+#include "pmSpan.h"
+#include "pmFootprint.h"
+#include "pmPeaks.h"
+#include "pmMoments.h"
+#include "pmGrowthCurve.h"
+#include "pmResiduals.h"
+#include "pmTrend2D.h"
+#include "pmPSF.h"
+#include "pmModel.h"
+#include "pmSource.h"
+#include "pmModelClass.h"
+#include "pmSourceIO.h"
+
+// panstarrs-style FITS table output (header + table in 1st extension)
+// this format consists of a header derived from the image header
+// followed by a zero-size matrix, followed by the table data
+
+// This version has elements intended for difference images & forced photometry:
+// diffStats entries (good for dipoles); flux + flux error (for insignificant detections)
+bool pmSourcesWrite_CMF_PS1_DV1 (psFits *fits, pmReadout *readout, psArray *sources,
+                                psMetadata *imageHeader, psMetadata *tableHeader, char *extname)
+{
+    PS_ASSERT_PTR_NON_NULL(fits, false);
+    PS_ASSERT_PTR_NON_NULL(sources, false);
+    PS_ASSERT_PTR_NON_NULL(extname, false);
+
+    int i;
+    psArray *table;
+    psMetadata *row;
+    psF32 *PAR, *dPAR;
+    psEllipseAxes axes;
+    psF32 xPos, yPos;
+    psF32 xErr, yErr;
+    psF32 errMag, chisq, apRadius;
+    psS32 nPix, nDOF;
+
+    pmChip *chip = readout->parent->parent;
+    pmFPA  *fpa  = chip->parent;
+
+    bool status1 = false;
+    bool status2 = false;
+    float magOffset = NAN;
+    float exptime   = psMetadataLookupF32 (&status1, fpa->concepts, "FPA.EXPOSURE");
+    float zeropt    = psMetadataLookupF32(&status2, fpa->concepts, "FPA.ZP");
+    if (!isfinite(zeropt)) {
+        zeropt    = psMetadataLookupF32 (&status2, imageHeader, "ZPT_OBS");
+    }
+    if (status1 && status2 && (exptime > 0.0)) {
+        magOffset = zeropt + 2.5*log10(exptime);
+    }
+    float zeroptErr = psMetadataLookupF32 (&status2, imageHeader, "ZPT_ERR");
+
+    // if the sequence is defined, write these in seq order; otherwise
+    // write them in S/N order:
+    if (sources->n > 0) {
+        pmSource *source = (pmSource *) sources->data[0];
+        if (source->seq == -1) {
+          // let's write these out in S/N order
+          sources = psArraySort (sources, pmSourceSortBySN);
+        } else {
+          sources = psArraySort (sources, pmSourceSortBySeq);
+        }
+    }
+
+    table = psArrayAllocEmpty (sources->n);
+
+    // we write out PSF-fits for all sources, regardless of quality.  the source flags tell us the state
+    // by the time we call this function, all values should be assigned.  let's use asserts to be sure in some cases.
+    for (i = 0; i < sources->n; i++) {
+        pmSource *source = (pmSource *) sources->data[i];
+
+        // If source->seq is -1, source was generated in this analysis.  If source->seq is
+        // not -1, source was read from elsewhere: in the latter case, preserve the source
+        // ID.  source.seq is used instead of source.id since the latter is a const
+        // generated on Alloc, and would thus be wrong for read in sources.
+        if (source->seq == -1) {
+          source->seq = i;
+        }
+
+        // no difference between PSF and non-PSF model
+        pmModel *model = source->modelPSF;
+
+        if (model != NULL) {
+            PAR = model->params->data.F32;
+            dPAR = model->dparams->data.F32;
+            xPos = PAR[PM_PAR_XPOS];
+            yPos = PAR[PM_PAR_YPOS];
+            if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
+              xErr = dPAR[PM_PAR_XPOS];
+              yErr = dPAR[PM_PAR_YPOS];
+            } else {
+              // in linear-fit mode, there is no error on the centroid
+              xErr = source->peak->dx;
+              yErr = source->peak->dy;
+            }
+            if (isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX]) && isfinite(PAR[PM_PAR_SXX])) {
+                axes = pmPSF_ModelToAxes (PAR, 20.0);
+            } else {
+                axes.major = NAN;
+                axes.minor = NAN;
+                axes.theta = NAN;
+            }
+            chisq = model->chisq;
+            nDOF = model->nDOF;
+            nPix = model->nPix;
+            apRadius = source->apRadius;
+            errMag = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
+        } else {
+            xPos = source->peak->xf;
+            yPos = source->peak->yf;
+            xErr = source->peak->dx;
+            yErr = source->peak->dy;
+            axes.major = NAN;
+            axes.minor = NAN;
+            axes.theta = NAN;
+            chisq = NAN;
+            nDOF = 0;
+            nPix = 0;
+            apRadius = NAN;
+            errMag = NAN;
+        }
+
+        float calMag = isfinite(magOffset) ? source->psfMag + magOffset : NAN;
+        float peakMag = (source->peak->flux > 0) ? -2.5*log10(source->peak->flux) : NAN;
+        psS16 nImageOverlap = 1;
+
+        psSphere ptSky = {0.0, 0.0, 0.0, 0.0};
+        float posAngle = 0.0;
+        float pltScale = 0.0;
+        pmSourceLocalAstrometry (&ptSky, &posAngle, &pltScale, chip, xPos, yPos);
+
+	pmSourceDiffStats diffStats;
+	pmSourceDiffStatsInit(&diffStats);
+	if (source->diffStats) {
+	    diffStats = *source->diffStats;
+	}
+
+        row = psMetadataAlloc ();
+        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
+        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF",            PS_DATA_F32, "PSF x coordinate",                           xPos);
+        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF",            PS_DATA_F32, "PSF y coordinate",                           yPos);
+        psMetadataAdd (row, PS_LIST_TAIL, "X_PSF_SIG",        PS_DATA_F32, "Sigma in PSF x coordinate",                  xErr); // XXX this is only measured for non-linear fits
+        psMetadataAdd (row, PS_LIST_TAIL, "Y_PSF_SIG",        PS_DATA_F32, "Sigma in PSF y coordinate",                  yErr); // XXX this is only measured for non-linear fits
+        psMetadataAdd (row, PS_LIST_TAIL, "POSANGLE",         PS_DATA_F32, "position angle at source (degrees)",         posAngle*PS_DEG_RAD);
+        psMetadataAdd (row, PS_LIST_TAIL, "PLTSCALE",         PS_DATA_F32, "plate scale at source (arcsec/pixel)",       pltScale*PS_DEG_RAD*3600.0);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG",     PS_DATA_F32, "PSF fit instrumental magnitude",             source->psfMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_MAG_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",        errMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX",     PS_DATA_F32, "PSF fit instrumental magnitude",            source->psfFlux);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_INST_FLUX_SIG", PS_DATA_F32, "Sigma of PSF instrumental magnitude",       source->psfFluxErr);
+        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG",           PS_DATA_F32, "magnitude in standard aperture",             source->apMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "AP_MAG_RADIUS",    PS_DATA_F32, "radius used for aperture mags",              apRadius);
+        psMetadataAdd (row, PS_LIST_TAIL, "PEAK_FLUX_AS_MAG", PS_DATA_F32, "Peak flux expressed as magnitude",           peakMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG",      PS_DATA_F32, "PSF Magnitude using supplied calibration",   calMag);
+        psMetadataAdd (row, PS_LIST_TAIL, "CAL_PSF_MAG_SIG",  PS_DATA_F32, "measured scatter of zero point calibration", zeroptErr);
+        psMetadataAdd (row, PS_LIST_TAIL, "RA_PSF",           PS_DATA_F64, "PSF RA coordinate (degrees)",                ptSky.r*PS_DEG_RAD);
+        psMetadataAdd (row, PS_LIST_TAIL, "DEC_PSF",          PS_DATA_F64, "PSF DEC coordinate (degrees)",               ptSky.d*PS_DEG_RAD);
+        psMetadataAdd (row, PS_LIST_TAIL, "SKY",              PS_DATA_F32, "Sky level",                                  source->sky);
+        psMetadataAdd (row, PS_LIST_TAIL, "SKY_SIGMA",        PS_DATA_F32, "Sigma of sky level",                         source->skyErr);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_CHISQ",        PS_DATA_F32, "Chisq of PSF-fit",                           chisq);
+        psMetadataAdd (row, PS_LIST_TAIL, "CR_NSIGMA",        PS_DATA_F32, "Nsigma deviations from PSF to CF",           source->crNsigma);
+        psMetadataAdd (row, PS_LIST_TAIL, "EXT_NSIGMA",       PS_DATA_F32, "Nsigma deviations from PSF to EXT",          source->extNsigma);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MAJOR",        PS_DATA_F32, "PSF width (major axis)",                     axes.major);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_MINOR",        PS_DATA_F32, "PSF width (minor axis)",                     axes.minor);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_THETA",        PS_DATA_F32, "PSF orientation angle",                      axes.theta);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_QF",           PS_DATA_F32, "PSF coverage/quality factor",                source->pixWeight);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NDOF",         PS_DATA_S32, "degrees of freedom",                         nDOF);
+        psMetadataAdd (row, PS_LIST_TAIL, "PSF_NPIX",         PS_DATA_S32, "number of pixels in fit",                    nPix);
+
+        // distinguish moments measure from window vs S/N > XX ??
+        float mxx = source->moments ? source->moments->Mxx : NAN;
+        float mxy = source->moments ? source->moments->Mxy : NAN;
+        float myy = source->moments ? source->moments->Myy : NAN;
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XX",       PS_DATA_F32, "second moments (X^2)",                      mxx);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_XY",       PS_DATA_F32, "second moments (X*Y)",                      mxy);
+        psMetadataAdd (row, PS_LIST_TAIL, "MOMENTS_YY",       PS_DATA_F32, "second moments (Y*Y)",                      myy);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NPOS",        PS_DATA_S32, "nPos (n pix > 3 sigma)",                    diffStats.nGood);
+        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_FRATIO",      PS_DATA_F32, "fPos / (fPos + fNeg)",                      diffStats.fRatio);
+        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_BAD",  PS_DATA_F32, "nPos / (nPos + nNeg)",                      diffStats.nRatioBad);
+        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_MASK", PS_DATA_F32, "nPos / (nPos + nMask)",                     diffStats.nRatioMask);
+        psMetadataAdd (row, PS_LIST_TAIL, "DIFF_NRATIO_ALL",  PS_DATA_F32, "nPos / (nGood + nMask + nBad)",             diffStats.nRatioAll);
+
+        psMetadataAdd (row, PS_LIST_TAIL, "FLAGS",            PS_DATA_U32, "psphot analysis flags",                      source->mode);
+
+        // XXX not sure how to get this : need to load Nimages with weight?
+        psMetadataAdd (row, PS_LIST_TAIL, "N_FRAMES",         PS_DATA_U16, "Number of frames overlapping source center", nImageOverlap);
+        psMetadataAdd (row, PS_LIST_TAIL, "PADDING",          PS_DATA_S16, "padding", 0);
+
+        psArrayAdd (table, 100, row);
+        psFree (row);
+
+        // EXT_NSIGMA will be NAN if: 1) contour ellipse is imaginary; 2) source is not
+        // subtracted
+
+        // CR_NSIGMA will be NAN if: 1) source is not subtracted; 2) source is on the image
+        // edge; 3) any pixels in the 3x3 peak region are masked;
+    }
+
+    psMetadata *header = psMetadataCopy(NULL, tableHeader);
+    pmSourceMasksHeader(header);
+
+    if (table->n == 0) {
+        psFitsWriteBlank(fits, header, extname);
+        psFree(table);
+        psFree(header);
+        return true;
+    }
+
+    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
+    if (!psFitsWriteTable(fits, header, table, extname)) {
+        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
+        psFree(table);
+        psFree(header);
+        return false;
+    }
+    psFree(table);
+    psFree(header);
+
+    return true;
+}
+
+// read in a readout from the fits file
+psArray *pmSourcesRead_CMF_PS1_DV1 (psFits *fits, psMetadata *header)
+{
+    PS_ASSERT_PTR_NON_NULL(fits, false);
+    PS_ASSERT_PTR_NON_NULL(header, false);
+
+    bool status;
+    psF32 *PAR, *dPAR;
+    psEllipseAxes axes;
+
+    // define PSF model type
+    // XXX need to carry the extra model parameters
+    int modelType = pmModelClassGetType ("PS_MODEL_GAUSS");
+
+    char *PSF_NAME = psMetadataLookupStr (&status, header, "PSF_NAME");
+    if (PSF_NAME != NULL) {
+        modelType = pmModelClassGetType (PSF_NAME);
+    }
+    assert (modelType > -1);
+
+    // We get the size of the table, and allocate the array of sources first because the table
+    // is large and ephemeral --- when the table gets blown away, whatever is allocated after
+    // the table is read blocks the free.  In fact, it's better to read the table row by row.
+    long numSources = psFitsTableSize(fits); // Number of sources in table
+    psArray *sources = psArrayAlloc(numSources); // Array of sources, to return
+
+    // convert the table to the pmSource entriesa
+    for (int i = 0; i < numSources; i++) {
+        psMetadata *row = psFitsReadTableRow(fits, i); // Table row
+
+        pmSource *source = pmSourceAlloc ();
+        pmModel *model = pmModelAlloc (modelType);
+        source->modelPSF  = model;
+        source->type = PM_SOURCE_TYPE_STAR; // XXX this should be added to the flags
+
+        // NOTE: A SEGV here because "model" is NULL is probably caused by not initialising the models.
+        PAR = model->params->data.F32;
+        dPAR = model->dparams->data.F32;
+
+        source->seq       = psMetadataLookupU32 (&status, row, "IPP_IDET");
+        PAR[PM_PAR_XPOS]  = psMetadataLookupF32 (&status, row, "X_PSF");
+        PAR[PM_PAR_YPOS]  = psMetadataLookupF32 (&status, row, "Y_PSF");
+        dPAR[PM_PAR_XPOS] = psMetadataLookupF32 (&status, row, "X_PSF_SIG");
+        dPAR[PM_PAR_YPOS] = psMetadataLookupF32 (&status, row, "Y_PSF_SIG");
+        axes.major        = psMetadataLookupF32 (&status, row, "PSF_MAJOR");
+        axes.minor        = psMetadataLookupF32 (&status, row, "PSF_MINOR");
+        axes.theta        = psMetadataLookupF32 (&status, row, "PSF_THETA");
+
+        PAR[PM_PAR_SKY]   = psMetadataLookupF32 (&status, row, "SKY");
+        dPAR[PM_PAR_SKY]  = psMetadataLookupF32 (&status, row, "SKY_SIGMA");
+        source->sky       = PAR[PM_PAR_SKY];
+        source->skyErr    = dPAR[PM_PAR_SKY];
+
+        // XXX use these to determine PAR[PM_PAR_I0]?
+        source->psfFlux    = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX");
+        source->psfFluxErr = psMetadataLookupF32 (&status, row, "PSF_INST_FLUX_SIG");
+
+        source->psfMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG");
+        source->errMag    = psMetadataLookupF32 (&status, row, "PSF_INST_MAG_SIG");
+        source->apMag     = psMetadataLookupF32 (&status, row, "AP_MAG");
+
+        // XXX this scaling is incorrect: does not include the 2 \pi AREA factor
+        PAR[PM_PAR_I0]    = (isfinite(source->psfMag)) ? pow(10.0, -0.4*source->psfMag) : NAN;
+        dPAR[PM_PAR_I0]   = (isfinite(source->psfMag)) ? PAR[PM_PAR_I0] * source->errMag : NAN;
+
+        pmPSF_AxesToModel (PAR, axes);
+
+        float peakMag     = psMetadataLookupF32 (&status, row, "PEAK_FLUX_AS_MAG");
+        float peakFlux    = (isfinite(peakMag)) ? pow(10.0, -0.4*peakMag) : NAN;
+
+        // recreate the peak to match (xPos, yPos) +/- (xErr, yErr)
+        source->peak = pmPeakAlloc(PAR[PM_PAR_XPOS], PAR[PM_PAR_YPOS], peakFlux, PM_PEAK_LONE);
+        source->peak->flux = peakFlux;
+        source->peak->dx   = dPAR[PM_PAR_XPOS];
+        source->peak->dy   = dPAR[PM_PAR_YPOS];
+        source->peak->SN   = sqrt(source->peak->flux); // XXX a proxy: various functions sort by peak S/N
+
+        source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");
+        source->crNsigma  = psMetadataLookupF32 (&status, row, "CR_NSIGMA");
+        source->extNsigma = psMetadataLookupF32 (&status, row, "EXT_NSIGMA");
+        source->apRadius  = psMetadataLookupS32 (&status, row, "AP_MAG_RADIUS");
+
+        // note that some older versions used PSF_PROBABILITY: this was not well defined.
+        model->chisq      = psMetadataLookupF32 (&status, row, "PSF_CHISQ");
+        model->nDOF       = psMetadataLookupS32 (&status, row, "PSF_NDOF");
+        model->nPix       = psMetadataLookupS32 (&status, row, "PSF_NPIX");
+
+        source->moments = pmMomentsAlloc ();
+        source->moments->Mxx = psMetadataLookupF32 (&status, row, "MOMENTS_XX");
+        source->moments->Mxy = psMetadataLookupF32 (&status, row, "MOMENTS_XY");
+        source->moments->Myy = psMetadataLookupF32 (&status, row, "MOMENTS_YY");
+
+        source->mode = psMetadataLookupU32 (&status, row, "FLAGS");
+        assert (status);
+
+	int nPos = psMetadataLookupS32 (&status, row, "DIFF_NPOS");
+	if (nPos) {
+	    source->diffStats = pmSourceDiffStatsAlloc();
+	    source->diffStats->nGood      = nPos;
+	    source->diffStats->fRatio     = psMetadataLookupF32 (&status, row, "DIFF_FRATIO");
+	    source->diffStats->nRatioBad  = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_BAD");
+	    source->diffStats->nRatioMask = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_MASK");
+	    source->diffStats->nRatioAll  = psMetadataLookupF32 (&status, row, "DIFF_NRATIO_ALL");
+	}
+
+        sources->data[i] = source;
+        psFree(row);
+    }
+
+    return sources;
+}
+
+// XXX this layout is still the same as PS1_DEV_1
+bool pmSourcesWrite_CMF_PS1_DV1_XSRC (psFits *fits, psArray *sources, char *extname, psMetadata *recipe)
+{
+
+    bool status;
+    psArray *table;
+    psMetadata *row;
+    psF32 *PAR, *dPAR;
+    psF32 xPos, yPos;
+    psF32 xErr, yErr;
+
+    // create a header to hold the output data
+    psMetadata *outhead = psMetadataAlloc ();
+
+    // write the links to the image header
+    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname);
+
+    // let's write these out in S/N order
+    sources = psArraySort (sources, pmSourceSortBySN);
+
+    table = psArrayAllocEmpty (sources->n);
+
+    // which extended source analyses should we perform?
+    // bool doPetrosian    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_PETROSIAN");
+    // bool doIsophotal    = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ISOPHOTAL");
+    // bool doAnnuli       = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANNULI");
+    // bool doKron         = psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_KRON");
+
+    psVector *radialBinsLower = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.LOWER");
+    psVector *radialBinsUpper = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
+    assert (radialBinsLower->n == radialBinsUpper->n);
+
+    // we write out all sources, regardless of quality.  the source flags tell us the state
+    for (int i = 0; i < sources->n; i++) {
+        // skip source if it is not a ext sourc
+        // XXX we have two places that extended source parameters are measured:
+        // psphotExtendedSources, which measures the aperture-like parameters and (potentially) the psf-convolved extended source models,
+        // psphotFitEXT, which does the simple extended source model fit (not psf-convolved)
+        // should we require both?
+
+        pmSource *source = sources->data[i];
+
+        // skip sources without measurements
+        if (source->extpars == NULL) continue;
+
+        // we require a PSF model fit (ignore the real crud)
+        pmModel *model = source->modelPSF;
+        if (model == NULL) continue;
+
+        // XXX I need to split the extended models from the extended aperture measurements
+        PAR = model->params->data.F32;
+        dPAR = model->dparams->data.F32;
+        xPos = PAR[PM_PAR_XPOS];
+        yPos = PAR[PM_PAR_YPOS];
+        xErr = dPAR[PM_PAR_XPOS];
+        yErr = dPAR[PM_PAR_YPOS];
+
+        row = psMetadataAlloc ();
+
+        // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
+        psMetadataAdd (row, PS_LIST_TAIL, "IPP_IDET",         PS_DATA_U32, "IPP detection identifier index",             source->seq);
+        psMetadataAdd (row, PS_LIST_TAIL, "X_EXT",            PS_DATA_F32, "EXT model x coordinate",                     xPos);
+        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT",            PS_DATA_F32, "EXT model y coordinate",                     yPos);
+        psMetadataAdd (row, PS_LIST_TAIL, "X_EXT_SIG",        PS_DATA_F32, "Sigma in EXT x coordinate",                  xErr);
+        psMetadataAdd (row, PS_LIST_TAIL, "Y_EXT_SIG",        PS_DATA_F32, "Sigma in EXT y coordinate",                  yErr);
+
+# if (0)
+        // Petrosian measurements
+        // XXX insert header data: petrosian ref radius, flux ratio
+        if (doPetrosian) {
+            pmSourcePetrosianValues *petrosian = source->extpars->petrosian;
+            if (petrosian) {
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       petrosian->mag);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", petrosian->magErr);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          petrosian->rad);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    petrosian->radErr);
+            } else {
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG",        PS_DATA_F32, "Petrosian Magnitude",       NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_MAG_ERR",    PS_DATA_F32, "Petrosian Magnitude Error", NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS",     PS_DATA_F32, "Petrosian Radius",          NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "PETRO_RADIUS_ERR", PS_DATA_F32, "Petrosian Radius Error",    NAN);
+            }
+        }
+
+        // Kron measurements
+        if (doKron) {
+            pmSourceKronValues *kron = source->extpars->kron;
+            if (kron) {
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       kron->mag);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", kron->magErr);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          kron->rad);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    kron->radErr);
+            } else {
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG",        PS_DATA_F32, "Kron Magnitude",       NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_MAG_ERR",    PS_DATA_F32, "Kron Magnitude Error", NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS",     PS_DATA_F32, "Kron Radius",          NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "KRON_RADIUS_ERR", PS_DATA_F32, "Kron Radius Error",    NAN);
+            }
+        }
+
+        // Isophot measurements
+        // XXX insert header data: isophotal level
+        if (doIsophotal) {
+            pmSourceIsophotalValues *isophot = source->extpars->isophot;
+            if (isophot) {
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       isophot->mag);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", isophot->magErr);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          isophot->rad);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    isophot->radErr);
+            } else {
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG",        PS_DATA_F32, "Isophot Magnitude",       NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_MAG_ERR",    PS_DATA_F32, "Isophot Magnitude Error", NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS",     PS_DATA_F32, "Isophot Radius",          NAN);
+                psMetadataAdd (row, PS_LIST_TAIL, "ISOPHOT_RADIUS_ERR", PS_DATA_F32, "Isophot Radius Error",    NAN);
+            }
+        }
+
+        // Flux Annuli
+        if (doAnnuli) {
+            pmSourceAnnuli *annuli = source->extpars->annuli;
+            if (annuli) {
+                psVector *fluxVal = annuli->flux;
+                psVector *fluxErr = annuli->fluxErr;
+                psVector *fluxVar = annuli->fluxVar;
+
+                for (int j = 0; j < fluxVal->n; j++) {
+                    char name[32];
+                    sprintf (name, "FLUX_VAL_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", fluxVal->data.F32[j]);
+                    sprintf (name, "FLUX_ERR_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", fluxErr->data.F32[j]);
+                    sprintf (name, "FLUX_VAR_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", fluxVar->data.F32[j]);
+                }
+            } else {
+                for (int j = 0; j < radialBinsLower->n; j++) {
+                    char name[32];
+                    sprintf (name, "FLUX_VAL_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux value in annulus", NAN);
+                    sprintf (name, "FLUX_ERR_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux error in annulus", NAN);
+                    sprintf (name, "FLUX_VAR_R_%02d", j);
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "flux stdev in annulus", NAN);
+                }
+            }
+        }
+
+# endif
+        psArrayAdd (table, 100, row);
+        psFree (row);
+    }
+
+    if (table->n == 0) {
+        psFitsWriteBlank (fits, outhead, extname);
+        psFree (outhead);
+        psFree (table);
+        return true;
+    }
+
+    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
+    if (!psFitsWriteTable (fits, outhead, table, extname)) {
+        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
+        psFree (outhead);
+        psFree(table);
+        return false;
+    }
+    psFree (outhead);
+    psFree (table);
+
+    return true;
+}
+
+// XXX this layout is still the same as PS1_DEV_1
+bool pmSourcesWrite_CMF_PS1_DV1_XFIT (psFits *fits, psArray *sources, char *extname)
+{
+
+    psArray *table;
+    psMetadata *row;
+    psF32 *PAR, *dPAR;
+    psEllipseAxes axes;
+    psF32 xPos, yPos;
+    psF32 xErr, yErr;
+    char name[64];
+
+    // create a header to hold the output data
+    psMetadata *outhead = psMetadataAlloc ();
+
+    // write the links to the image header
+    psMetadataAddStr (outhead, PS_LIST_TAIL, "EXTNAME", PS_META_REPLACE, "xsrc table extension", extname);
+
+    // let's write these out in S/N order
+    sources = psArraySort (sources, pmSourceSortBySN);
+
+    // we are writing one row per model; we need to write out same number of columns for each row: find the max Nparams
+    int nParamMax = 0;
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+        if (source->modelFits == NULL) continue;
+        for (int j = 0; j < source->modelFits->n; j++) {
+            pmModel *model = source->modelFits->data[j];
+            assert (model);
+            nParamMax = PS_MAX (nParamMax, model->params->n);
+        }
+    }
+
+    table = psArrayAllocEmpty (sources->n);
+
+    // we write out all sources, regardless of quality.  the source flags tell us the state
+    for (int i = 0; i < sources->n; i++) {
+
+        pmSource *source = sources->data[i];
+
+        // XXX if no model fits are saved, write out modelEXT?
+        if (source->modelFits == NULL) continue;
+
+        // We have multiple sources : need to flag the one used to subtract the light (the 'best' model)
+        for (int j = 0; j < source->modelFits->n; j++) {
+
+            // choose the convolved EXT model, if available, otherwise the simple one
+            pmModel *model = source->modelFits->data[j];
+            assert (model);
+
+            PAR = model->params->data.F32;
+            dPAR = model->dparams->data.F32;
+            xPos = PAR[PM_PAR_XPOS];
+            yPos = PAR[PM_PAR_YPOS];
+            xErr = dPAR[PM_PAR_XPOS];
+            yErr = dPAR[PM_PAR_YPOS];
+
+            axes = pmPSF_ModelToAxes (PAR, 20.0);
+
+            row = psMetadataAlloc ();
+
+            // XXX we are not writing out the mode (flags) or the type (psf, ext, etc)
+            psMetadataAddU32 (row, PS_LIST_TAIL, "IPP_IDET",         0, "IPP detection identifier index",             source->seq);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT",            0, "EXT model x coordinate",                     xPos);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT",            0, "EXT model y coordinate",                     yPos);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "X_EXT_SIG",        0, "Sigma in EXT x coordinate",                  xErr);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "Y_EXT_SIG",        0, "Sigma in EXT y coordinate",                  yErr);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG",     0, "EXT fit instrumental magnitude",             model->mag);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_INST_MAG_SIG", 0, "Sigma of PSF instrumental magnitude",        model->magErr);
+
+            psMetadataAddF32 (row, PS_LIST_TAIL, "NPARAMS",          0, "number of model parameters",                 model->params->n);
+            psMetadataAddStr (row, PS_LIST_TAIL, "MODEL_TYPE",       0, "name of model",                              pmModelClassGetName (model->type));
+
+            // XXX these should be major and minor, not 'x' and 'y'
+            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MAJ",    0, "EXT width in x coordinate",                  axes.major);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_WIDTH_MIN",    0, "EXT width in y coordinate",                  axes.minor);
+            psMetadataAddF32 (row, PS_LIST_TAIL, "EXT_THETA",        0, "EXT orientation angle",                      axes.theta);
+
+            // write out the other generic parameters
+            for (int k = 0; k < nParamMax; k++) {
+                if (k == PM_PAR_I0) continue;
+                if (k == PM_PAR_SKY) continue;
+                if (k == PM_PAR_XPOS) continue;
+                if (k == PM_PAR_YPOS) continue;
+                if (k == PM_PAR_SXX) continue;
+                if (k == PM_PAR_SXY) continue;
+                if (k == PM_PAR_SYY) continue;
+
+                snprintf (name, 64, "EXT_PAR_%02d", k);
+
+                if (k < model->params->n) {
+                    psMetadataAdd (row, PS_LIST_TAIL, name, PS_DATA_F32, "", model->params->data.F32[k]);
+                } else {
+                    psMetadataAddF32 (row, PS_LIST_TAIL, name, PS_DATA_F32, "", NAN);
+                }
+            }
+
+            // XXX other parameters which may be set.
+            // XXX flag / value to define the model
+            // XXX write out the model type, fit status flags
+
+            psArrayAdd (table, 100, row);
+            psFree (row);
+        }
+    }
+
+    if (table->n == 0) {
+        psFitsWriteBlank (fits, outhead, extname);
+        psFree (outhead);
+        psFree (table);
+        return true;
+    }
+
+    psTrace ("pmFPAfile", 5, "writing ext data %s\n", extname);
+    if (!psFitsWriteTable (fits, outhead, table, extname)) {
+        psError(PS_ERR_IO, false, "writing ext data %s\n", extname);
+        psFree (outhead);
+        psFree(table);
+        return false;
+    }
+    psFree (outhead);
+    psFree (table);
+    return true;
+}
Index: trunk/psModules/src/objects/pmSourceMoments.c
===================================================================
--- trunk/psModules/src/objects/pmSourceMoments.c	(revision 27530)
+++ trunk/psModules/src/objects/pmSourceMoments.c	(revision 27531)
@@ -163,8 +163,10 @@
     }
 
-    // if we have less than (1/2) of the possible pixels, force a retry
+    // if we have less than (1/4) of the possible pixels (in circle or box), force a retry
+    int minPixels = PS_MIN(0.75*R2, source->pixels->numCols*source->pixels->numRows/4.0);
+
     // XXX EAM - the limit is a bit arbitrary.  make it user defined?
-    if ((numPixels < 0.75*R2) || (Sum <= 0)) {
-	psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, (int)(0.75*R2), Sum);
+    if ((numPixels < minPixels) || (Sum <= 0)) {
+	psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", numPixels, minPixels, Sum);
 	return (false);
     }
Index: trunk/psModules/src/objects/pmSourcePhotometry.c
===================================================================
--- trunk/psModules/src/objects/pmSourcePhotometry.c	(revision 27530)
+++ trunk/psModules/src/objects/pmSourcePhotometry.c	(revision 27531)
@@ -109,7 +109,10 @@
 	psAssert (isfinite(fluxScale), "how can the flux scale be invalid? source at %d, %d\n", source->peak->x, source->peak->y);
 	psAssert (fluxScale > 0.0, "how can the flux scale be negative? source at %d, %d\n", source->peak->x, source->peak->y);
-	source->psfMag = -2.5*log10(fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0]);
+	source->psfFlux = fluxScale * source->modelPSF->params->data.F32[PM_PAR_I0];
+	source->psfFluxErr = fluxScale * source->modelPSF->dparams->data.F32[PM_PAR_I0];
+	source->psfMag = -2.5*log10(source->psfFlux);
     } else {
-        status = pmSourcePhotometryModel (&source->psfMag, source->modelPSF);
+        status = pmSourcePhotometryModel (&source->psfMag, &source->psfFlux, source->modelPSF);
+	source->psfFluxErr = source->psfFlux * (source->modelPSF->dparams->data.F32[PM_PAR_I0] / source->modelPSF->params->data.F32[PM_PAR_I0]);
     }
 
@@ -119,5 +122,5 @@
 	for (int i = 0; i < source->modelFits->n; i++) {
 	    pmModel *model = source->modelFits->data[i];
-	    status = pmSourcePhotometryModel (&model->mag, model);
+	    status = pmSourcePhotometryModel (&model->mag, NULL, model);
 	    if (model == source->modelEXT) foundEXT = true;
 	}
@@ -125,9 +128,9 @@
 	    source->extMag = source->modelEXT->mag;
 	} else {
-	    status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
+	    status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
 	}
     } else {
 	if (source->modelEXT) {
-	    status = pmSourcePhotometryModel (&source->extMag, source->modelEXT);
+	    status = pmSourcePhotometryModel (&source->extMag, NULL, source->modelEXT);
 	}
     }
@@ -143,4 +146,9 @@
     if (mode & PM_SOURCE_PHOT_WEIGHT) {
         pmSourcePixelWeight (&source->pixWeight, model, source->maskObj, maskVal);
+    }
+
+    // measure the contribution of included pixels
+    if (mode & PM_SOURCE_PHOT_DIFFSTATS) {
+        pmSourceMeasureDiffStats (source, maskVal);
     }
 
@@ -217,21 +225,26 @@
 
 // return source model magnitude
-bool pmSourcePhotometryModel (float *fitMag, pmModel *model)
-{
-    PS_ASSERT_PTR_NON_NULL(fitMag, false);
-    if (model == NULL) {
-        return false;
-    }
-
-    float fitSum = 0;
-    *fitMag = NAN;
+bool pmSourcePhotometryModel (float *fitMag, float *fitFlux, pmModel *model)
+{
+    psAssert (fitMag || fitFlux, "at least one of magnitude or flux must be requested (not NULL)");
+    if (model == NULL) return false;
+
+    float mag  = NAN;
+    float flux = NAN;
 
     // measure fitMag
-    fitSum = model->modelFlux (model->params);
-    if (fitSum <= 0)
-        return false;
-    if (!isfinite(fitSum))
-        return false;
-    *fitMag = -2.5*log10(fitSum);
+    flux = model->modelFlux (model->params);
+    if (flux > 0) {
+	mag = -2.5*log10(flux);
+    }
+    if (fitMag) {
+	*fitMag = mag;
+    }
+    if (fitFlux) {
+	*fitFlux = flux;
+    }
+
+    if (flux <= 0) return false;
+    if (!isfinite(flux)) return false;
 
     return (true);
@@ -356,4 +369,55 @@
 
     *pixWeight = validSum / modelSum;
+    return (true);
+}
+
+# define FLUX_LIMIT 3.0
+
+// return source aperture magnitude
+bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal)
+{
+    PS_ASSERT_PTR_NON_NULL(source, false);
+
+    if (source->diffStats == NULL) {
+	source->diffStats = pmSourceDiffStatsAlloc();
+    }
+
+    float fGood = 0.0;
+    float fBad  = 0.0;
+    int   nGood = 0;
+    int   nMask = 0;
+    int   nBad  = 0;
+    
+    psImage *flux     = source->pixels;
+    psImage *variance = source->variance;
+    psImage *mask     = source->maskObj;
+
+    for (int iy = 0; iy < flux->numRows; iy++) {
+	for (int ix = 0; ix < flux->numCols; ix++) {
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) {
+		nMask ++;
+                continue;
+	    }
+
+	    float SN = flux->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
+
+	    if (SN > +FLUX_LIMIT) { 
+		nGood ++;
+		fGood += fabs(flux->data.F32[iy][ix]);
+	    }
+
+	    if (SN < -FLUX_LIMIT) { 
+		nBad ++;
+		fBad += fabs(flux->data.F32[iy][ix]);
+	    }
+	}
+    }
+
+    source->diffStats->nGood      = nGood;
+    source->diffStats->fRatio     = (fGood + fBad         == 0.0) ? NAN : fGood / (fGood + fBad);	   
+    source->diffStats->nRatioBad  = (nGood + nBad         == 0)   ? NAN : nGood / (float) (nGood + nBad);	   
+    source->diffStats->nRatioMask = (nGood + nMask        == 0)   ? NAN : nGood / (float) (nGood + nMask);	   
+    source->diffStats->nRatioAll  = (nGood + nMask + nBad == 0)   ? NAN : nGood / (float) (nGood + nMask + nBad);
+
     return (true);
 }
Index: trunk/psModules/src/objects/pmSourcePhotometry.h
===================================================================
--- trunk/psModules/src/objects/pmSourcePhotometry.h	(revision 27530)
+++ trunk/psModules/src/objects/pmSourcePhotometry.h	(revision 27531)
@@ -29,13 +29,15 @@
 
 typedef enum {
-    PM_SOURCE_PHOT_NONE   = 0x0000,
-    PM_SOURCE_PHOT_GROWTH = 0x0001,
-    PM_SOURCE_PHOT_APCORR = 0x0002,
-    PM_SOURCE_PHOT_WEIGHT = 0x0004,
-    PM_SOURCE_PHOT_INTERP = 0x0008,
+    PM_SOURCE_PHOT_NONE      = 0x0000,
+    PM_SOURCE_PHOT_GROWTH    = 0x0001,
+    PM_SOURCE_PHOT_APCORR    = 0x0002,
+    PM_SOURCE_PHOT_WEIGHT    = 0x0004,
+    PM_SOURCE_PHOT_INTERP    = 0x0008,
+    PM_SOURCE_PHOT_DIFFSTATS = 0x0010,
 } pmSourcePhotometryMode;
 
 bool pmSourcePhotometryModel(
     float *fitMag,                      ///< integrated fit magnitude
+    float *fitFlux,                     ///< integrated fit magnitude
     pmModel *model                      ///< model used for photometry
 );
@@ -54,4 +56,5 @@
 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal, const float covarFactor);
 
+bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal);
 
 double pmSourceDataDotModel (const pmSource *Mi, const pmSource *Mj, const bool unweighted_sum, const float covarFactor);
Index: trunk/psModules/src/psmodules.h
===================================================================
--- trunk/psModules/src/psmodules.h	(revision 27530)
+++ trunk/psModules/src/psmodules.h	(revision 27531)
@@ -115,4 +115,6 @@
 #include <pmDetections.h>
 #include <pmMoments.h>
+#include <pmSourceExtendedPars.h>
+#include <pmSourceDiffStats.h>
 #include <pmResiduals.h>
 #include <pmGrowthCurve.h>