Index: /branches/eam_rel8_b1/psModules/src/objects/pmPSFtry.c
===================================================================
--- /branches/eam_rel8_b1/psModules/src/objects/pmPSFtry.c	(revision 5599)
+++ /branches/eam_rel8_b1/psModules/src/objects/pmPSFtry.c	(revision 5600)
@@ -183,6 +183,6 @@
 
     // XXX this function wants aperture radius for pmSourcePhotometry
-    pmPSFtryMetric (try
-                    , RADIUS);
+    pmPSFtryMetric_Alt (try
+                        , RADIUS);
     psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
               modelName, try
@@ -277,10 +277,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);
@@ -304,12 +298,8 @@
     }
 
-    // linear fit to rfBin, daBin
-    psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1);
-
-    // XXX EAM : this is the intended API (cycle 7? cycle 8?)
-    poly = psVectorFitPolynomial1D(poly, maskB, 1, daBin, NULL, rfBin);
-
-    // XXX EAM : replace this when the above version is implemented
-    // poly = psVectorFitPolynomial1DOrd(poly, maskB, rfBin, daBin, NULL);
+    // linear clipped fit to rfBin, daBin
+    psPolynomial1D *poly = psPolynomial1DAlloc(1, PS_POLYNOMIAL_ORD);
+    psStats *fitstat = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    poly = psVectorClipFitPolynomial1D (poly, fitstat, maskB, 1, daBin, NULL, rfBin);
 
     psVector *daBinFit = psPolynomial1DEvalVector (poly, rfBin);
@@ -334,5 +324,76 @@
     psFree (poly);
     psFree (stats);
+    psFree (fitstat);
 
     return true;
 }
+
+bool pmPSFtryMetric_Alt (pmPSFtry *try
+                         , float RADIUS)
+{
+
+    // the measured (aperture - fit) magnitudes (dA == try->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 (try
+                                     ->sources->n, PS_TYPE_F64);
+    for (int i = 0; i < try
+                ->sources->n; i++) {
+            if (try
+                    ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
+            rflux->data.F64[i] = pow(10.0, 0.4*try
+                                     ->fitMag->data.F64[i]);
+        }
+
+    // XXX EAM : try 3hi/1lo sigma clipping on the rflux vs metric fit
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    // XXX EAM
+    stats->min = 1.0;
+    stats->max = 3.0;
+    stats->clipIter = 3;
+
+    // linear clipped fit to rfBin, daBin
+    psPolynomial1D *poly = psPolynomial1DAlloc (1, PS_POLYNOMIAL_ORD);
+    poly = psVectorClipFitPolynomial1D (poly, stats, try
+                                            ->mask, PSFTRY_MASK_ALL, try
+                                                ->metric, NULL, rflux);
+    fprintf (stderr, "fit stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
+
+    try
+        ->psf->ApResid = poly->coeff[0];
+    try
+        ->psf->dApResid = stats->sampleStdev;
+    try
+        ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+
+    FILE *f;
+    f = fopen ("apresid.dat", "w");
+    if (f == NULL)
+        psAbort ("pmPSFtry", "can't open output file");
+
+    for (int i = 0; i < try
+                ->sources->n; i++) {
+            fprintf (f, "%3d %8.4f %12.5e %8.4f %3d\n", i, try
+                         ->fitMag->data.F64[i], rflux->data.F64[i], try
+                             ->metric->data.F64[i], try
+                                 ->mask->data.U8[i]);
+        }
+    fclose (f);
+
+    psFree (rflux);
+    psFree (poly);
+    psFree (stats);
+
+    // psFree (daFit);
+    // psFree (daResid);
+
+    return true;
+}