Index: /branches/eam_branch_20070830/psphot/src/psphotApResid.c
===================================================================
--- /branches/eam_branch_20070830/psphot/src/psphotApResid.c	(revision 14785)
+++ /branches/eam_branch_20070830/psphot/src/psphotApResid.c	(revision 14786)
@@ -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,8 +51,8 @@
 
     psVector *mask    = psVectorAllocEmpty (300, PS_TYPE_U8);
-    psVector *xPos    = psVectorAllocEmpty (300, PS_TYPE_F64);
-    psVector *yPos    = psVectorAllocEmpty (300, PS_TYPE_F64);
-    psVector *apResid = psVectorAllocEmpty (300, PS_TYPE_F64);
-    psVector *dMag    = psVectorAllocEmpty (300, PS_TYPE_F64);
+    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;
 
@@ -97,16 +94,16 @@
         }
 
-        apResid->data.F64[Npsf] = dap
-        xPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_XPOS];
-        yPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_YPOS];
+        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];
+        dMag->data.F32[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], 
+		     xPos->data.F32[Npsf], yPos->data.F32[Npsf], 
+		     source->psfMag, dMag->data.F32[Npsf], apResid->data.F32[Npsf], 
 		     mask->data.U8[Npsf]);
 	}
@@ -135,212 +132,20 @@
 
     // XXX deprecating the old code which allowed the ApResid to be fitted as a function of flux and r^2/flux
-    // 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 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);
+    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN);
     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)
+    // XXX set the pmTrend2D mode, nXtrend, nYtrend based on the user inputs
+    // XXX psf->ApTrend = pmTrend2DAlloc (mode, readout->image->numCols, readout->image->numRows, nXtrend, nYtrend, 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
+    psf->ApTrend = pmTrend2DAlloc (PM_TREND_MAP, readout->image, 5, 5, stats);
+    pmTrend2DFit (psf->ApTrend, xPos, yPos, apResid, dMag);
+    
+    psphotSaveImage (NULL, psf->ApTrend->map->map, "apresid.map.fits");
 
     // construct the fitted values and the residuals
-    psVector *fit   = psPolynomial4DEvalVector (psf->ApTrend, xPos, yPos, r2rflux, flux);
+    psVector *fit   = pmTrend2DEvalVector (psf->ApTrend, xPos, yPos);
     psVector *resid = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) fit);
 
@@ -349,5 +154,5 @@
     psStats *residStats = psStatsAlloc (PS_STAT_SAMPLE_STDEV);
     for (int i = 0; i < dMag->n; i++) {
-        if (dMag->data.F64[i] > 0.01) {
+        if (dMag->data.F32[i] > 0.01) {
             mask->data.U8[i] |= 0x02;
         }
@@ -357,15 +162,11 @@
 
     // 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->ApResid  = pmTrend2DEval (psf->ApTrend, 0.0, 0.0); // XXX use chip center?
     psf->dApResid = residStats->sampleStdev;
-    psf->ApTrend->coeff[0][0][1][0] = 0;
-    psf->ApTrend->coeff[0][0][0][1] = 0;
     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);
@@ -373,20 +174,10 @@
     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 (mask);
     psFree (xPos);
     psFree (yPos);
-    psFree (flux);
-    psFree (r2rflux);
     psFree (apResid);
     psFree (dMag);