Index: trunk/psphot/src/psphotApResid.c
===================================================================
--- trunk/psphot/src/psphotApResid.c	(revision 14943)
+++ trunk/psphot/src/psphotApResid.c	(revision 14944)
@@ -1,9 +1,5 @@
 # include "psphotInternal.h"
-// XXXX this code fails if there are too few sources to measure the aperture residual
-// the larger problem is that the rules for accepting more polynomial terms are weak.
-
-static pmPSFApTrendOptions DEFAULT_OPTION = PM_PSF_APTREND_SKYBIAS;
-
-// measure the aperture residual statistics
+
+// measure the aperture residual statistics and 2D variations
 bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psMaskType maskVal, psMaskType mark)
 {
@@ -16,5 +12,5 @@
 
     PS_ASSERT_PTR_NON_NULL(psf, false);
-    PS_ASSERT_PTR_NON_NULL(psf->ApTrend, false);
+    // XXX drop this? PS_ASSERT_PTR_NON_NULL(psf->ApTrend, false);
     PS_ASSERT_PTR_NON_NULL(readout, false);
     PS_ASSERT_PTR_NON_NULL(sources, false);
@@ -30,6 +26,7 @@
     bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
     bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP");
-    int NSTAR_APERTURE_CORRECTION_MIN =
-        psMetadataLookupS32(&status, recipe, "NSTAR_APERTURE_CORRECTION_MIN");
+
+    // XXX is this still needed?  the pmTrend2D stuff should be auto-adjusting...
+    int NSTAR_APERTURE_CORRECTION_MIN = psMetadataLookupS32(&status, recipe, "NSTAR_APERTURE_CORRECTION_MIN");
     if (!status) {
         NSTAR_APERTURE_CORRECTION_MIN = 5;
@@ -54,16 +51,10 @@
 
     psVector *mask    = psVectorAllocEmpty (300, PS_TYPE_U8);
-    psVector *xPos    = psVectorAllocEmpty (300, PS_TYPE_F64);
-    psVector *yPos    = psVectorAllocEmpty (300, PS_TYPE_F64);
-    psVector *flux    = psVectorAllocEmpty (300, PS_TYPE_F64);
-    psVector *r2rflux = psVectorAllocEmpty (300, PS_TYPE_F64);
-    psVector *apResid = psVectorAllocEmpty (300, PS_TYPE_F64);
-    psVector *dMag    = psVectorAllocEmpty (300, PS_TYPE_F64);
+    psVector *mag     = psVectorAllocEmpty (300, PS_TYPE_F32);
+    psVector *xPos    = psVectorAllocEmpty (300, PS_TYPE_F32);
+    psVector *yPos    = psVectorAllocEmpty (300, PS_TYPE_F32);
+    psVector *apResid = psVectorAllocEmpty (300, PS_TYPE_F32);
+    psVector *dMag    = psVectorAllocEmpty (300, PS_TYPE_F32);
     Npsf = 0;
-
-    FILE *dumpFile = NULL;
-    if (psTraceGetLevel("psphot") >= 5) {
-	dumpFile = fopen ("apresid.dat", "w");
-    }
 
     // select all good PM_SOURCE_TYPE_STAR entries
@@ -90,41 +81,28 @@
         }
 
-        // sanity clipping : if the model is so discrepant, but your expectation in the recipe
-        apResid->data.F64[Npsf] = source->apMag - source->psfMag;
-
-        if ((MAX_AP_OFFSET > 0) && (fabs(apResid->data.F64[Npsf]) > MAX_AP_OFFSET)) {
+	// aperture residual for this source
+	float dap = source->apMag - source->psfMag;
+
+        // sanity clipping : if the model is very discrepant, put your expectation in the recipe
+        if ((MAX_AP_OFFSET > 0) && (fabs(dap) > MAX_AP_OFFSET)) {
             Nfail ++;
             continue;
         }
 
-        xPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_XPOS];
-        yPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_YPOS];
-
-        flux->data.F64[Npsf]    = pow(10.0, -0.4*source->psfMag);
-        r2rflux->data.F64[Npsf] = PS_SQR(model->radiusFit) / flux->data.F64[Npsf];
+	mag->data.F32[Npsf]     = source->psfMag;
+        apResid->data.F32[Npsf] = dap;
+        xPos->data.F32[Npsf]    = model->params->data.F32[PM_PAR_XPOS];
+        yPos->data.F32[Npsf]    = model->params->data.F32[PM_PAR_YPOS];
 
         mask->data.U8[Npsf] = 0;
 
-        dMag->data.F64[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
-
-	if (psTraceGetLevel("psphot") >= 5) {
-	    fprintf (dumpFile, "%f %f  %f %f %f  %x\n", 
-		     xPos->data.F64[Npsf], yPos->data.F64[Npsf], 
-		     source->psfMag, dMag->data.F64[Npsf], apResid->data.F64[Npsf], 
-		     mask->data.U8[Npsf]);
-	}
+        dMag->data.F32[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
 
         psVectorExtend (mask,    100, 1);
         psVectorExtend (xPos,    100, 1);
         psVectorExtend (yPos,    100, 1);
-        psVectorExtend (flux,    100, 1);
-        psVectorExtend (r2rflux, 100, 1);
         psVectorExtend (dMag,    100, 1);
         psVectorExtend (apResid, 100, 1);
         Npsf ++;
-    }
-
-    if (psTraceGetLevel("psphot") >= 5) {
-	fclose (dumpFile);
     }
     
@@ -139,236 +117,64 @@
     }
 
-    // APTREND options : NONE CONSTANT SKYBIAS XY_LIN XY_QUAD SKY_XY_LIN FULL
-    // APTREND options are used in the switch block below
-    pmPSFApTrendOptions ApTrendOption = DEFAULT_OPTION;
-    char *optionName = psMetadataLookupPtr (&status, recipe, "APTREND");
-    if (status) ApTrendOption = pmPSFApTrendOptionFromName (optionName);
-    if (ApTrendOption == PM_PSF_APTREND_ERROR) {
-        psError(PSPHOT_ERR_APERTURE, true, "invalid aperture residual trend %s", optionName);
-        return false;
-    }
-
+    // XXX deprecating the old code which allowed the ApResid to be fitted as a function of flux and r^2/flux
+    // XXX is this asymmetric clipping still needed?  this analysis should come after neighbors are subtracted...
     // 3hi/1lo sigma clipping on the rflux vs metric fit
-    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    stats->min = 2.0;
-    stats->max = 3.0;
-
-#define P_APTREND_SWITCH_CLEANUP        /* Cleanup memory on error in ApTrendOption switch */ \
-    psFree(psf->growth); psf->growth = NULL; \
-    psFree(mask); \
-    psFree(xPos); \
-    psFree(yPos); \
-    psFree(flux); \
-    psFree(r2rflux); \
-    psFree(apResid); \
-    psFree(dMag); \
-    psFree(stats)
-
-    // no correction
-    switch (ApTrendOption) {
-      case PM_PSF_APTREND_NONE:
-        // remove ApTrend fit from pmPSFtry
-        psf->ApTrend->coeff[0][0][0][0] = 0;
-        break;
-      case PM_PSF_APTREND_CONSTANT:
-        stats->clipIter = 2;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
-        if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        // apply the fit
-        stats->clipIter = 3;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
-        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        break;
-      case PM_PSF_APTREND_SKYBIAS:
-        // first clip out objects which are too far from the median
-        stats->clipIter = 2;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
-        if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        // apply the fit
-        stats->clipIter = 3;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
-        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        break;
-      case PM_PSF_APTREND_SKYSAT:
-        // first clip out objects which are too far from the median
-        stats->clipIter = 2;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
-        if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        // apply the fit
-        stats->clipIter = 2;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
-        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        // apply the fit
-        stats->clipIter = 3;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT);
-        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        break;
-      case PM_PSF_APTREND_XY_LIN:
-        // first clip out objects which are too far from the median
-        stats->clipIter = 2;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
-        if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        // apply the fit
-        stats->clipIter = 3;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN);
-        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        break;
-      case PM_PSF_APTREND_XY_QUAD:
-        // first clip out objects which are too far from the median
-        stats->clipIter = 2;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
-        if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        // apply the fit
-        stats->clipIter = 3;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD);
-        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        break;
-      case PM_PSF_APTREND_SKY_XY_LIN:
-        // first clip out objects which are too far from the median
-        stats->clipIter = 2;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
-        if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        // apply the fit
-        stats->clipIter = 3;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN);
-        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        break;
-      case PM_PSF_APTREND_SKYSAT_XY_LIN:
-        // first clip out objects which are too far from the median
-        stats->clipIter = 2;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
-        if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        // apply the fit
-        stats->clipIter = 3;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
-        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        // apply the fit
-        stats->clipIter = 3;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN);
-        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        break;
-      case PM_PSF_APTREND_ALL:
-        // first clip out objects which are too far from the median
-        stats->clipIter = 2;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
-        if (!psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        // fit just SkyBias and clip out objects which are too far from the median
-        stats->clipIter = 2;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
-        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        // finally, fit x, y, SkyBias and clip out objects which are too far from the median
-        stats->clipIter = 3;
-        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL);
-        if (!psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux)) {
-            psError(PSPHOT_ERR_PHOTOM, false, "fitting all");
-            P_APTREND_SWITCH_CLEANUP;
-            return false;
-        }
-        break;
-      default:
-        psError(PSPHOT_ERR_PHOTOM, true, "Unknown APTREND value: %s", optionName);
-        return false;
-    }
-#undef P_APTREND_SWITCH_CLEANUP
+    // systematic error information
+    float errorScale = 0.0;
+    float errorFloor = 0.0;
+
+    float errorFloorMin = FLT_MAX;
+    int entryMin = -1;
+
+    // *** iterate over spatial scale until error Floor increases
+    // *** fit out the dap vs mag trend?
+    // *** stop if Npsf / (Nx * Ny) < 3
+    for (int i = 1; i < 10; i++) {
+
+	if (!psphotApResidTrend (readout, psf, Npsf, i, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag)) {
+	    break;
+	}
+
+	// store the resulting errorFloor values and the scales, redo the best
+	if (errorFloor < errorFloorMin) {
+	    errorFloorMin = errorFloor;
+	    entryMin = i;
+	}
+    }
+    if (entryMin == -1) {
+	psAbort ("failed on ApResid Trend");
+    }
+
+    psphotApResidTrend (readout, psf, Npsf, entryMin, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag);
 
     // construct the fitted values and the residuals
-    psVector *fit   = psPolynomial4DEvalVector (psf->ApTrend, xPos, yPos, r2rflux, flux);
-    psVector *resid = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) fit);
-
-    // measure scatter for sources with dMag < 0.01 (S/N = 100)
-    int Nkeep = 0;
-    psStats *residStats = psStatsAlloc (PS_STAT_SAMPLE_STDEV);
-    for (int i = 0; i < dMag->n; i++) {
-        if (dMag->data.F64[i] > 0.01) {
-            mask->data.U8[i] |= 0x02;
-        }
-        if (! mask->data.U8[i]) Nkeep ++;
-    }
-    psVectorStats (residStats, resid, NULL, mask, 0x03);
+    psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, xPos, yPos);
+    psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit);
+    psVector *dMagSys = (psVector *) psBinaryOp (NULL, (void *) dMag, "*", (void *) psScalarAlloc(errorScale, PS_TYPE_F32));
+
+    psphotSaveImage (NULL, psf->ApTrend->map->map, "apresid.map.fits");
+
+    if (psTraceGetLevel("psphot") >= 5) {
+	FILE *dumpFile = fopen ("apresid.dat", "w");
+	for (int i = 0; i < xPos->n; i++) {
+	    fprintf (dumpFile, "%f %f  %f %f %f  %f %f  %x\n", 
+		     xPos->data.F32[i], yPos->data.F32[i], 
+		     mag->data.F32[i], dMag->data.F32[i], dMagSys->data.F32[i],
+		     apResid->data.F32[i], apResidRes->data.F32[i],
+		     mask->data.U8[i]);
+	}
+	fclose (dumpFile);
+    }
 
     // apply ApTrend results
-    psf->skyBias  = psf->ApTrend->coeff[0][0][1][0];
-    psf->skySat   = psf->ApTrend->coeff[0][0][0][1];
-    psf->ApResid  = psf->ApTrend->coeff[0][0][0][0];
-    psf->dApResid = residStats->sampleStdev;
-    psf->ApTrend->coeff[0][0][1][0] = 0;
-    psf->ApTrend->coeff[0][0][0][1] = 0;
+    float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
+    float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
+    psf->ApResid  = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
+    psf->dApResid = errorFloor;
     psf->nApResid = Npsf;
 
     // save results for later output
-    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skyBias);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   psf->skySat);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   0.0);
+    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   0.0);
     psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
     psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
@@ -376,27 +182,14 @@
     psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "ap residual scatter", psf->nApResid);
 
-    psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d of %d objects: %f sec\n",
-              Nkeep, Npsf, psTimerMark ("psphot"));
-    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f : %f bias, %f skysat\n",
-              psf->ApResid, psf->dApResid, psf->skyBias, psf->skySat);
-    psLogMsg ("psphot.apresid", PS_LOG_MINUTIA, "apresid trends: %f %f %f %f %f\n",
-              1e3*psf->ApTrend->coeff[1][0][0][0],
-              1e6*psf->ApTrend->coeff[2][0][0][0],
-              1e6*psf->ApTrend->coeff[1][1][0][0],
-              1e3*psf->ApTrend->coeff[0][1][0][0],
-              1e6*psf->ApTrend->coeff[0][2][0][0]);
-
+    // psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d of %d objects: %f sec\n", Nkeep, Npsf, psTimerMark ("psphot"));
+    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);
+
+    psFree (mag);
     psFree (mask);
     psFree (xPos);
     psFree (yPos);
-    psFree (flux);
-    psFree (r2rflux);
     psFree (apResid);
     psFree (dMag);
 
-    psFree (fit);
-    psFree (resid);
-    psFree (stats);
-    psFree (residStats);
     return true;
 }
@@ -408,2 +201,111 @@
 */
 
+bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, int nGroup) {
+
+    psStats *statsS = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    psStats *statsM = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
+
+    // measure the trend in bins with 10 values each; if < 10 total, use them all
+    int nBin = PS_MAX (dMag->n / nGroup, 1);
+
+    // output vectors for ApResid trend 
+    psVector *dSo = psVectorAlloc (nBin, PS_TYPE_F32);
+    psVector *dMo = psVectorAlloc (nBin, PS_TYPE_F32);
+    psVector *dRo = psVectorAlloc (nBin, PS_TYPE_F32);
+
+    // use dMag to group the dMag and dap vectors
+    psVector *index = psVectorSortIndex (NULL, dMag);
+
+    // subset vectors for dMag and dap values within the given range
+    psVector *dMSubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
+    psVector *dASubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
+
+    int n = 0;
+    for (int i = 0; i < dMo->n; i++) {
+	int j;
+	for (j = 0; (j < nGroup) && (n < dMag->n); j++, n++) {
+	    int N = index->data.U32[n];
+	    dMSubset->data.F32[j] = dMag->data.F32[N];
+	    dASubset->data.F32[j] = dap->data.F32[N];
+	}
+	dMSubset->n = j;
+	dASubset->n = j;
+
+	psStatsInit (statsS);
+	psStatsInit (statsM);
+
+	psVectorStats (statsS, dASubset, NULL, NULL, 0xff);
+	psVectorStats (statsM, dMSubset, NULL, NULL, 0xff);
+
+	dSo->data.F32[i] = statsS->robustStdev;
+	dMo->data.F32[i] = statsM->sampleMean;
+
+	dRo->data.F32[i] = statsS->robustStdev / statsM->sampleMean;
+    }
+    psFree (dMSubset);
+    psFree (dASubset);
+    
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+    psVectorStats (stats, dRo, NULL, NULL, 0);
+
+    *errorScale = stats->sampleMedian;
+    *errorFloor = dSo->data.F32[0];
+
+    psFree (stats);
+    psFree (index);
+
+    psFree (dRo);
+    psFree (dMo);
+    psFree (dSo);
+    
+    psFree (statsS);
+    psFree (statsM);
+
+    return true;
+}
+
+bool psphotApResidTrend (pmReadout *readout, pmPSF *psf, int Npsf, int scale, float *errorScale, float *errorFloor, psVector *mask, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag) {
+
+    int Nx, Ny;
+	
+    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    stats->min = 2.0;
+    stats->max = 3.0;
+
+    if (readout->image->numCols > readout->image->numRows) {
+	Nx = scale;
+	Ny = PS_MAX (1, Nx * (readout->image->numRows / readout->image->numCols));
+    } else {
+	Ny = scale;
+	Nx = PS_MAX (1, Ny * (readout->image->numCols / readout->image->numRows));
+    }
+
+    if (Npsf < 3*Nx*Ny) {
+	return false;
+    }
+
+    // measure Trend2D for the current spatial scale
+    psFree (psf->ApTrend);
+    psf->ApTrend = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats);
+    pmTrend2DFit (psf->ApTrend, mask, 0xff, xPos, yPos, apResid, dMag);
+    
+    // construct the fitted values and the residuals
+    psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, xPos, yPos);
+    psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit);
+
+    // measure systematic errorFloor & systematic / photon scale factor
+    // XXX this is a bit arbitrary, but it forces ~3 stars from the bright bin per spatial bin
+    int nGroup = PS_MAX (3*Nx*Ny, 10);
+    psphotMagErrorScale (errorScale, errorFloor, dMag, apResidRes, nGroup);
+
+    psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup);
+    psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic error / photon error: %f\n", *errorScale);
+    psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic scatter floor: %f\n", *errorFloor);
+
+    psFree (stats);
+    psFree (apResidFit);
+    psFree (apResidRes);
+
+    return true;
+}
+