Index: /trunk/psphot/src/apply_psf_model.c
===================================================================
--- /trunk/psphot/src/apply_psf_model.c	(revision 4115)
+++ /trunk/psphot/src/apply_psf_model.c	(revision 4116)
@@ -1,6 +1,8 @@
 # include "psphot.h"
 
-// fit psf model to all objects not already fitted (modelPSF == NULL)
+// fit psf model to all objects 
+// PSFSTAR objects will be refitted
 // run this function to a specific flux limit?
+
 
 bool apply_psf_model (psImage *image, psMetadata *config, psArray *sources, pmPSF *psf, psStats *sky) 
@@ -35,8 +37,4 @@
 	psSource *source = sources->data[i];
 
-	// skip the existing PSF stars 
-	// XXX drop this -- refit?
-	if (source->modelPSF != NULL) continue;
-
 	// skip non-astronomical objects (very likely defects)
 	if (source->type == PS_SOURCE_DEFECT) continue;
@@ -54,4 +52,5 @@
 	if (isnan(model->radius)) {
 	  fprintf (stderr, "error in radius\n");
+	  continue;
 	}
 	    
Index: /trunk/psphot/src/mark_psf_source.c
===================================================================
--- /trunk/psphot/src/mark_psf_source.c	(revision 4115)
+++ /trunk/psphot/src/mark_psf_source.c	(revision 4116)
@@ -22,5 +22,5 @@
     if (source->modelPSF == NULL) return (false);
 
-    if (source->type == PS_SOURCE_PSFSTAR) return (true);
+    // SATSTAR objects are fitted with PSF regardless
     if (source->type == PS_SOURCE_SATSTAR) return (true);
 
@@ -36,8 +36,13 @@
     // assign PS_SOURCE_BRIGHTSTAR to bright objects within PSF region of dparams[]
     if ((fabs(nSx) < shapeNsigma) && (fabs(nSy) < shapeNsigma)) {
+      if (source->type == PS_SOURCE_PSFSTAR) return (true);
       source->type = PS_SOURCE_BRIGHTSTAR;
       return (true);
     }
     
+    if (source->type == PS_SOURCE_PSFSTAR) {
+      psLogMsg ("psphot", 3, "PSFSTAR demoted based on dSx, dSy\n");
+    }
+
     if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) {
       source->type = PS_SOURCE_GALAXY;
Index: /trunk/psphot/src/psPolynomials.c
===================================================================
--- /trunk/psphot/src/psPolynomials.c	(revision 4115)
+++ /trunk/psphot/src/psPolynomials.c	(revision 4116)
@@ -453,30 +453,37 @@
 			    const psVector* dz)
 {
+    PS_VECTOR_CHECK_NULL(mask, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_NULL(z, NULL);
+    PS_VECTOR_CHECK_NULL(dz, NULL);
+
     psVector *zFit   = NULL;
-    psVector *zResid = NULL;
-    psStats  *stats  = NULL;
+    psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64);
+    psStats  *stats  = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
 
     for (int N = 0; N < 3; N++) {
 	poly   = VectorFitPolynomial2DOrd_EAM (poly, mask, x, y, z, dz);
 	zFit   = Polynomial2DEvalVectorD (poly, x, y);
-	zResid = (psVector *) psBinaryOp (NULL, (void *) z, "-", (void *) zFit);
-
-	stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
+	zResid = (psVector *) psBinaryOp (zResid, (void *) z, "-", (void *) zFit);
+
 	stats  = psVectorStats (stats, zResid, NULL, mask, 1);
-	psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g (%d pts)\n", 
-		 stats->clippedMean, stats->clippedStdev, stats->clippedNvalues);
+	psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n", 
+		 stats->sampleMean, stats->sampleStdev);
 
 	// set mask if pts are not valid
-	int n = 0;
+	// we are masking out any point which is out of range 
+	// recovery is not allowed with this scheme
 	for (int i = 0; i < zResid->n; i++) {
-	    if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) {
-	      mask->data.U8[i] = 1;
-	      continue;
-	    } else {
-	      mask->data.U8[i] = 0;
-	    }	      
-	    n++;
-	}
-    }
+	  if (mask->data.U8[i]) continue;
+	  if (fabs(zResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) {
+	    mask->data.U8[i] = 1;
+	    continue;
+	  }       
+	}
+	psFree (zFit);
+    }
+    psFree (zResid);
+    psFree (stats);
     return (poly);
 }
Index: /trunk/psphot/src/pspsf.c
===================================================================
--- /trunk/psphot/src/pspsf.c	(revision 4115)
+++ /trunk/psphot/src/pspsf.c	(revision 4116)
@@ -172,5 +172,5 @@
     // we are doing a robust fit.  after each pass, we drop points which are 
     // more deviant than three sigma. 
-    // mask is currently updated for each pass. this is inconsistent
+    // mask is currently updated for each pass. this is inconsistent?
 
     for (int i = 0; i < psf->params->n; i++) {
Index: /trunk/psphot/src/subtract_psf_source.c
===================================================================
--- /trunk/psphot/src/subtract_psf_source.c	(revision 4115)
+++ /trunk/psphot/src/subtract_psf_source.c	(revision 4116)
@@ -9,4 +9,5 @@
   psSource *source = sources->data[i];
 
+  // non-stellar sources are ignored
   if (source->type == PS_SOURCE_GALAXY) return (false);
   if (source->type == PS_SOURCE_DEFECT) return (false);