Index: /trunk/psphot/Makefile
===================================================================
--- /trunk/psphot/Makefile	(revision 4374)
+++ /trunk/psphot/Makefile	(revision 4375)
@@ -41,9 +41,21 @@
 $(SRC)/find_defects.$(ARCH).o
 
+FITSOURCE = \
+$(SRC)/fitsource.$(ARCH).o \
+$(SRC)/onesource.$(ARCH).o \
+$(SRC)/fs_args.$(ARCH).o \
+$(SRC)/psphot-utils.$(ARCH).o \
+$(SRC)/psPolynomials.$(ARCH).o \
+$(SRC)/psUtils.$(ARCH).o \
+$(SRC)/setup.$(ARCH).o \
+$(SRC)/LocalSky.$(ARCH).o
+
 psphot: $(BIN)/psphot.$(ARCH)
+$(BIN)/psphot.$(ARCH) : $(PSPHOT)
+$(PSPHOT): $(SRC)/psphot.h
 
-$(BIN)/psphot.$(ARCH) : $(PSPHOT)
-
-$(PSPHOT): $(SRC)/psphot.h
+fitsource: $(BIN)/fitsource.$(ARCH)
+$(BIN)/fitsource.$(ARCH) : $(FITSOURCE)
+$(FITSOURCE): $(SRC)/psphot.h
 
 INSTALL = psphot
Index: /trunk/psphot/src/LocalSky.c
===================================================================
--- /trunk/psphot/src/LocalSky.c	(revision 4374)
+++ /trunk/psphot/src/LocalSky.c	(revision 4375)
@@ -13,7 +13,8 @@
     srcRegion = psRegionForImage (srcRegion, imdata->image, srcRegion);
 
-//  psImage *subImage = psImageSubset(imdata->image, srcRegion);
-//  psImage *subImageMask = psImageSubset(imdata->mask, srcRegion);
-//  psImage *subImageNoise = psImageSubset(imdata->noise, srcRegion);
+    //  use these when psImageSubset is updated
+    //  psImage *subImage = psImageSubset(imdata->image, srcRegion);
+    //  psImage *subImageMask = psImageSubset(imdata->mask, srcRegion);
+    //  psImage *subImageNoise = psImageSubset(imdata->noise, srcRegion);
 
     psImage *subImage = psImageSubset(imdata->image, srcRegion->x0, srcRegion->y0, srcRegion->x1, srcRegion->y1);
@@ -47,7 +48,9 @@
     psImageMaskRegion (mask, srcRegion, AND, 0x7f);
 
-    source->moments = pmMomentsAlloc();
     psF64 tmpF64;
     p_psGetStatValue(myStats, &tmpF64);
+    if (isnan(tmpF64)) return (false);
+
+    source->moments = pmMomentsAlloc();
     source->moments->Sky = (psF32) tmpF64;
 
Index: /trunk/psphot/src/apply_psf_model.c
===================================================================
--- /trunk/psphot/src/apply_psf_model.c	(revision 4374)
+++ /trunk/psphot/src/apply_psf_model.c	(revision 4375)
@@ -11,4 +11,5 @@
     float y;
     int   Nfit = 0;
+    int   Nsub = 0;
     int   Niter = 0;
 
@@ -71,7 +72,10 @@
 
 	mark_psf_source (source, shapeNsigma);
-	subtract_psf_source (source);
+	if (subtract_psf_source (source)) {
+	  Nsub ++;
+	}
     }
     psLogMsg ("psphot", 3, "fit PSF models: %f sec for %d objects (%d total iterations)\n", psTimerMark ("psphot"), Nfit, Niter);
+    psLogMsg ("psphot", 3, "subtracted %d PSF objects\n", Nsub);
     return (true);
 }
Index: /trunk/psphot/src/by_SN.c
===================================================================
--- /trunk/psphot/src/by_SN.c	(revision 4374)
+++ /trunk/psphot/src/by_SN.c	(revision 4375)
@@ -9,4 +9,6 @@
     psF32 fA = (A->moments == NULL) ? 0 : A->moments->SN;
     psF32 fB = (B->moments == NULL) ? 0 : B->moments->SN;
+    if (isnan (fA)) fA = 0;
+    if (isnan (fB)) fB = 0;
 
     psF32 diff = fA - fB;
Index: /trunk/psphot/src/choose_psf_model.c
===================================================================
--- /trunk/psphot/src/choose_psf_model.c	(revision 4374)
+++ /trunk/psphot/src/choose_psf_model.c	(revision 4375)
@@ -3,5 +3,5 @@
 // test PSF models and select best option
 
-pmPSF *choose_psf_model (psMetadata *config, psArray *sources) 
+pmPSF *choose_psf_model (psMetadata *config, psArray *sources, psStats *skystats) 
 { 
     bool        status;
@@ -34,4 +34,25 @@
 	modelName = psMetadataLookupPtr (&status, config, key);
 	tests->data[i] = pmPSF_TestModel (stars, modelName, RADIUS);
+
+	// dump psf test data to file:
+	{
+	  sprintf (key, "metric.%d.dat", i);
+	  FILE *f = fopen (key, "w");
+	  test = tests->data[i];
+
+
+	  for (int j = 0; j < test->sources->n; j++) {
+	    psSource *source = test->sources->data[j];
+	    psModel  *model  = test->modelPSF->data[j];	    
+	    float sky = 0;
+
+	    if (model != NULL) sky = model->params->data.F32[0];
+
+	    fprintf (f, "%3d %6.1f %6.1f  %7.4f %7.4f  %8.4f\n", 
+		     j, source->moments->x, source->moments->y,
+		     test->metric->data.F64[j], test->fitMag->data.F64[j], sky);
+	  }
+	  fclose (f);
+	}
     }
 
@@ -40,8 +61,8 @@
     test = tests->data[0];
     int   bestN = 0;
-    float bestM = test->metricStats->clippedStdev;
+    float bestM = test->dApResid;
     for (int i = 1; i < Ntest; i++) {
 	test = tests->data[i];
-	int M = test->metricStats->clippedStdev;
+	int M = test->dApResid;
 	if (M < bestM) {
 	    bestM = M;
@@ -51,5 +72,5 @@
     test = tests->data[bestN];
     modelName = psModelGetType (test->modelType);
-    psLogMsg ("psphot.pspsf", 3, "selected psf model %s, metric: %f +/- %f\n", modelName, test->metricStats->clippedMean, test->metricStats->clippedStdev);
+    psLogMsg ("psphot.pspsf", 3, "selected psf model %s, ApResid: %f +/- %f\n", modelName, test->ApResid, test->dApResid);
 
     // set source->model based on best psf model fit
Index: /trunk/psphot/src/fitsource.c
===================================================================
--- /trunk/psphot/src/fitsource.c	(revision 4375)
+++ /trunk/psphot/src/fitsource.c	(revision 4375)
@@ -0,0 +1,21 @@
+# include "psphot.h"
+
+int main (int argc, char **argv) {
+
+    psMetadata  *config  = NULL;		// user-provided configuration information
+    psImageData *imdata  = NULL;
+    float x, y;
+
+    config = fs_args (&argc, argv);
+
+    // load input data (image and config)
+    // create or load mask and noise images
+    imdata = setup (config);
+
+    x = atof (argv[4]);
+    y = atof (argv[5]);
+
+    onesource (imdata, config, argv[3], x, y);
+
+    exit (0);
+}
Index: /trunk/psphot/src/fs_args.c
===================================================================
--- /trunk/psphot/src/fs_args.c	(revision 4375)
+++ /trunk/psphot/src/fs_args.c	(revision 4375)
@@ -0,0 +1,49 @@
+# include "psphot.h"
+
+psMetadata *fs_args (int *argc, char **argv) {
+
+  int N, Nfail;
+
+  // basic pslib options
+  psLogSetFormat ("M");
+  psLogArguments (argc, argv);
+  psTraceArguments (argc, argv);
+
+  // identify options in args list - these go on config 
+
+  char *mask = NULL;
+  if ((N = get_argument (*argc, argv, "--mask"))) {
+    remove_argument (N, argc, argv);
+    mask = psStringCopy (argv[N]);
+    remove_argument (N, argc, argv);
+  }
+
+  char *noise = NULL;
+  if ((N = get_argument (*argc, argv, "--noise"))) {
+    remove_argument (N, argc, argv);
+    noise = psStringCopy (argv[N]);
+    remove_argument (N, argc, argv);
+  }
+
+  if (*argc != 6) fs_usage ();
+
+  psMetadata *config = psMetadataParseConfig (NULL, &Nfail, argv[2], FALSE);
+
+  psMetadataAddStr (config, PS_LIST_HEAD, "INPUT",  "", argv[1]);
+  if (mask != NULL) {
+    psMetadataAddStr (config, PS_LIST_HEAD, "MASK", "", mask);
+  }
+  if (noise != NULL) {
+    psMetadataAddStr (config, PS_LIST_HEAD, "NOISE", "", noise);
+  }
+  return (config);
+}
+
+int fs_usage () {
+
+    fprintf (stderr, "USAGE: psphot (image.fits) (config) (model) (x) (y)\n");
+    fprintf (stderr, "options: \n");
+    fprintf (stderr, "  --mask  (filename)\n");
+    fprintf (stderr, "  --noise (filename)\n");
+    exit (2);
+}
Index: /trunk/psphot/src/load_args.c
===================================================================
--- /trunk/psphot/src/load_args.c	(revision 4374)
+++ /trunk/psphot/src/load_args.c	(revision 4375)
@@ -45,6 +45,6 @@
     fprintf (stderr, "USAGE: psphot (image.fits) (output.fits) (config)\n");
     fprintf (stderr, "options: \n");
-    fprintf (stderr, "  --mask  (filename)");
-    fprintf (stderr, "  --noise (filename)");
+    fprintf (stderr, "  --mask  (filename)\n");
+    fprintf (stderr, "  --noise (filename)\n");
     exit (2);
 }
Index: /trunk/psphot/src/mark_psf_source.c
===================================================================
--- /trunk/psphot/src/mark_psf_source.c	(revision 4374)
+++ /trunk/psphot/src/mark_psf_source.c	(revision 4375)
@@ -14,4 +14,5 @@
 // any object which meets the criterion is marked as 
 // PS_SOURCE_BRIGHTSTAR 
+# define MIN_DS 0.01
 bool mark_psf_source (psSource *source, float shapeNsigma)
 { 
@@ -30,6 +31,9 @@
     dSY = source->modelPSF->dparams->data.F32[5];
 
-    nSx = dSX * SX * SN;
-    nSy = dSY * SY * SN;
+    nSx = dSX / hypot (MIN_DS, 1 / (SX * SN));
+    nSy = dSY / hypot (MIN_DS, 1 / (SY * SN));
+    // sigma = 1 / SX
+    // dsx = 1 / (SX * SN)
+    // dsx_o = hypot (1/SX*SN + MIN_DSX)
 
     // assign PS_SOURCE_BRIGHTSTAR to bright objects within PSF region of dparams[]
Index: /trunk/psphot/src/onesource.c
===================================================================
--- /trunk/psphot/src/onesource.c	(revision 4375)
+++ /trunk/psphot/src/onesource.c	(revision 4375)
@@ -0,0 +1,41 @@
+# include "psphot.h"
+
+bool onesource (psImageData *imdata, psMetadata *config, char *modelName, float x, float y) {
+
+    bool status;
+    int modelType;
+    float sky, obsMag, fitMag;
+
+    modelType   = psModelSetType (modelName);
+    if (modelType < 0) psAbort ("fitsource", "unknown model %s", modelName);
+
+    float INNER  = psMetadataLookupF32 (&status, config, "INNER_RADIUS");
+    float OUTER  = psMetadataLookupF32 (&status, config, "OUTER_RADIUS");
+    float MRAD   = psMetadataLookupF32 (&status, config, "PSF_MOMENTS_RADIUS");
+    float RADIUS = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
+
+    psSource *source = pmSourceAlloc();
+    source->peak = pmPeakAlloc (x, y, 0, 0);
+    pmSourceDefinePixels (source, imdata, x, y, OUTER);
+
+    status = pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER);
+    status = pmSourceMoments (source, MRAD);
+    source->peak->counts = source->moments->Peak;
+
+    psModel *model = pmSourceModelGuess (source, modelType); 
+
+    psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
+    status = pmSourceFitModel (source, model, false);
+
+    pmSourcePhotometry (&fitMag, &obsMag, model, source->pixels, source->mask);
+    fprintf (stderr, "ap: %f, fit: %f, apmifit: %f\n", obsMag, fitMag, obsMag - fitMag);
+
+    // subtract object, leave local sky
+    sky = model->params->data.F32[0];
+    model->params->data.F32[0] = 0;
+    pmSourceSubModel (source->pixels, source->mask, model, false);
+    model->params->data.F32[0] = sky;
+    
+    DumpImage (source->pixels, "resid.fits");
+    DumpImage (source->mask, "mask.fits");
+}
Index: /trunk/psphot/src/psPolynomials.c
===================================================================
--- /trunk/psphot/src/psPolynomials.c	(revision 4374)
+++ /trunk/psphot/src/psPolynomials.c	(revision 4375)
@@ -1,65 +1,42 @@
 # include "psphot.h"
 
-// XXX EAM : this version uses myPoly->nX as Norder, not Nterms
-psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly,
-				 const psVector *x,
-				 const psVector *y)
-
-{
-    PS_POLY_CHECK_NULL(myPoly, NULL);
-    PS_VECTOR_CHECK_NULL(x, NULL);
-    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
-    PS_VECTOR_CHECK_NULL(y, NULL);
-    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
-
-    psVector *tmp;
-    psS32 vecLen=x->n;
-
-    // Determine the length of the output vector to by the minimum of the x,y vectors
-    if (y->n < vecLen) {
-        vecLen = y->n;
-    }
-
-    // Create output vector to return
-    tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
-
-    // Evaluate the polynomial at the specified points
-    for (psS32 i=0; i<vecLen; i++) {
-        tmp->data.F32[i] = Polynomial2DEval(myPoly,x->data.F32[i],y->data.F32[i]);
-    }
-
-    // Return output vector
-    return(tmp);
-}
-
-// XXX EAM : this version uses the F64 vectors
-psVector *Polynomial2DEvalVectorD(const psPolynomial2D *myPoly,
-				  const psVector *x,
-				  const psVector *y)
-
+// write out the terms of the given 1D polynomial
+void psPolynomial1DDump (psPolynomial1D *poly) {
+
+    for (int i = 0; i < poly->n + 1; i++) {
+	fprintf (stderr, "x^%d : %g +/- %g\n", i, poly->coeff[i], poly->coeffErr[i]);
+    }
+}    
+
+psF32 Polynomial1DEval_EAM(psF32 x, const psPolynomial1D* myPoly)
+{
+    psS32 loop_x = 0;
+    psF32 polySum = 0.0;
+    psF32 xSum = 1.0;
+
+    for (loop_x = 0; loop_x < myPoly->n + 1; loop_x++) {
+        if (myPoly->mask[loop_x] == 0) {
+            polySum += xSum * myPoly->coeff[loop_x];
+        }
+        xSum *= x;
+    }
+
+    return(polySum);
+}
+
+psVector *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly,
+				       const psVector *x)
 {
     PS_POLY_CHECK_NULL(myPoly, NULL);
     PS_VECTOR_CHECK_NULL(x, NULL);
     PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
-    PS_VECTOR_CHECK_NULL(y, NULL);
-    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
 
     psVector *tmp;
-    psS32 vecLen=x->n;
-
-    // Determine the length of the output vector to by the minimum of the x,y vectors
-    if (y->n < vecLen) {
-        vecLen = y->n;
-    }
-
-    // Create output vector to return
-    tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
-
-    // Evaluate the polynomial at the specified points
-    for (psS32 i=0; i<vecLen; i++) {
-        tmp->data.F64[i] = Polynomial2DEval(myPoly,x->data.F64[i],y->data.F64[i]);
-    }
-
-    // Return output vector
+
+    tmp = psVectorAlloc(x->n, PS_TYPE_F64);
+    for (psS32 i=0;i<x->n;i++) {
+        tmp->data.F64[i] = Polynomial1DEval_EAM(x->data.F64[i], myPoly);
+    }
+
     return(tmp);
 }
@@ -93,44 +70,4 @@
 }
 
-// XXX EAM : use Nterm = Norder + 1 definition  
-// the user requests a polynomial of order Norder
-psPolynomial2D* Polynomial2DAlloc(psS32 nXorder, psS32 nYorder,
-                                    psPolynomialType type)
-{
-    PS_INT_CHECK_NON_NEGATIVE(nXorder, NULL);
-    PS_INT_CHECK_NON_NEGATIVE(nYorder, NULL);
-
-    psS32 x = 0;
-    psS32 y = 0;
-    psS32 nXterm = nXorder + 1;
-    psS32 nYterm = nYorder + 1;
-    psPolynomial2D* newPoly = NULL;
-
-    newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D));
-    // p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial2DFree);
-    // XXX EAM : me, being lazy
-
-    newPoly->type = type;
-    newPoly->nX = nXorder;
-    newPoly->nY = nYorder;
-
-    newPoly->coeff = (psF32 **)psAlloc(nXterm * sizeof(psF32 *));
-    newPoly->coeffErr = (psF32 **)psAlloc(nXterm * sizeof(psF32 *));
-    newPoly->mask = (psU8 **)psAlloc(nXterm * sizeof(psU8 *));
-    for (x = 0; x < nXterm; x++) {
-        newPoly->coeff[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32));
-        newPoly->coeffErr[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32));
-        newPoly->mask[x] = (psU8 *)psAlloc(nYterm * sizeof(psU8));
-    }
-    for (x = 0; x < nXterm; x++) {
-        for (y = 0; y < nYterm; y++) {
-            newPoly->coeff[x][y] = 0.0;
-            newPoly->coeffErr[x][y] = 0.0;
-            newPoly->mask[x][y] = 0;
-        }
-    }
-    return(newPoly);
-}
-
 // XXX EAM : my alternate BuildSums1D
 static psVector *BuildSums1D(psVector* sums, 
@@ -153,36 +90,4 @@
         sums->data.F64[i] = xSum;
         xSum *= x;
-    }
-    return (sums);
-}
-
-// XXX EAM : BuildSums2D in analogy with BuildSums1D
-static psImage *BuildSums2D(psImage* sums,
-			    psF64 x,      psF64 y,
-			    psS32 nXterm, psS32 nYterm)
-{
-    psS32 nXsum = 0;
-    psS32 nYsum = 0;
-    psF64 xSum = 1.0;
-    psF64 ySum = 1.0;
-
-    nXsum = 2*nXterm;
-    nYsum = 2*nYterm;
-    if (sums == NULL) {
-        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
-    }
-    if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) {
-	psFree (sums);
-        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
-    }
-
-    ySum = 1.0;
-    for (int j = 0; j < nYsum; j++) {
-	xSum = ySum;
-	for (int i = 0; i < nXsum; i++) {
-	    sums->data.F64[i][j] = xSum;
-	    xSum *= x;
-	}
-	ySum *= y;
     }
     return (sums);
@@ -291,4 +196,142 @@
 }
 
+// ********************** 2D polynomial functions ******************
+
+// XXX EAM : this version uses myPoly->nX as Norder, not Nterms
+psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly,
+				 const psVector *x,
+				 const psVector *y)
+
+{
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
+
+    psVector *tmp;
+    psS32 vecLen=x->n;
+
+    // Determine the length of the output vector to by the minimum of the x,y vectors
+    if (y->n < vecLen) {
+        vecLen = y->n;
+    }
+
+    // Create output vector to return
+    tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
+
+    // Evaluate the polynomial at the specified points
+    for (psS32 i=0; i<vecLen; i++) {
+        tmp->data.F32[i] = Polynomial2DEval(myPoly,x->data.F32[i],y->data.F32[i]);
+    }
+
+    // Return output vector
+    return(tmp);
+}
+
+// XXX EAM : this version uses the F64 vectors
+psVector *Polynomial2DEvalVectorD(const psPolynomial2D *myPoly,
+				  const psVector *x,
+				  const psVector *y)
+
+{
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
+
+    psVector *tmp;
+    psS32 vecLen=x->n;
+
+    // Determine the length of the output vector to by the minimum of the x,y vectors
+    if (y->n < vecLen) {
+        vecLen = y->n;
+    }
+
+    // Create output vector to return
+    tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
+
+    // Evaluate the polynomial at the specified points
+    for (psS32 i=0; i<vecLen; i++) {
+        tmp->data.F64[i] = Polynomial2DEval(myPoly,x->data.F64[i],y->data.F64[i]);
+    }
+
+    // Return output vector
+    return(tmp);
+}
+
+// XXX EAM : use Nterm = Norder + 1 definition  
+// the user requests a polynomial of order Norder
+psPolynomial2D* Polynomial2DAlloc(psS32 nXorder, psS32 nYorder,
+                                    psPolynomialType type)
+{
+    PS_INT_CHECK_NON_NEGATIVE(nXorder, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nYorder, NULL);
+
+    psS32 x = 0;
+    psS32 y = 0;
+    psS32 nXterm = nXorder + 1;
+    psS32 nYterm = nYorder + 1;
+    psPolynomial2D* newPoly = NULL;
+
+    newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D));
+    // p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial2DFree);
+    // XXX EAM : me, being lazy
+
+    newPoly->type = type;
+    newPoly->nX = nXorder;
+    newPoly->nY = nYorder;
+
+    newPoly->coeff = (psF32 **)psAlloc(nXterm * sizeof(psF32 *));
+    newPoly->coeffErr = (psF32 **)psAlloc(nXterm * sizeof(psF32 *));
+    newPoly->mask = (psU8 **)psAlloc(nXterm * sizeof(psU8 *));
+    for (x = 0; x < nXterm; x++) {
+        newPoly->coeff[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32));
+        newPoly->coeffErr[x] = (psF32 *)psAlloc(nYterm * sizeof(psF32));
+        newPoly->mask[x] = (psU8 *)psAlloc(nYterm * sizeof(psU8));
+    }
+    for (x = 0; x < nXterm; x++) {
+        for (y = 0; y < nYterm; y++) {
+            newPoly->coeff[x][y] = 0.0;
+            newPoly->coeffErr[x][y] = 0.0;
+            newPoly->mask[x][y] = 0;
+        }
+    }
+    return(newPoly);
+}
+
+// XXX EAM : BuildSums2D in analogy with BuildSums1D
+static psImage *BuildSums2D(psImage* sums,
+			    psF64 x,      psF64 y,
+			    psS32 nXterm, psS32 nYterm)
+{
+    psS32 nXsum = 0;
+    psS32 nYsum = 0;
+    psF64 xSum = 1.0;
+    psF64 ySum = 1.0;
+
+    nXsum = 2*nXterm;
+    nYsum = 2*nYterm;
+    if (sums == NULL) {
+        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
+    }
+    if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) {
+	psFree (sums);
+        sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64);
+    }
+
+    ySum = 1.0;
+    for (int j = 0; j < nYsum; j++) {
+	xSum = ySum;
+	for (int i = 0; i < nXsum; i++) {
+	    sums->data.F64[i][j] = xSum;
+	    xSum *= x;
+	}
+	ySum *= y;
+    }
+    return (sums);
+}
+
 // XXX EAM : test version of 2d fitting
 psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly,
@@ -387,12 +430,4 @@
 }    
 
-// write out the terms of the given 1D polynomial
-void psPolynomial1DDump (psPolynomial1D *poly) {
-
-    for (int i = 0; i < poly->n + 1; i++) {
-	fprintf (stderr, "x^%d : %g +/- %g\n", i, poly->coeff[i], poly->coeffErr[i]);
-    }
-}    
-
 psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly,
 			    const psVector* x,
Index: /trunk/psphot/src/psUtils.c
===================================================================
--- /trunk/psphot/src/psUtils.c	(revision 4374)
+++ /trunk/psphot/src/psUtils.c	(revision 4375)
@@ -448,2 +448,14 @@
 }
 
+psVector *psVectorCreate (double lower, double upper, double delta, psElemType type) {
+
+  int nBin = (upper - lower) / delta;
+
+  psVector *out = psVectorAlloc (nBin, type);
+  
+  for (int i = 0; i < nBin; i++) {
+    out->data.F64[i] = lower + i * delta;
+  }
+
+  return (out);
+}
Index: /trunk/psphot/src/psphot-utils.c
===================================================================
--- /trunk/psphot/src/psphot-utils.c	(revision 4374)
+++ /trunk/psphot/src/psphot-utils.c	(revision 4375)
@@ -77,14 +77,86 @@
 	params = model->params;
 	dparams = model->dparams;
-	fprintf (f, "%7.1f %7.1f  %5.1f %7.1f  %7.4f  ", 
+	fprintf (f, "%7.1f %7.1f  %5.1f %7.3f  %7.4f  ", 
 		 params[0].data.F32[2], params[0].data.F32[3], 
-		 params[0].data.F32[0], params[0].data.F32[1], (dparams[0].data.F32[1]/params[0].data.F32[1]));
-	for (j = 0; j < model->params->n - 4; j++) {
-	    fprintf (f, "%9.3g ", params[0].data.F32[j+4]);
-	}
-	for (j = 0; j < model->params->n - 4; j++) {
-	    fprintf (f, "%9.3g ", dparams[0].data.F32[j+4]);
+		 params[0].data.F32[0], -2.5*log10(params[0].data.F32[1]), (dparams[0].data.F32[1]/params[0].data.F32[1]));
+	for (j = 0; j < model->params->n - 4; j++) {
+	    fprintf (f, "%9.6f ", params[0].data.F32[j+4]);
+	}
+	for (j = 0; j < model->params->n - 4; j++) {
+	    fprintf (f, "%9.6f ", dparams[0].data.F32[j+4]);
 	}
 	fprintf (f, ": %2d %7.4f %7.4f %7.2f (%d %d)\n", source[0].type, log10(model[0].chisq), source[0].moments->SN, model[0].radius, model[0].nDOF, model[0].nIter);
+    }
+    fclose (f);
+    return true;
+}
+
+// dump the sources to an output file
+bool DumpModelFLT (psArray *sources, char *filename) 
+{
+
+    int i, j;
+    FILE *f;
+    psVector *params;
+    psVector *dparams;
+    psModel  *model;
+
+    f = fopen (filename, "w");
+    if (f == NULL) {
+	psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename);
+	return false;
+    }
+
+    // write sources with models first 
+    for (i = 0; i < sources->n; i++) {
+	psSource *source = (psSource *) sources->data[i];
+	model = (psModel  *) source->modelFLT;
+	if (model == NULL) continue;
+	params = model->params;
+	dparams = model->dparams;
+	fprintf (f, "%7.1f %7.1f  %5.1f %7.3f  %7.4f  ", 
+		 params[0].data.F32[2], params[0].data.F32[3], 
+		 params[0].data.F32[0], -2.5*log10(params[0].data.F32[1]), (dparams[0].data.F32[1]/params[0].data.F32[1]));
+	for (j = 0; j < model->params->n - 4; j++) {
+	    fprintf (f, "%9.6f ", params[0].data.F32[j+4]);
+	}
+	for (j = 0; j < model->params->n - 4; j++) {
+	    fprintf (f, "%9.6f ", dparams[0].data.F32[j+4]);
+	}
+	fprintf (f, ": %2d %7.4f %7.4f %7.2f (%d %d)\n", source[0].type, log10(model[0].chisq), source[0].moments->SN, model[0].radius, model[0].nDOF, model[0].nIter);
+    }
+    fclose (f);
+    return true;
+}
+
+// dump the sources to an output file
+bool DumpModelNULL (psArray *sources, char *filename) 
+{
+
+    int i, j;
+    FILE *f;
+    psVector *params;
+    psVector *dparams;
+    psModel  *model;
+
+    f = fopen (filename, "w");
+    if (f == NULL) {
+	psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename);
+	return false;
+    }
+
+    // write sources with models first 
+    for (i = 0; i < sources->n; i++) {
+	psSource *source = (psSource *) sources->data[i];
+	if (source->modelFLT != NULL) continue;
+	if (source->modelPSF != NULL) continue;
+	if (source->moments == NULL) continue;
+	fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %8.1f %7.1f %7.1f %7.1f  %4d %2d\n", 
+		 source->peak->x, source->peak->y, source->peak->counts, 
+		 source->moments->x, source->moments->y, 
+		 source->moments->Sx, source->moments->Sy, 
+		 source->moments->Sum, source->moments->Peak, 
+		 source->moments->Sky, source->moments->SN, 
+		 source->moments->nPixels, source->type); 
     }
     fclose (f);
@@ -118,8 +190,8 @@
 		 params[0].data.F32[0], params[0].data.F32[1]);
 	for (j = 0; j < model->params->n - 4; j++) {
-	    fprintf (f, "%9.3g ", params[0].data.F32[j+4]);
-	}
-	for (j = 0; j < model->params->n - 4; j++) {
-	    fprintf (f, "%9.3g ", dparams[0].data.F32[j+4]);
+	    fprintf (f, "%9.6f ", params[0].data.F32[j+4]);
+	}
+	for (j = 0; j < model->params->n - 4; j++) {
+	    fprintf (f, "%9.6f ", dparams[0].data.F32[j+4]);
 	}
 	fprintf (f, ":  %7.1f (%d %d)\n", model[0].chisq, model[0].nDOF, model[0].nIter);
@@ -138,4 +210,28 @@
 }
 
+bool pmSourcePhotometry (float *fitMag, float *obsMag, psModel *model, psImage *image, psImage *mask) {
+
+    float obsSum = 0;
+    float fitSum = 0;
+    float sky = model->params->data.F32[0];
+
+    psModelFlux modelFluxFunc = psModelFlux_GetFunction (model->type);
+    fitSum = modelFluxFunc (model->params);
+
+    for (int ix = 0; ix < image->numCols; ix++) {
+	for (int iy = 0; iy < image->numRows; iy++) {
+	    if (mask->data.U8[iy][ix]) continue;
+	    obsSum += image->data.F32[iy][ix] - sky;
+	    // fitSum += psModelEval (model, image, ix, iy) - sky;
+	}
+    }
+    if (obsSum <= 0) return false;
+    if (fitSum <= 0) return false;
+
+    *fitMag = -2.5*log10(fitSum);
+    *obsMag = -2.5*log10(obsSum);
+    return (true);
+} 
+
 # if (0)
 void pmSourceMaskSaturated (psSource *source, float saturate) {
Index: /trunk/psphot/src/psphot.c
===================================================================
--- /trunk/psphot/src/psphot.c	(revision 4374)
+++ /trunk/psphot/src/psphot.c	(revision 4375)
@@ -29,8 +29,8 @@
 
     // use stellar objects SN > PSF_SN_LIM
-    psf = choose_psf_model (config, sources);
+    psf = choose_psf_model (config, sources, sky);
 
     // identify defects based on 1-pixel wide x and y direction 2nd moments
-    find_defects (zap, sources, config, psf);
+    // find_defects (zap, sources, config, psf);
 
     // test PSF on all except SATURATE and DEFECT (artifacts)
@@ -46,5 +46,9 @@
  
     DumpImage (imdata->image, argv[2]);
-    DumpModelPSF (sources, "sources.dat");
+    DumpModelPSF (sources, "psfsources.dat");
+    DumpModelFLT (sources, "fltsources.dat");
+    DumpModelNULL (sources, "nullsources.dat");
+    DumpMoments (sources, "moments.dat");
+
     exit (0);
 }
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 4374)
+++ /trunk/psphot/src/psphot.h	(revision 4375)
@@ -5,4 +5,7 @@
 # include <unistd.h>
 # include <stdlib.h>
+# include <math.h>
+
+# define  M_PI 3.14159265358979323846264338328      /* pi */
 
 typedef struct {
@@ -22,14 +25,41 @@
   psVector   *mask;
   psVector   *metric;
-  psStats    *metricStats;
+  psVector   *fitMag;
+  float       ApResid;
+  float       dApResid;
+  float       skyBias;
 } pmPSF_Test;
 
-bool psTimerStart (char *name);
-psF64 psTimerMark (char *name);
+// used by mask operations
+# define AND true
+# define OR false
 
-bool psImageInit (psImage *image,...);
-void psImageSmooth (psImage *image, float sigma, float Nsigma);
-psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in);
 
+// top-level psphot functions
+psImageData *setup (psMetadata *config);
+psStats *image_stats (psImageData *imdata, psMetadata *config);
+psArray *find_peaks (psImageData *imdata, psMetadata *config, psStats *sky);
+psArray *source_moments (psImageData *imdata, psMetadata *config, psArray *allpeaks);
+pmPSF *choose_psf_model (psMetadata *config, psArray *sources, psStats *sky);
+bool apply_psf_model (psImageData *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
+bool test_psf_scatter (psArray *sources);
+bool mark_psf_sources (psArray *sources, psMetadata *config);
+bool subtract_psf_sources (psArray *sources);
+bool fit_galaxies (psImage *image, psMetadata *config, psArray *sources);
+bool subtract_galaxies (psArray *sources, psMetadata *config);
+bool subtract_psf_source (psSource *source);
+bool mark_psf_source (psSource *source, float shapeNsigma);
+bool find_defects (psImage *zapmask, psArray *sources, psMetadata *config, pmPSF *psf);
+
+// psf utilities
+pmPSF *pmPSFAlloc (psModelType type);
+pmPSF_Test *pmPSF_TestAlloc (psArray *stars, char *modelName);
+pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float radius);
+bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask);
+psModel *psModelFromPSF (psModel *model, pmPSF *psf);
+bool pmSourcePhotometry (float *fitMag, float *obsMag, psModel *model, psImage *image, psImage *mask);
+bool pmPSFMetricModel (pmPSF_Test *test, float RADIUS);
+
+// output options
 bool DumpPeaks (psArray *sources, char *filename);
 bool DumpMoments (psArray *sources, char *filename);
@@ -38,61 +68,17 @@
 bool DumpImage (psImage *image, char *filename);
 
+// psphot utilities
 int usage ();
-int get_argument (int argc, char **argv, char *arg);
-int remove_argument (int N, int *argc, char **argv);
 psMetadata *load_args (int *argc, char **argv);
-
+int by_SN (const void **a, const void **b);
 psArray *pmPeaksSubset(psArray *peaks, psF32 maxValue, const psRegion *valid);
 psArray *SelectPSFStars (psArray *sources);
-
-psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name);
-
-psF32 Polynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y);
-psImage *psBuildSums2D(psImage* sums,psF64 x,psF64 y,psS32 nXterm, psS32 nYterm);
-psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* yErr);
-psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
-psPolynomial2D* Polynomial2DAlloc(psS32 nXorder, psS32 nYorder, psPolynomialType type);
-void psPolynomial2DDump (psPolynomial2D *poly);
-psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
-psPolynomial2D* RobustFit2D(psPolynomial2D* poly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
-psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly, const psVector *x,const psVector *y);
-psVector *psBuildSums1D(psVector* sums, psF64 x,psS32 nTerm);
-void psPolynomial1DDump (psPolynomial1D *poly);
-
-psS32 psLogArguments (int *argc, char **argv);
-psS32 psTraceArguments (int *argc, char **argv);
-psVector *p_psGetRowVectorFromImage(psImage *image, psU32 row);
-psVector *psGetRowVectorFromImage(psImage *image, psU32 row);
-
-pmPSF *pmPSFAlloc (psModelType type);
-pmPSF_Test *pmPSF_TestAlloc (psArray *stars, char *modelName);
-pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float radius);
-
-bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask);
-psModel *psModelFromPSF (psModel *model, pmPSF *psf);
-void pmSourceMaskRegion (psSource *source, psRegion *region);
-void pmSourceMaskSaturated (psSource *source, float saturate);
-
-// top-level psphot functions
-psImageData *setup (psMetadata *config);
-psStats *image_stats (psImageData *imdata, psMetadata *config);
-psArray *find_peaks (psImageData *imdata, psMetadata *config, psStats *sky);
-psArray *source_moments (psImageData *imdata, psMetadata *config, psArray *allpeaks);
-pmPSF *choose_psf_model (psMetadata *config, psArray *sources);
-bool apply_psf_model (psImageData *imdata, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky);
-bool test_psf_scatter (psArray *sources);
-bool mark_psf_sources (psArray *sources, psMetadata *config);
-bool subtract_psf_sources (psArray *sources);
-bool fit_galaxies (psImage *image, psMetadata *config, psArray *sources);
-bool subtract_galaxies (psArray *sources, psMetadata *config);
-
-int by_SN (const void **a, const void **b);
-bool subtract_psf_source (psSource *source);
-bool mark_psf_source (psSource *source, float shapeNsigma);
-bool find_defects (psImage *zapmask, psArray *sources, psMetadata *config, pmPSF *psf);
-
 bool pmSourceDefinePixels(psSource *mySource, const psImageData *imdata, psF32 x, psF32 y, psF32 Radius);
 bool pmSourceLocalSky_EAM (psSource *source, psStatsOptions statsOptions, psF32 Radius);
-psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius);
+
+// fitsource utilities
+bool onesource (psImageData *imdata, psMetadata *config, char *modelName, float x, float y);
+psMetadata *fs_args (int *argc, char **argv);
+int fs_usage ();
 
 // image mask functions
@@ -102,4 +88,42 @@
 void psImageKeepCircle (psImage *image, double x, double y, double radius, bool logical_and, int maskValue);
 
-# define AND true
-# define OR false
+// basic image functions
+bool psImageInit (psImage *image,...);
+void psImageSmooth (psImage *image, float sigma, float Nsigma);
+psRegion *psRegionForImage (psRegion *out, psImage *image, psRegion *in);
+psRegion *psRegionSquare (psF32 x, psF32 y, psF32 radius);
+
+// misc global functions
+psVector *p_psGetRowVectorFromImage(psImage *image, psU32 row);
+psVector *psGetRowVectorFromImage(psImage *image, psU32 row);
+
+// utilities
+bool psTimerStart (char *name);
+psF64 psTimerMark (char *name);
+psS32 psLogArguments (int *argc, char **argv);
+psS32 psTraceArguments (int *argc, char **argv);
+int get_argument (int argc, char **argv, char *arg);
+int remove_argument (int N, int *argc, char **argv);
+psF32 pmConfigLookupF32 (bool *status, psMetadata *config, psMetadata *header, char *name);
+psVector *psVectorCreate (double lower, double upper, double delta, psElemType type);
+
+// polynomial functions
+psF32 Polynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y);
+psImage *psBuildSums2D(psImage* sums,psF64 x,psF64 y,psS32 nXterm, psS32 nYterm);
+psPolynomial2D* VectorFitPolynomial2DOrd_EAM(psPolynomial2D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
+psPolynomial2D* Polynomial2DAlloc(psS32 nXorder, psS32 nYorder, psPolynomialType type);
+void psPolynomial2DDump (psPolynomial2D *poly);
+psPolynomial2D* RobustFit2D_nomask(psPolynomial2D* poly, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
+psPolynomial2D* RobustFit2D(psPolynomial2D* poly, psVector* mask, const psVector* x, const psVector* y, const psVector* z, const psVector* zErr);
+psVector *Polynomial2DEvalVector(const psPolynomial2D *myPoly, const psVector *x,const psVector *y);
+
+psVector *psBuildSums1D(psVector* sums, psF64 x,psS32 nTerm);
+void psPolynomial1DDump (psPolynomial1D *poly);
+psF32 Polynomial1DEval_EAM(psF32 x, const psPolynomial1D* myPoly);
+psVector *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly, const psVector *x);
+psPolynomial1D* Polynomial1DAlloc(psS32 nOrder, psPolynomialType type);
+psPolynomial1D* VectorFitPolynomial1DOrd_EAM(psPolynomial1D* myPoly, psVector* mask, const psVector* x, const psVector* y, const psVector* yErr);
+
+// deprecated
+void pmSourceMaskRegion (psSource *source, psRegion *region);
+void pmSourceMaskSaturated (psSource *source, float saturate);
Index: /trunk/psphot/src/pspsf.c
===================================================================
--- /trunk/psphot/src/pspsf.c	(revision 4374)
+++ /trunk/psphot/src/pspsf.c	(revision 4375)
@@ -34,25 +34,36 @@
     test->modelFLT    = psArrayAlloc (sources->n);
     test->modelPSF    = psArrayAlloc (sources->n);
-    test->metricStats = NULL;
     test->metric      = psVectorAlloc (sources->n, PS_TYPE_F64);
+    test->fitMag      = psVectorAlloc (sources->n, PS_TYPE_F64);
     test->mask        = psVectorAlloc (sources->n, PS_TYPE_U8);
+    test->ApResid     = 0;
+    test->dApResid    = 0;
+    test->skyBias     = 0;
 
     for (int i = 0; i < test->modelFLT->n; i++) {
+	test->mask->data.U8[i]  = 0;
 	test->modelFLT->data[i] = NULL;
 	test->modelPSF->data[i] = NULL;
-	test->mask->data.U8[i]  = 0;
+	test->metric->data.F64[i] = 0;
+	test->fitMag->data.F64[i] = 0;
     }	
     return (test);
 }
 
+// test->mask values indicate reason source was rejected:
+// 1: outlier in psf polynomial fit
+// 2: flt model failed to converge
+// 3: psf model failed to converge
+// 4: invalid source photometry
 pmPSF_Test *pmPSF_TestModel (psArray *sources, char *modelName, float RADIUS) 
 {
-  
+    bool status;
+    float obsMag;
+    float fitMag;
     float x;
     float y;
-    bool status;
     int Nflt = 0;
     int Npsf = 0;
-    psVector *metricMask = psVectorAlloc (sources->n, PS_TYPE_U8);
+
     pmPSF_Test *test = pmPSF_TestAlloc (sources, modelName);
 
@@ -66,11 +77,12 @@
 
 	// set temporary object mask and fit object
-	// fit model as FLT, not PSF (mask & skip poor fits)
+	// fit model as FLT, not PSF
 	psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
 	status = pmSourceFitModel (source, model, false);
 	psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
 
+	// exclude the poor fits
 	if (!status) {
-	  test->mask->data.U8[i] = 1;
+	  test->mask->data.U8[i] = 2;
 	  continue;
 	}
@@ -84,26 +96,16 @@
     // stage 2: construct a psf (pmPSF) from this collection of model fits
     pmPSFFromModels (test->psf, test->modelFLT, test->mask);
-    
-    // count valid sources
-    int Nkeep = 0;
-    for (int i = 0; i < sources->n; i++) {
-      if (test->mask->data.U8[i]) continue;
-      Nkeep++;
-    }
-    psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates\n", Nkeep, sources->n);
 
     // stage 3: refit with fixed shape parameters
     psTimerStart ("fit");
     for (int i = 0; i < test->sources->n; i++) {
+	// masked for: bad model fit, outlier in parameters
+	if (test->mask->data.U8[i]) continue;
+
 	psSource *source = test->sources->data[i];
 	psModel  *modelFLT = test->modelFLT->data[i];
 
-	// masked for: bad model fit, outlier in parameters
-	if (test->mask->data.U8[i]) continue;
-	if (modelFLT == NULL) {
-	  psLogMsg ("psphot.psftest", 2, "programming error: missing model not masked\n");
-	  continue;
-	}
-	psModel  *modelPSF = psModelFromPSF (modelFLT, test->psf); // set shape for this model
+	// set shape for this model based on PSF
+	psModel *modelPSF = psModelFromPSF (modelFLT, test->psf); 
 	x = source->peak->x;
 	y = source->peak->y;
@@ -111,13 +113,29 @@
 	psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
 	status = pmSourceFitModel (source, modelPSF, true);
-	psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
 
 	// skip poor fits
 	if (!status) {
-	  test->mask->data.U8[i] = 1;
-	  continue;
+	    test->mask->data.U8[i] = 3;
+	    goto next_source;
 	}
+
+	// otherwise, save the resulting model
 	test->modelPSF->data[i] = modelPSF;
+
+	// XXX : use a different aperture radius from the fit radius?
+	// XXX : use a different estimator for the local sky?
+	// XXX : pass 'source' as input?
+	if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
+	    test->mask->data.U8[i] = 4;
+	    goto next_source;
+	}	    
+
+	test->metric->data.F64[i] = obsMag - fitMag;
+	test->fitMag->data.F64[i] = fitMag;
 	Npsf ++;
+
+    next_source:
+	psImageKeepCircle (source->mask, x, y, RADIUS, AND, 0x7f);
+
     }
     psLogMsg ("psphot.psftest", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
@@ -125,51 +143,8 @@
     DumpModelFits (test->modelPSF, "modelsPSF.dat");
 
-    // stage 4: measure PSF metric (Ap-Fit)
-    // should I be saving the aperture photometry for all sources?
-    // save the metric for all stars
-    for (int i = 0; i < test->sources->n; i++) {
-
-	float obsSum = 0;
-	float fitSum = 0;
-      
-	// is this metricMask redundant with test->mask ?
-	metricMask->data.U8[i] = 1;
-	test->metric->data.F64[i] = 0;
-
-	if (test->mask->data.U8[i]) continue;
-
-	psModel *model = test->modelPSF->data[i];
-	if (model == NULL) {
-	  psLogMsg ("psphot.psftest", 2, "programming error: missing model not masked\n");
-	  continue;
-	}
-
-        psImage *image = ((psSource *)test->sources->data[i])->pixels;
-        psImage *mask  = ((psSource *)test->sources->data[i])->mask;
-
-	// this metric is Ap-Fit
-	float sky = model->params->data.F32[0];
-	for (int ix = 0; ix < image->numCols; ix++) {
-	    for (int iy = 0; iy < image->numRows; iy++) {
-		if (mask->data.U8[iy][ix]) continue;
-		obsSum += image->data.F32[iy][ix] - sky;
-		fitSum += psModelEval (model, image, ix, iy) - sky;
-	    }
-	}
-	// fprintf (stderr, "%d %f %f  %f\n", i, obsSum, fitSum, test->metric->data.F64[i]);
-	// keep the good metrics
-	if ((fitSum > 0) && (obsSum > 0)) {
-	    metricMask->data.U8[i] = 0;
-	    test->metric->data.F64[i] = -2.5*log10(obsSum/fitSum);
-	} else {
-	  test->mask->data.U8[i] = 1;
-	}
-    }
-
-    // XXX use robust stats, not clipped stats
-    psStats *stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
-    psVectorStats (stats, test->metric, NULL, metricMask, 1);
-    test->metricStats = stats;
-    psLogMsg ("psphot.pspsf", 3, "test model %s, metric: %f +/- %f\n", modelName, stats->clippedMean, stats->clippedStdev);
+    // XXX this function wants aperture radius from pmSourcePhotometry
+    pmPSFMetricModel (test, RADIUS);
+    psLogMsg ("psphot.pspsf", 3, "test model %s, ap-fit: %f +/- %f, sky bias: %f\n", 
+	      modelName, test->ApResid, test->dApResid, test->skyBias);
     return (test);
 }
@@ -181,7 +156,7 @@
     // construct the fit vectors from the collection of objects
     // use the mask to ignore missing fits 
-    psVector *x = psVectorAlloc (models->n, PS_TYPE_F64);
-    psVector *y = psVectorAlloc (models->n, PS_TYPE_F64);
-    psVector *z = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *x  = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *y  = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *z  = psVectorAlloc (models->n, PS_TYPE_F64);
     psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64);
 
@@ -208,17 +183,8 @@
 	    // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4
 	}
-	// if (n != x->n) psAbort ("pmPSFFromModels", "programming error: mismatched number of NULL models\n");
 
 	psf->params->data[i] = RobustFit2D (psf->params->data[i], mask, x, y, z, dz);
+	// psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
 	// psPolynomial2DDump (psf->params->data[i]);
-
-	// count valid sources
-	int Nkeep = 0;
-	for (int j = 0; j < mask->n; j++) {
-	  if (mask->data.U8[j]) continue;
-	  Nkeep++;
-	}
-	psTrace ("psphot.psftest", 3, "keeping %d of %d PSF candidates (PSF param %d)\n", Nkeep, mask->n, i);
-
     }
     return (true);
@@ -241,2 +207,104 @@
 }
 
+bool pmPSFMetricModel (pmPSF_Test *test, float RADIUS) {
+
+  float dBin;
+  int   nKeep, nSkip;
+  
+
+  // the measured (aperture - fit) magnitudes (dA == test->metric)
+  //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
+  //     dA = dAo + dsky/flux
+  //   where flux is the flux of the star
+  // we fit this trend to find the infinite flux aperture correction (dAo),
+  //   the nominal sky bias (dsky), and the error on dAo
+  // the values of dA are contaminated by stars with close neighbors in the aperture
+  //   we use an outlier rejection to avoid this bias
+
+  // rflux = ten(0.4*fitMag);
+  psVector *rflux = psVectorAlloc (test->sources->n, PS_TYPE_F64);
+  for (int i = 0; i < test->sources->n; i++) {
+    if (test->mask->data.U8[i]) continue;
+    rflux->data.F64[i] = pow(10.0, 0.4*test->fitMag->data.F64[i]);
+  }
+  // psScalar *t1 = psScalar(0.4);
+  // psVector *v1 = psBinaryOp (NULL, test->fitMag, "*", t1);
+  // psVector *rflux = psUnaryOp (NULL, v1, "ten");
+  // psFree (t1);
+  // psFree (v1);
+
+  // find min and max of (1/flux):
+  psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
+  psVectorStats (stats, rflux, NULL, test->mask, 0xff);
+  
+  // build binned versions of rflux, metric
+  dBin = (stats->max - stats->min) / 10.0;
+  psVector *rfBin = psVectorCreate (stats->min, stats->max, dBin, PS_TYPE_F64);
+  psVector *daBin = psVectorAlloc (rfBin->n, PS_TYPE_F64);
+  psVector *maskB = psVectorAlloc (rfBin->n, PS_TYPE_U8);
+  psFree (stats);
+
+  psTrace ("psphot.metricmodel", 3, "rflux max: %g, min: %g, delta: %g\n", stats->max, stats->min, dBin);
+
+  for (int i = 0; i < daBin->n; i++) {
+
+    psVector *tmp = psVectorAlloc (test->sources->n, PS_TYPE_F64);
+    tmp->n = 0;
+
+    // accumulate data within bin range
+    for (int j = 0; j < test->sources->n; j++) {
+      // masked for: bad model fit, outlier in parameters
+      if (test->mask->data.U8[j]) continue;
+    
+      // skip points with extreme dA values
+      if (fabs(test->metric->data.F64[j]) > 0.5) continue;
+
+      // skip points outside of this bin
+      if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin) continue;
+      if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin) continue;
+
+      tmp->data.F64[tmp->n] = test->metric->data.F64[j];
+      tmp->n ++;
+    }
+
+    // is this a valid point?
+    maskB->data.U8[i] = 0;
+    if (tmp->n < 2) {
+      maskB->data.U8[i] = 1;
+      continue;
+    } 
+
+    // dA values are contaminated with low outliers 
+    // measure statistics only on upper 50% of points
+    psVectorSort (tmp, tmp);
+    nKeep = 0.5*tmp->n;
+    nSkip = tmp->n - nKeep;
+
+    psVector *tmp2 = psVectorAlloc (nKeep, PS_TYPE_F64);
+    for (int j = 0; j < tmp2->n; j++) {
+      tmp2->data.F64[j] = tmp->data.F64[j + nSkip];
+    }
+
+    stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+    psVectorStats (stats, tmp2, NULL, NULL, 0);
+    psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
+
+    daBin->data.F64[i] = stats->sampleMedian;
+  }
+
+  // linear fit to rfBin, daBin
+  psPolynomial1D *poly = Polynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
+  poly = VectorFitPolynomial1DOrd_EAM (poly, maskB, rfBin, daBin, NULL);
+
+  psVector *daBinFit = Polynomial1DEvalVector_EAM (poly, rfBin);
+  psVector *daResid  = (psVector *) psBinaryOp (NULL, (void *) daBin, "-", (void *) daBinFit);
+
+  stats = psStatsAlloc (PS_STAT_CLIPPED_STDEV);
+  stats = psVectorStats (stats, daResid, NULL, maskB, 1);
+
+  test->ApResid = poly->coeff[0];
+  test->dApResid = stats->clippedStdev;
+  test->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+
+  return true;
+}
Index: /trunk/psphot/src/setup.c
===================================================================
--- /trunk/psphot/src/setup.c	(revision 4374)
+++ /trunk/psphot/src/setup.c	(revision 4375)
@@ -18,5 +18,5 @@
     psTimerStart ("psphot");
 
-    // setup header template, specify COMMENT and HISTORY as multis
+    // setup header template, specify COMMENT and HISTORY as MULTI
     header = psMetadataAlloc ();
     psMetadataAdd (header, PS_LIST_HEAD, "COMMENT", PS_META_MULTI, "folder for comment", NULL);
@@ -89,4 +89,5 @@
 	}
     }
+    DumpImage (mask, "mask.fits");
 
     psLogMsg ("psphot", 3, "load data: %f sec\n", psTimerMark ("psphot"));
Index: /trunk/psphot/src/source_moments.c
===================================================================
--- /trunk/psphot/src/source_moments.c	(revision 4374)
+++ /trunk/psphot/src/source_moments.c	(revision 4375)
@@ -5,4 +5,5 @@
     bool     status  = false;
     psArray *sources = NULL;
+    pmPSFClump psfClump;
 
     psTimerStart ("psphot");
@@ -25,14 +26,19 @@
 
 	// this should use ROBUST not SAMPLE median, but it is broken
-	pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER);
-	pmSourceMoments (source, RADIUS);
+	status = pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER);
+	if (!status) continue;
+
+	status = pmSourceMoments (source, RADIUS);
+	if (!status) continue;
+	
 	psArrayAdd (sources, 100, source);
     }
 
     // group into STAR, COSMIC, GALAXY, SATURATED
-    pmSourceRoughClass (sources, config);
+    psfClump = pmSourcePSFClump (sources, config);
+    pmSourceRoughClass (sources, config, psfClump);
 
     // make this optional
-    DumpMoments (sources, "moments.dat");
+    // DumpMoments (sources, "moments.dat");
 
     psLogMsg ("psphot", 3, "moments: %f sec\n", psTimerMark ("psphot"));
@@ -41,2 +47,4 @@
     return (sources);
 }
+
+// filter out bad peaks (eg, no valid pixels available for sky)