Index: trunk/psphot/src/psPolynomials.c
===================================================================
--- trunk/psphot/src/psPolynomials.c	(revision 4977)
+++ trunk/psphot/src/psPolynomials.c	(revision 5068)
@@ -497,4 +497,52 @@
 
 // XXX EAM : be careful here with F32 vs F64 vectors 
+psPolynomial1D* VectorClipFitPolynomial1D_EAM(psPolynomial1D* poly,
+					      psStats *stats,
+					      psVector* mask,
+					      psMaskType maskValue,
+					      const psVector* z,
+					      const psVector* dz,
+					      const psVector* x)
+{
+    psVector *zFit   = NULL;
+    psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64);
+
+    // XXX EAM : use SAMPLE_MEAN and SAMPLE_STDEV for stats:
+    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    for (int N = 0; N < 3; N++) {
+	int Nkeep = 0;
+
+	poly   = VectorFitPolynomial1D_EAM (poly, mask, maskValue, z, dz, x);
+	zFit   = Polynomial1DEvalVector_EAM (poly, x);
+	zResid = (psVector *) psBinaryOp (zResid, (void *) z, "-", (void *) zFit);
+
+	stats  = psVectorStats (stats, zResid, NULL, mask, maskValue);
+	psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n", 
+		 stats->sampleMedian, stats->sampleStdev);
+
+	// set mask if pts are not valid
+	// 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 (mask->data.U8[i]) continue;
+	    if (fabs(zResid->data.F64[i] - stats->sampleMedian) > 2*stats->sampleStdev) {
+		mask->data.U8[i] |= 0x01;
+		continue;
+	    }       
+	    Nkeep ++;
+	}
+
+	psTrace (".psphot.RobustFit", 4, "keeping %d of %d pts for fit\n", 
+		 Nkeep, x->n);
+
+	psFree (zFit);
+    }
+    psFree (zResid);
+    return (poly);
+}
+
+
+// XXX EAM : be careful here with F32 vs F64 vectors 
 psPolynomial2D* VectorClipFitPolynomial2D_EAM(psPolynomial2D* poly,
 					      psStats *stats,
@@ -510,7 +558,9 @@
 
     // XXX EAM : use SAMPLE_MEAN and SAMPLE_STDEV for stats:
-    stats->options |= (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
 
     for (int N = 0; N < 3; N++) {
+	int Nkeep = 0;
+
 	poly   = VectorFitPolynomial2D_EAM (poly, mask, maskValue, z, dz, x, y);
 	zFit   = Polynomial2DEvalVectorD_EAM (poly, x, y);
@@ -519,5 +569,5 @@
 	stats  = psVectorStats (stats, zResid, NULL, mask, maskValue);
 	psTrace (".psphot.RobustFit", 4, "residual stats for robust fit:  %g +/- %g\n", 
-		 stats->sampleMean, stats->sampleStdev);
+		 stats->sampleMedian, stats->sampleStdev);
 
 	// set mask if pts are not valid
@@ -526,9 +576,20 @@
 	for (int i = 0; i < zResid->n; i++) {
 	    if (mask->data.U8[i]) continue;
-	    if (fabs(zResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) {
-		mask->data.U8[i] &= 0x01;
+	    if (fabs(zResid->data.F64[i] - stats->sampleMedian) > 3*stats->sampleStdev) {
+		mask->data.U8[i] |= 0x01;
 		continue;
 	    }       
-	}
+	    Nkeep ++;
+	}
+
+	psTrace (".psphot.RobustFit", 4, "keeping %d of %d pts for fit\n", 
+		 Nkeep, x->n);
+
+//	psTrace (".psphot.RobustFit", 4, "model pars: %f %f %f %f\n", 
+//		 poly->coeff[0][0], poly->coeff[1][0], 
+//		 poly->coeff[0][1], poly->coeff[1][1]);
+
+	psTrace (".psphot.RobustFit", 4, "model pars: %f\n", poly->coeff[0][0]);
+
 	psFree (zFit);
     }