Index: /trunk/psphot/src/pmPSF.c
===================================================================
--- /trunk/psphot/src/pmPSF.c	(revision 5067)
+++ /trunk/psphot/src/pmPSF.c	(revision 5068)
@@ -40,5 +40,5 @@
 	// XXX EAM : make this a user-defined value?
 	// XXX EAM : future version (0.7.0?) psf->params->data[i] = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 1, 1);
-	psf->params->data[i] = Polynomial2DAlloc_EAM(PS_POLYNOMIAL_ORD, 1, 1);
+	psf->params->data[i] = Polynomial2DAlloc_EAM(PS_POLYNOMIAL_ORD, 0, 0);
     }
 
Index: /trunk/psphot/src/pmPSFtry.c
===================================================================
--- /trunk/psphot/src/pmPSFtry.c	(revision 5067)
+++ /trunk/psphot/src/pmPSFtry.c	(revision 5068)
@@ -183,4 +183,8 @@
   //   we use an outlier rejection to avoid this bias
 
+  FILE *f;
+  f = fopen ("apresid.dat", "w");
+  if (f == NULL) psAbort ("pmPSFtry", "can't open output file");
+
   // rflux = ten(0.4*fitMag);
   psVector *rflux = psVectorAlloc (try->sources->n, PS_TYPE_F64);
@@ -188,5 +192,7 @@
     if (try->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
     rflux->data.F64[i] = pow(10.0, 0.4*try->fitMag->data.F64[i]);
+    fprintf (f, "%3d %8.4f %12.5e %8.4f\n", i, try->fitMag->data.F64[i], rflux->data.F64[i], try->metric->data.F64[i]);
   }
+  fclose (f);
 
   // find min and max of (1/flux):
@@ -236,10 +242,4 @@
     // measure statistics only on upper 50% of points
     // this would be easier if we could sort in reverse:
-    //
-    // psVectorSort (tmp, tmp);
-    // tmp->n = 0.5*tmp->n;
-    // stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-    // psVectorStats (stats, tmp, NULL, NULL, 0);
-    // psTrace ("psphot.metricmodel", 4, "rfBin %d (%g): %d pts, %g\n", i, rfBin->data.F64[i], tmp->n, stats->sampleMedian);
 
     psVectorSort (tmp, tmp);
@@ -267,5 +267,8 @@
 
   // XXX EAM : this is the intended API (cycle 7? cycle 8?) 
-  poly = VectorFitPolynomial1D_EAM (poly, maskB, 1, daBin, NULL, rfBin);
+  psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+  poly = VectorClipFitPolynomial1D_EAM (poly, fitstat, maskB, 1, daBin, NULL, rfBin);
+
+  // poly = VectorFitPolynomial1D_EAM (poly, maskB, 1, daBin, NULL, rfBin);
 
   // XXX EAM : replace this when the above version is implemented
@@ -290,4 +293,5 @@
   psFree (poly);
   psFree (stats);
+  psFree (fitstat);
 
   return true;
Index: unk/psphot/src/pmSourceUtils.c
===================================================================
--- /trunk/psphot/src/pmSourceUtils.c	(revision 5067)
+++ 	(revision )
@@ -1,1 +1,0 @@
-# include "psphot.h"
Index: /trunk/psphot/src/psLibUtils.h
===================================================================
--- /trunk/psphot/src/psLibUtils.h	(revision 5067)
+++ /trunk/psphot/src/psLibUtils.h	(revision 5068)
@@ -62,4 +62,5 @@
 psVector       *Polynomial1DEvalVector_EAM(const psPolynomial1D *myPoly, const psVector *x);
 psPolynomial1D *VectorFitPolynomial1D_EAM(psPolynomial1D* myPoly, psVector* mask, psMaskType maskValue, const psVector* y, const psVector* yErr, const psVector* x);
+psPolynomial1D *VectorClipFitPolynomial1D_EAM(psPolynomial1D* poly, psStats *stats, psVector* mask, psMaskType maskValue, const psVector* z, const psVector* zErr, const psVector* x);
 
 psPolynomial2D *Polynomial2DAlloc_EAM(psPolynomialType type, psS32 nXorder, psS32 nYorder);
Index: /trunk/psphot/src/psPolynomials.c
===================================================================
--- /trunk/psphot/src/psPolynomials.c	(revision 5067)
+++ /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);
     }
Index: /trunk/psphot/src/psphot.h
===================================================================
--- /trunk/psphot/src/psphot.h	(revision 5067)
+++ /trunk/psphot/src/psphot.h	(revision 5068)
@@ -42,5 +42,5 @@
 bool 	     pmSourcesWriteCMP (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
 bool 	     pmSourcesWriteCMF (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
-bool 	     pmSourcesWriteSX (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
+bool 	     pmSourcesWriteSX (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *sky);
 int  	     pmSourcesDophotType (pmSource *source);
 bool 	     pmPeaksWriteText (psArray *sources, char *filename);
Index: /trunk/psphot/src/psphotOutput.c
===================================================================
--- /trunk/psphot/src/psphotOutput.c	(revision 5067)
+++ /trunk/psphot/src/psphotOutput.c	(revision 5068)
@@ -31,5 +31,5 @@
   
     if (!strcasecmp (outputMode, "SX")) {
-	pmSourcesWriteSX (imdata, outputFile, sources, psf, sky);
+	pmSourcesWriteSX (imdata, config, outputFile, sources, psf, sky);
 	return;
     }
@@ -159,6 +159,6 @@
 }
 
-// elixir/sextractor-style output list with fixed line width
-bool pmSourcesWriteSX (eamReadout *imdata, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) {
+// elixir-mode / sextractor-style output list with fixed line width
+bool pmSourcesWriteSX (eamReadout *imdata, psMetadata *config, char *filename, psArray *sources, pmPSF *psf, psStats *skyStats) {
 
     int i, status, flags;
@@ -166,6 +166,8 @@
     psF32 *PAR, *dPAR;
     float sky, dmag, apMag, fitMag;
-
-    psLine *line = psLineAlloc (110);  // 104 is dophot-defined line length
+    float x, y, rflux;
+    bool result;
+
+    psLine *line = psLineAlloc (110);  // 110 is sextractor line length
 
     FILE *f = fopen (filename, "w");
@@ -174,4 +176,7 @@
 	return false;
     }
+
+    // aperture radius for ap magnitude
+    float RADIUS = psMetadataLookupF32 (&result, config, "PSF_FIT_RADIUS");
 
     // write sources with models first 
@@ -214,4 +219,12 @@
 	dPAR = model->dparams->data.F32;
 
+	x = source->peak->x;
+	y = source->peak->y;
+
+	// we have already (psphotApplyPSF) defined pixels at least OUTER_RADIUS from source
+	// we need to mask pixels to measure the aperture magnitude
+	// set aperture mask circle of PSF_FIT_RADIUS 
+	psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
+
 	// save local sky for later
 	sky = model->params->data.F32[0];
@@ -224,5 +237,9 @@
 
 	// measure object photometry
-	status = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask);
+	status  = pmSourcePhotometry (&fitMag, &apMag, model, source->pixels, source->mask);
+
+	// correct both apMag and fitMag to same system, consistent with infinite flux star in aperture RADIUS
+	rflux   = pow (10.0, 0.4*fitMag);
+	apMag  -= rflux * psf->skyBias * (M_PI * PS_SQR(RADIUS));
 	fitMag += psf->ApResid;
 
@@ -232,4 +249,7 @@
 	pmSourceSubModel (source->pixels, source->mask, model, false);
 	model->params->data.F32[0] = sky;
+
+	// unmask aperture
+	psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
 
 	if (status == FALSE) continue;