Index: trunk/psphot/src/psphotApResid.c
===================================================================
--- trunk/psphot/src/psphotApResid.c	(revision 5952)
+++ trunk/psphot/src/psphotApResid.c	(revision 5980)
@@ -1,5 +1,114 @@
 # include "psphot.h"
 
-bool psphotApResid (psArray *sources, psMetadata *config, pmPSF *psf) { 
+psPolynomial4D *psVectorChiClipFitPolynomial4D(
+    psPolynomial4D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z,
+    const psVector *t)
+{
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_PTR_NON_NULL(stats, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(z, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(t, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(t, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(fErr, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL);
+
+    // clipping range defined by min and max and/or clipSigma
+    float minClipSigma;
+    float maxClipSigma;
+    if (isfinite(stats->max)) {
+        maxClipSigma = +fabs(stats->max);
+    } else {
+        maxClipSigma = +fabs(stats->clipSigma);
+    }
+    if (isfinite(stats->min)) {
+        minClipSigma = -fabs(stats->min);
+    } else {
+        minClipSigma = -fabs(stats->clipSigma);
+    }
+    psVector *fit   = NULL;
+    psVector *resid = psVectorAlloc (x->n, PS_TYPE_F64);
+
+    // eventual expansion: user supplies one of various stats option pairs,
+    // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to
+    // evaluate the clipping sigma
+    // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
+    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    for (int N = 0; N < stats->clipIter; N++) {
+        int Nkeep = 0;
+
+        poly = psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t);
+        fit = psPolynomial4DEvalVector (poly, x, y, z, t);
+        resid = (psVector *) psBinaryOp (resid, (void *) f, "-", (void *) fit);
+
+        stats  = psVectorStats (stats, resid, NULL, mask, maskValue);
+        psTrace (".psphot.VectorClipFit", 5, "resid stats: %f +/- %f\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 < resid->n; i++) {
+            if ((mask != NULL) && (mask->data.U8[i] & maskValue)) {
+                continue;
+            }
+	    float sigma = hypot (psVectorGet (fErr, i), stats->sampleStdev);
+            if (resid->data.F64[i] - stats->sampleMedian > sigma*maxClipSigma) {
+                if (mask != NULL) {
+                    mask->data.U8[i] |= 0x01;
+                }
+                continue;
+            }
+            if (resid->data.F64[i] - stats->sampleMedian < sigma*minClipSigma) {
+                if (mask != NULL) {
+                    mask->data.U8[i] |= 0x01;
+                }
+                continue;
+            }
+            Nkeep ++;
+        }
+
+        psTrace (".psphot.VectorClipFit", 4, "keeping %d of %d pts for fit\n",
+                 Nkeep, x->n);
+
+        stats->clippedNvalues = Nkeep;
+        psFree (fit);
+    }
+    // Free local temporary variables
+    psFree (resid);
+
+    if (poly == NULL) {
+        psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data.  Returning NULL.\n");
+        return(NULL);
+    }
+    return(poly);
+}
+
+// measure the aperture residual statistics
+bool psphotApResid (eamReadout *imdata, psArray *sources, psMetadata *config, pmPSF *psf) { 
 
     int Npsf;
@@ -8,22 +117,23 @@
     pmSource *source;
 
-    // XXX EAM : check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS
-    float RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS");
-
     psTimerStart ("psphot");
 
-    // XXX EAM : drop this if we know the list is sorted
-    sources = psArraySort (sources, psphotSortBySN);
-
+    // measure the aperture loss as a function of radius for PSF
+    float REF_RADIUS = psMetadataLookupF32 (&status, config, "PSF_REF_RADIUS");
+    psf->growth = pmGrowthCurveAlloc (3.0, REF_RADIUS, 0.1);
+    psphotGrowthCurve (imdata, psf);
+    
     psVector *mask    = psVectorAlloc (300, PS_TYPE_U8);
     psVector *xPos    = psVectorAlloc (300, PS_TYPE_F64);
     psVector *yPos    = psVectorAlloc (300, PS_TYPE_F64);
-    psVector *rflux   = psVectorAlloc (300, PS_TYPE_F64);
+    psVector *flux    = psVectorAlloc (300, PS_TYPE_F64);
+    psVector *r2rflux = psVectorAlloc (300, PS_TYPE_F64);
     psVector *apResid = psVectorAlloc (300, PS_TYPE_F64);
-    mask->n = xPos->n = yPos->n = rflux->n = apResid->n = 0;
+    psVector *dMag    = psVectorAlloc (300, PS_TYPE_F64);
+    mask->n = xPos->n = yPos->n = flux->n = r2rflux->n = apResid->n = dMag->n = 0;
     Npsf = 0;
 
-    // select the NNN brightest, non-saturated sources, or just select PSFSTARs?
-    for (int i = 0; (i < sources->n) && (Npsf < 300); i++) {
+    // select all good PM_SOURCE_STAR entries
+    for (int i = 0; i < sources->n; i++) {
 	source = sources->data[i];
 
@@ -34,6 +144,6 @@
 	if (source->mode &  PM_SOURCE_POOR) continue;
 
-	// get magnitudes, uncorrected for (x, y, rflux)
-	model = pmSourceMagnitudes (source, NULL, RADIUS);
+	// get uncorrected magnitudes in scaled apertures
+	model = pmSourceMagnitudes (source, NULL, 0);
 	if (model == NULL) continue;
 
@@ -41,45 +151,199 @@
 	xPos->data.F64[Npsf] = model->params->data.F32[2];
 	yPos->data.F64[Npsf] = model->params->data.F32[3];
-	rflux->data.F64[Npsf] = pow(10.0, 0.4*source->fitMag);
-	apResid->data.F64[Npsf] = source->apMag - source->fitMag;
-
-	psVectorExtend (mask, 100, 1);
-	psVectorExtend (xPos, 100, 1);
-	psVectorExtend (yPos, 100, 1);
-	psVectorExtend (rflux, 100, 1);
+
+	flux->data.F64[Npsf] = pow(10.0, -0.4*source->fitMag);
+	r2rflux->data.F64[Npsf] = PS_SQR(model->radius) / flux->data.F64[Npsf];
+	
+	apResid->data.F64[Npsf] = source->apMag + pmGrowthCurveCorrect (psf->growth, model->radius) - source->fitMag ;
+
+	// XXX sanity clip?
+	// XXX need to see if all data were tossed?
+	// XXX need to subtract median?
+	if (fabs(apResid->data.F64[Npsf]) > 0.2) continue;
+
+	dMag->data.F64[Npsf] = model->dparams->data.F32[1] / model->params->data.F32[1];
+
+	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 ++;
     }
-    psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec\n", psTimerMark ("psphot"));
+    psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec for %d objects\n", psTimerMark ("psphot"), Npsf);
+
+    // APTREND options : NONE SKYBIAS XY_LIN XY_QUAD SKY_XY_LIN FULL 
+    char *ApTrendOption = psMetadataLookupPtr (&status, config, "APTREND");
+    if (!status) ApTrendOption = psStringCopy ("SKYBIAS");
 
     // 3hi/1lo sigma clipping on the rflux vs metric fit
     psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
-    stats->min = 1.0;
+    stats->min = 3.0;
     stats->max = 3.0;
-    stats->clipIter = 3;
-
-    // first clip out objects which are too far from the median 
-    stats->clipIter = 1;
-    maskToConstant (psf->ApTrend);
-    psf->ApTrend  = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, xPos, yPos, rflux);
-
-    // next, fit just SkyBias and clip out objects which are too far from the median 
-    stats->clipIter = 2;
-    maskToSkyBias (psf->ApTrend);
-    psf->ApTrend  = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, xPos, yPos, rflux);
-
-    // finally, fit x, y, SkyBias and clip out objects which are too far from the median 
-    stats->clipIter = 2;
-    maskToDefault (psf->ApTrend);
-    psf->ApTrend  = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, xPos, yPos, rflux);
-
-    // linear clipped fit of apResid to rflux, xPos, yPos
-    # if (0)
-    psf->ApTrend  = psVectorClipFitPolynomial3D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, NULL, xPos, yPos, rflux);
-    psf->skyBias  = psf->ApTrend->coeff[0][0][1] / (M_PI * PS_SQR(RADIUS));
-    psf->ApResid  = psf->ApTrend->coeff[0][0][0];
-    psf->dApResid = stats->sampleStdev;
-    psf->ApTrend->coeff[0][0][1] = 0;
-    # endif
+
+    // constant only
+    if (!strcasecmp (ApTrendOption, "CONSTANT")) {
+	stats->clipIter = 2;
+	pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
+	psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+
+	stats->clipIter = 3;
+	pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
+	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+    }
+
+    // constant and skybias only
+    if (!strcasecmp (ApTrendOption, "SKYBIAS")) {
+	// first clip out objects which are too far from the median 
+	stats->clipIter = 2;
+	pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
+	psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+
+	// apply the fit
+	stats->clipIter = 3;
+	pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS);
+	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+    }
+
+    if (!strcasecmp (ApTrendOption, "SKYSAT")) {
+	// first clip out objects which are too far from the median 
+	stats->clipIter = 2;
+	pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
+	psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+
+	// apply the fit
+	stats->clipIter = 2;
+	pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS);
+	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+
+	// apply the fit
+	stats->clipIter = 3;
+	pmPSF_MaskApTrend (psf, PM_PSF_SKYSAT);
+	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+    }
+
+    // constant and linear X,Y only
+    if (!strcasecmp (ApTrendOption, "XY_LIN")) {
+	// first clip out objects which are too far from the median 
+	stats->clipIter = 2;
+	pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
+	psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+
+	// apply the fit
+	stats->clipIter = 3;
+	pmPSF_MaskApTrend (psf, PM_PSF_XY_LIN);
+	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+    }
+
+    // constant and quadratic X,Y only
+    if (!strcasecmp (ApTrendOption, "XY_QUAD")) {
+	// first clip out objects which are too far from the median 
+	stats->clipIter = 2;
+	pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
+	psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+
+	// apply the fit
+	stats->clipIter = 3;
+	pmPSF_MaskApTrend (psf, PM_PSF_XY_QUAD);
+	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+    }
+
+     // constant and sky, linear X,Y only
+    if (!strcasecmp (ApTrendOption, "SKY_XY_LIN")) {
+	// first clip out objects which are too far from the median 
+	stats->clipIter = 2;
+	pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
+	psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+
+	// apply the fit
+	stats->clipIter = 3;
+	pmPSF_MaskApTrend (psf, PM_PSF_SKY_XY_LIN);
+	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+    }
+
+     // constant and sky, linear X,Y only
+    if (!strcasecmp (ApTrendOption, "SKYSAT_XY_LIN")) {
+	// first clip out objects which are too far from the median 
+	stats->clipIter = 2;
+	pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
+	psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+
+	// apply the fit
+	stats->clipIter = 3;
+	pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS);
+	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+
+	// apply the fit
+	stats->clipIter = 3;
+	pmPSF_MaskApTrend (psf, PM_PSF_SKYSAT_XY_LIN);
+	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+    }
+
+    if (!strcasecmp (ApTrendOption, "ALL")) {
+	// first clip out objects which are too far from the median 
+	stats->clipIter = 2;
+	pmPSF_MaskApTrend (psf, PM_PSF_CONSTANT);
+	psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+
+	// fit just SkyBias and clip out objects which are too far from the median 
+	stats->clipIter = 2;
+	pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS);
+	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+
+	// finally, fit x, y, SkyBias and clip out objects which are too far from the median 
+	stats->clipIter = 3;
+	pmPSF_MaskApTrend (psf, PM_PSF_ALL);
+	psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
+    }
+
+# if (1)
+    psPolynomial4D *poly = psf->ApTrend;
+    for (int nt = 0; nt <= poly->nT; nt++) {
+	for (int nz = 0; nz <= poly->nZ; nz++) {
+	    for (int ny = 0; ny <= poly->nY; ny++) {
+		for (int nx = 0; nx <= poly->nX; nx++) {
+		    if (poly->mask[nx][ny][nz][nt]) continue;
+		    fprintf (stderr, "%d %d %d %d : %22.15g\n", nx, ny, nz, nt, poly->coeff[nx][ny][nz][nt]);
+		}
+	    }
+	}
+    }
+# endif
+
+    // 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);
+
+# if (0)
+    FILE *fout = fopen ("resid.dat", "w");
+    for (int i = 0; i < resid->n; i++) {
+	fprintf (fout, "%d %f %f %f %f  %f %f %f %f %d\n", 
+		 i, 
+		 (float) psVectorGet(xPos, i), (float) psVectorGet(yPos, i), (float) psVectorGet(r2rflux, i), (float) psVectorGet(flux, i), 
+		 (float) psVectorGet(apResid, i), (float) psVectorGet(fit, i), (float) psVectorGet(resid, i), (float) psVectorGet(dMag, i),
+		 mask->data.U8[i]);
+    }
+    fclose (fout);
+# endif
+
+    // measure scatter for sources with dMag < 0.01 (S/N = 100)
+    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;
+	}
+    }
+    residStats  = psVectorStats (residStats, resid, NULL, mask, 0x03);
+
+    // 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;
+    psf->nApResid = residStats->clippedNvalues;
 
     /*
@@ -89,23 +353,19 @@
     */
 
-    # if (0)
-    psPolynomial3D *poly = psf->ApTrend;
-    for (int nz = 0; nz <= poly->nZ; nz++) {
-	for (int ny = 0; ny <= poly->nY; ny++) {
-	    for (int nx = 0; nx <= poly->nX; nx++) {
-		fprintf (stderr, "%d %d %d : %22.15g\n", nx, ny, nz, poly->coeff[nx][ny][nz]);
-	    }
-	}
-    }
-    # endif
-
     psLogMsg ("psphot.apresid", 3, "measure full-frame aperture residual: %f sec\n", psTimerMark ("psphot"));
     psLogMsg ("psphot.apresid", 4, "aperture residual: %f +/- %f : %f bias\n", psf->ApResid, psf->dApResid, psf->skyBias);
-
-    psFree (stats);
-    psFree (mask);
-    psFree (rflux);
-    psFree (apResid);
+    psLogMsg ("psphot.apresid", 4, "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]);
+
+    // psFree (stats);
+    // psFree (mask);
+    // psFree (rflux);
+    // psFree (apResid);
 
     return true;
 }
+