Index: trunk/psModules/src/objects/pmPSFtry.c
===================================================================
--- trunk/psModules/src/objects/pmPSFtry.c	(revision 19474)
+++ trunk/psModules/src/objects/pmPSFtry.c	(revision 19906)
@@ -5,6 +5,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.61 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-09-11 00:38:23 $
+ *  @version $Revision: 1.62 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-10-06 13:05:13 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -150,5 +150,5 @@
     }
     psLogMsg ("psphot.psftry", 4, "fit ext:   %f sec for %d of %ld sources\n", psTimerMark ("fit"), Next, sources->n);
-    psTrace ("psphot.psftry", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n);
+    psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n);
 
     if (Next == 0) {
@@ -215,5 +215,5 @@
         psfTry->metricErr->data.F32[i] = source->errMag;
 
-        psTrace ("psphot.psftry", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
+        psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
         Npsf ++;
     }
@@ -221,5 +221,5 @@
 
     psLogMsg ("psphot.psftry", 4, "fit psf:   %f sec for %d of %ld sources\n", psTimerMark ("fit"), Npsf, sources->n);
-    psTrace ("psphot.psftry", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n);
+    psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n);
 
     if (Npsf == 0) {
@@ -568,5 +568,5 @@
         // check the fit residuals and increase Nx,Ny until the error is minimized
         // pmPSFParamTrend increases the number along the longer of x or y
-        for (int i = 1; i < PS_MAX (psf->trendNx, psf->trendNy); i++) {
+        for (int i = 1; i <= PS_MAX (psf->trendNx, psf->trendNy); i++) {
 
             psVectorInit (mask, 0);
@@ -578,5 +578,7 @@
             // use the bright end scatter from the constant fit to set the scale for the
             // versions with 2D variation.  potentially scale by poisson errors...
-            if (i == 1) {
+            // if (i == 1) {
+	    // XXX let's drop this for the moment: relies on valid result from errorFloor
+	    if (0) {
                 // if (psf->poissonErrorsParams) do something else..
                 dz = psVectorAlloc (sources->n, PS_TYPE_F32);
@@ -629,12 +631,30 @@
         for (int i = 0; i < psf->psfTrendStats->clipIter; i++) {
 
+	    psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
+	    psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
+
+	    pmTrend2D *trend = NULL;
+	    float mean, stdev;
+
             // XXX we are using the same stats structure on each pass: do we need to re-init it?
             bool status = true;
-            status &= pmTrend2DFit (psf->params->data[PM_PAR_E0], srcMask, 0xff, x, y, e0, NULL);
-            psTrace ("psModules.pmPSFtry", 4, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0->n);
-            status &= pmTrend2DFit (psf->params->data[PM_PAR_E1], srcMask, 0xff, x, y, e1, NULL);
-            psTrace ("psModules.pmPSFtry", 4, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1->n);
-            status &= pmTrend2DFit (psf->params->data[PM_PAR_E2], srcMask, 0xff, x, y, e2, NULL);
-            psTrace ("psModules.pmPSFtry", 4, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2->n);
+
+	    trend = psf->params->data[PM_PAR_E0];
+            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL);
+	    mean = psStatsGetValue (trend->stats, meanOption);
+	    stdev = psStatsGetValue (trend->stats, stdevOption);
+            psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n);
+
+	    trend = psf->params->data[PM_PAR_E1];
+            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL);
+	    mean = psStatsGetValue (trend->stats, meanOption);
+	    stdev = psStatsGetValue (trend->stats, stdevOption);
+            psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n);
+
+	    trend = psf->params->data[PM_PAR_E2];
+            status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL);
+	    mean = psStatsGetValue (trend->stats, meanOption);
+	    stdev = psStatsGetValue (trend->stats, stdevOption);
+            psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n);
 
             if (!status) {
@@ -676,8 +696,12 @@
     if (psf->fieldNx > psf->fieldNy) {
         Nx = scale;
-        Ny = PS_MAX (1, Nx * (psf->fieldNy / psf->fieldNx));
+	float AR = psf->fieldNy / (float) psf->fieldNx;
+	Ny = (int) (Nx * AR + 0.5);
+        Ny = PS_MAX (1, Ny);
     } else {
         Ny = scale;
-        Nx = PS_MAX (1, Ny * (psf->fieldNx / psf->fieldNy));
+	float AR = psf->fieldNx / (float) psf->fieldNy;
+	Nx = (int) (Ny * AR + 0.5);
+        Nx = PS_MAX (1, Nx);
     }
 
@@ -716,10 +740,38 @@
     for (int i = 0; i < nIter; i++) {
         // XXX we are using the same stats structure on each pass: do we need to re-init it?
-        pmTrend2DFit (psf->params->data[PM_PAR_E0], mask, 0xff, x, y, e0obs, dz);
-        psTrace ("psModules.objects", 5, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0obs->n);
-        pmTrend2DFit (psf->params->data[PM_PAR_E1], mask, 0xff, x, y, e1obs, dz);
-        psTrace ("psModules.objects", 5, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1obs->n);
-        pmTrend2DFit (psf->params->data[PM_PAR_E2], mask, 0xff, x, y, e2obs, dz);
-        psTrace ("psModules.objects", 5, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2obs->n);
+        // pmTrend2DFit (psf->params->data[PM_PAR_E0], mask, 0xff, x, y, e0obs, dz);
+        // psTrace ("psModules.objects", 5, "clipped E0 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e0obs->n);
+        // pmTrend2DFit (psf->params->data[PM_PAR_E1], mask, 0xff, x, y, e1obs, dz);
+        // psTrace ("psModules.objects", 5, "clipped E1 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e1obs->n);
+        // pmTrend2DFit (psf->params->data[PM_PAR_E2], mask, 0xff, x, y, e2obs, dz);
+        // psTrace ("psModules.objects", 5, "clipped E2 : keeping %ld of %ld\n", psf->psfTrendStats->clippedNvalues, e2obs->n);
+
+	psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options);
+	psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options);
+
+	pmTrend2D *trend = NULL;
+	float mean, stdev;
+
+	// XXX we are using the same stats structure on each pass: do we need to re-init it?
+	bool status = true;
+
+	trend = psf->params->data[PM_PAR_E0];
+	status &= pmTrend2DFit (trend, mask, 0xff, x, y, e0obs, dz);
+	mean = psStatsGetValue (trend->stats, meanOption);
+	stdev = psStatsGetValue (trend->stats, stdevOption);
+	psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs->n);
+
+	trend = psf->params->data[PM_PAR_E1];
+	status &= pmTrend2DFit (trend, mask, 0xff, x, y, e1obs, dz);
+	mean = psStatsGetValue (trend->stats, meanOption);
+	stdev = psStatsGetValue (trend->stats, stdevOption);
+	psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs->n);
+
+	trend = psf->params->data[PM_PAR_E2];
+	status &= pmTrend2DFit (trend, mask, 0xff, x, y, e2obs, dz);
+	mean = psStatsGetValue (trend->stats, meanOption);
+	stdev = psStatsGetValue (trend->stats, stdevOption);
+	psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n);
+
     }
     psf->psfTrendStats->clipIter = nIter; // restore default setting
@@ -734,6 +786,33 @@
     psVector *e2res = (psVector *) psBinaryOp (NULL, (void *) e2obs, "-", (void *) e2fit);
 
+// XXX this code determines the formal error on the map, not the scatter
+# if (0)
     // measure errors on the map pixels
     pmTrend2D *trend;
+    psStatsOptions meanOption;
+    psStatsOptions stdevOption;
+    float mapErrorSum = 0.0;
+    float mapError = 0.0;
+
+    trend = psf->params->data[PM_PAR_E0];
+    meanOption = psStatsMeanOption (trend->stats->options);
+    stdevOption = psStatsStdevOption (trend->stats->options);
+    mapError = PS_SQR(psStatsGetValue (trend->stats, stdevOption));
+    mapErrorSum += mapError;
+    psTrace ("psModules.objects", 4, "E0 error: %f\n", sqrt(mapError));
+
+    trend = psf->params->data[PM_PAR_E1];
+    meanOption = psStatsMeanOption (trend->stats->options);
+    stdevOption = psStatsStdevOption (trend->stats->options);
+    mapError = PS_SQR(psStatsGetValue (trend->stats, stdevOption));
+    mapErrorSum += mapError;
+    psTrace ("psModules.objects", 4, "E1 error: %f\n", sqrt(mapError));
+
+    trend = psf->params->data[PM_PAR_E2];
+    meanOption = psStatsMeanOption (trend->stats->options);
+    stdevOption = psStatsStdevOption (trend->stats->options);
+    mapError = PS_SQR(psStatsGetValue (trend->stats, stdevOption));
+    mapErrorSum += mapError;
+    psTrace ("psModules.objects", 4, "E2 error: %f\n", sqrt(mapError));
 
     trend = psf->params->data[PM_PAR_E0];
@@ -769,4 +848,5 @@
 
     psTrace ("psModules.objects", PS_LOG_INFO, "total map error: %f\n", sqrt(mapErrorSum));
+# endif
 
     // measure systematic errorFloor & systematic / photon scale factor
@@ -776,5 +856,6 @@
                             psStatsStdevOption(psf->psfTrendStats->options));
 
-    *errorTotal = sqrt(PS_SQR(*errorFloor) + mapErrorSum);
+    // XXX Bogus: *errorTotal = sqrt(PS_SQR(*errorFloor) + mapErrorSum);
+    *errorTotal = *errorFloor;
 
     psLogMsg ("psphot.psftry", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup);
@@ -815,4 +896,5 @@
     for (int i = 0; i < nBin; i++) {
         int j;
+	int nValid = 0;
         for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) {
             int N = index->data.U32[n];
@@ -822,5 +904,8 @@
 
             mkSubset->data.U8[j]   = mask->data.U8[N];
-        }
+	    if (!mask->data.U8[N]) nValid ++;
+        }
+	if (nValid < 3) continue;
+
         dE0subset->n = j;
         dE1subset->n = j;
@@ -849,4 +934,6 @@
         }
     }
+    *errorFloor = min;
+
     psFree (dE0subset);
     psFree (dE1subset);