Index: trunk/psphot/src/pspsf.c
===================================================================
--- trunk/psphot/src/pspsf.c	(revision 4114)
+++ trunk/psphot/src/pspsf.c	(revision 4115)
@@ -36,8 +36,10 @@
     test->metricStats = NULL;
     test->metric      = psVectorAlloc (sources->n, PS_TYPE_F64);
+    test->mask        = psVectorAlloc (sources->n, PS_TYPE_U8);
 
     for (int i = 0; i < test->modelFLT->n; i++) {
 	test->modelFLT->data[i] = NULL;
 	test->modelPSF->data[i] = NULL;
+	test->mask->data.U8[i]  = 0;
     }	
     return (test);
@@ -57,6 +59,9 @@
 	psModel  *model  = pmSourceModelGuess (source, test->modelType); 
 
-	// fit model as FLT, not PSF (skip poor fits)
-	if (!pmSourceFitModel (source, model, false)) continue;
+	// fit model as FLT, not PSF (mask & skip poor fits)
+	if (!pmSourceFitModel (source, model, false)) {
+	  test->mask->data.U8[i] = 1;
+	  continue;
+	}
 	test->modelFLT->data[i] = model;
 	Nflt ++;
@@ -66,5 +71,5 @@
 
     // stage 2: construct a psf (pmPSF) from this collection of model fits
-    pmPSFFromModels (test->psf, test->modelFLT);
+    pmPSFFromModels (test->psf, test->modelFLT, test->mask);
 
     // stage 3: refit with fixed shape parameters
@@ -73,9 +78,18 @@
 	psSource *source = test->sources->data[i];
 	psModel  *modelFLT = test->modelFLT->data[i];
-	if (modelFLT == NULL) continue;
+
+	// 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
 
 	// fit model as PSF, not FLT (skip poor fits)
-	if (!pmSourceFitModel (source, modelPSF, true)) continue;
+	if (!pmSourceFitModel (source, modelPSF, true)) {
+	  test->mask->data.U8[i] = 1;
+	  continue;
+	}
 	test->modelPSF->data[i] = modelPSF;
 	Npsf ++;
@@ -92,15 +106,21 @@
 	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) continue;
+	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++) {
@@ -110,10 +130,12 @@
 	    }
 	}
-	// this metric is Ap-Fit
 	// 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;
+	}
     }
 
@@ -128,13 +150,15 @@
 // input: an array of psModels, pre-allocated psf
 // some of the array entries may be NULL, ignore them
-bool pmPSFFromModels (pmPSF *psf, psArray *models) {
+bool pmPSFFromModels (pmPSF *psf, psArray *models, psVector *mask) {
 
     int n;
 
-    // construct the x and y vectors from the collection of objects
-    // count the number of valid model measurements
+    // 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);
-    n = 0;
+    psVector *z = psVectorAlloc (models->n, PS_TYPE_F64);
+    psVector *dz = psVectorAlloc (models->n, PS_TYPE_F64);
+
     for (int i = 0; i < models->n; i++) {
 	psModel *model = models->data[i];
@@ -142,30 +166,24 @@
 
 	// XXX EAM : this is fragile: x and y must be stored consistently at 2,3
-	x->data.F64[n] = model->params->data.F32[2];
-	y->data.F64[n] = model->params->data.F32[3];
-	n++;
-    }
-    x->n = y->n = n;
+	x->data.F64[i] = model->params->data.F32[2];
+	y->data.F64[i] = model->params->data.F32[3];
+    }
 
     // we are doing a robust fit.  after each pass, we drop points which are 
-    // more deviant than three sigma.
-    psVector *z = psVectorAlloc (x->n, PS_TYPE_F64);
-    psVector *dz = psVectorAlloc (x->n, PS_TYPE_F64);
+    // more deviant than three sigma. 
+    // mask is currently updated for each pass. this is inconsistent
+
     for (int i = 0; i < psf->params->n; i++) {
-	n = 0;
 	for (int j = 0; j < models->n; j++) {
 	    psModel *model = models->data[j];
 	    if (model == NULL) continue;
-	    z->data.F64[n] = model->params->data.F32[i + 4];
-	    dz->data.F64[n] = 1;
-	    n++;
+	    z->data.F64[j] = model->params->data.F32[i + 4];
+	    dz->data.F64[j] = 1;
 	    // XXX EAM : need to use actual errors?
 	    // XXX EAM : this is fragile: psf(Nparams) = model(Nparams) - 4
 	}
-	if (n != x->n) psAbort ("pmPSFFromModels", "programming error: mismatched number of NULL models\n");
-	z->n = n;
-	dz->n = n;
-
-	psf->params->data[i] = RobustFit2D (psf->params->data[i], x, y, z, dz);
+	// 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);
 	// psPolynomial2DDump (psf->params->data[i]);
     }