Index: /branches/eam_branches/20091113/psModules/src/astrom/pmAstrometryObjects.c
===================================================================
--- /branches/eam_branches/20091113/psModules/src/astrom/pmAstrometryObjects.c	(revision 26229)
+++ /branches/eam_branches/20091113/psModules/src/astrom/pmAstrometryObjects.c	(revision 26230)
@@ -37,4 +37,7 @@
 #include "pmKapaPlots.h"
 #include "pmAstrometryVisual.h"
+
+// XXX this is defined in pmPSFtry.h, which makes no sense
+float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction);
 
 #define PM_ASTROMETRYOBJECTS_DEBUG 1
@@ -283,5 +286,6 @@
             return results;
         }
-        psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", results->xStats->clippedMean, results->xStats->clippedStdev, results->xStats->clippedNvalues, x->n);
+        // psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", results->xStats->clippedMean, results->xStats->clippedStdev, results->xStats->clippedNvalues, x->n);
+        psTrace ("psModules.astrom", 3, "x resid: %f +/- %f (%ld of %ld)\n", results->xStats->robustMedian, results->xStats->robustStdev, results->xStats->clippedNvalues, x->n);
 
         if (!psVectorClipFitPolynomial2D (map->y, results->yStats, mask, 0xff, y, wt, X, Y)) {
@@ -296,18 +300,73 @@
             return results;
         }
-        psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", results->yStats->clippedMean, results->yStats->clippedStdev, results->yStats->clippedNvalues, y->n);
+        // psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", results->yStats->clippedMean, results->yStats->clippedStdev, results->yStats->clippedNvalues, y->n);
+        psTrace ("psModules.astrom", 3, "y resid: %f +/- %f (%ld of %ld)\n", results->yStats->robustMedian, results->yStats->robustStdev, results->yStats->clippedNvalues, y->n);
     }
     results->xStats->clipIter = stats->clipIter;
     results->yStats->clipIter = stats->clipIter;
 
-# if (0)
-    // XXX calculate the 90%-ile and the systematic scatter for each direction.
+    // *** calculate the 90%-ile and the systematic scatter for each direction.
+
+    // generate the X residual vector
     psVector *xFit = psPolynomial2DEvalVector (map->x, X, Y);
     if (!xFit) abort();
-    psVector *xRes = psBinaryOp (x, "-", xFit);
+    psVector *xRes = (psVector *) psBinaryOp (NULL, x, "-", xFit);
     if (!xRes) abort();
-
-    float dXsys = psVectorSystematicError (xRes, xErr, 0.05);
-# endif    
+    psFree (xFit);
+    
+    psVector *yFit = psPolynomial2DEvalVector (map->y, X, Y);
+    if (!yFit) abort();
+    psVector *yRes = (psVector *) psBinaryOp (NULL, y, "-", yFit);
+    if (!yRes) abort();
+    psFree (yFit);
+
+    // extract a high-quality subset (unmasked, S/N > XXX) and position errors
+    // XXX for now, generate a position error based on the magnitude error 
+    psVector *xErr     = psVectorAllocEmpty (match->n, PS_TYPE_F32);
+    psVector *yErr     = psVectorAllocEmpty (match->n, PS_TYPE_F32);
+    psVector *xResGood = psVectorAllocEmpty (match->n, PS_TYPE_F32);
+    psVector *yResGood = psVectorAllocEmpty (match->n, PS_TYPE_F32);
+
+    for (int i = 0; i < match->n; i++) {
+	if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
+        pmAstromMatch *pair = match->data[i];
+        pmAstromObj *rawStar = raw->data[pair->raw];
+        if (!isfinite(rawStar->dMag)) continue;
+        if (rawStar->dMag > 0.02) continue;
+
+	// two likely failure values: NAN or 0.0 --> use dMag in this case
+	float xErrValue, yErrValue;
+	if (isfinite(rawStar->chip->xErr) && (rawStar->chip->xErr > 0.0)) {
+	    xErrValue = rawStar->chip->xErr;
+	} else {
+	    xErrValue = PS_MAX(0.005, rawStar->dMag);
+	}
+	if (isfinite(rawStar->chip->yErr) && (rawStar->chip->yErr > 0.0)) {
+	    yErrValue = rawStar->chip->yErr;
+	} else {
+	    yErrValue = PS_MAX(0.005, rawStar->dMag);
+	}
+
+	psVectorAppend (xErr,     xErrValue);
+	psVectorAppend (yErr,     yErrValue);
+	psVectorAppend (xResGood, xRes->data.F32[i]);
+	psVectorAppend (yResGood, yRes->data.F32[i]);
+    }
+
+    results->dXsys = psVectorSystematicError (xResGood, xErr, 0.05);
+    results->dYsys = psVectorSystematicError (yResGood, yErr, 0.05);
+
+    results->dXrange = pmAstromVectorRange (xResGood, 0.1, 0.9, results->xStats->clippedStdev);
+    results->dYrange = pmAstromVectorRange (yResGood, 0.1, 0.9, results->yStats->clippedStdev);
+
+    psTrace ("psModules.astrom", 3, "dXsys: %f, dXrange: %f\n", results->dXsys, results->dXrange);
+    psTrace ("psModules.astrom", 3, "dYsys: %f, dYrange: %f\n", results->dYsys, results->dYrange);
+
+    psFree (xErr);
+    psFree (yErr);
+    psFree (xRes);
+    psFree (yRes);
+    psFree (xResGood);
+    psFree (yResGood);
 
     psFree (x);
@@ -321,4 +380,138 @@
 }
 
+// set the bin closest to the corresponding value.  if USE_END is +/- 1,
+// out-of-range saturates on lower/upper bin REGARDLESS of actual value
+#define PS_BIN_FOR_VALUE(RESULT, VECTOR, VALUE, USE_END) { \
+        psVectorBinaryDisectResult result; \
+        psScalar tmpScalar; \
+        tmpScalar.type.type = PS_TYPE_F32; \
+        tmpScalar.data.F32 = (VALUE); \
+        RESULT = psVectorBinaryDisect (&result, VECTOR, &tmpScalar); \
+        switch (result) { \
+          case PS_BINARY_DISECT_PASS: \
+            break; \
+          case PS_BINARY_DISECT_OUTSIDE_RANGE: \
+            psTrace("psModules.astrom", 6, "selected bin outside range"); \
+            if (USE_END == -1) { RESULT = 0; } \
+            if (USE_END == +1) { RESULT = VECTOR->n - 1; } \
+            break; \
+          case PS_BINARY_DISECT_INVALID_INPUT: \
+          case PS_BINARY_DISECT_INVALID_TYPE: \
+            psAbort ("programming error"); \
+            break; \
+        } }
+
+# define PS_BIN_INTERPOLATE(RESULT, VECTOR, BOUNDS, BIN, VALUE) { \
+        float dX, dY, Xo, Yo, Xt; \
+        if (BIN == BOUNDS->n - 1) { \
+            dX = 0.5*(BOUNDS->data.F32[BIN+1] - BOUNDS->data.F32[BIN-1]); \
+            dY = VECTOR->data.F32[BIN] - VECTOR->data.F32[BIN-1]; \
+            Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \
+            Yo = VECTOR->data.F32[BIN]; \
+        } else { \
+            dX = 0.5*(BOUNDS->data.F32[BIN+2] - BOUNDS->data.F32[BIN]); \
+            dY = VECTOR->data.F32[BIN+1] - VECTOR->data.F32[BIN]; \
+            Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \
+            Yo = VECTOR->data.F32[BIN]; \
+        } \
+        if (dY != 0) { \
+            Xt = (VALUE - Yo)*dX/dY + Xo; \
+        } else { \
+            Xt = Xo; \
+        } \
+        Xt = PS_MIN (BOUNDS->data.F32[BIN+1], PS_MAX(BOUNDS->data.F32[BIN], Xt)); \
+        psTrace("psModules.astrom", 6, "(Xo, Yo, dX, dY, Xt, Yt) is (%.2f %.2f %.2f %.2f %.2f %.2f)\n", \
+                Xo, Yo, dX, dY, Xt, VALUE); \
+        RESULT = Xt; }
+
+float pmAstromVectorRange (psVector *myVector, float minFrac, float maxFrac, float stdevGuess) {
+
+    psStats *stats = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
+    psHistogram *histogram = NULL;      // Histogram of the data
+    psHistogram *cumulative = NULL;     // Cumulative histogram of the data
+    float min = NAN, max = NAN;         // Mimimum and maximum values
+
+    // Get the minimum and maximum values
+    if (!psVectorStats(stats, myVector, NULL, NULL, 0)) {
+	psFree(stats);
+	return NAN;
+    }
+    min = stats->min;
+    max = stats->max;
+    if (isnan(min) || isnan(max)) {
+	psFree(stats);
+	return NAN;
+    }
+
+    psTrace("psModules.astrom", 5, "Data min/max is (%.2f, %.2f)\n", min, max);
+
+    // If all data points have the same value, then we set the appropriate members of stats and return.
+    if (fabs(max - min) <= FLT_EPSILON) {
+	psFree (stats);
+	return 0.0;
+    }
+    
+    // Define the histogram bin size.  
+    float binSize = 0.001;
+    long numBins = PS_MIN(100000, (max - min) / binSize); // Number of bins
+    psTrace("psModules.astrom", 5, "Numbins is %ld\n", numBins);
+    psTrace("psModules.astrom", 5, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max);
+
+    // allocate the histogram containers
+    histogram = psHistogramAlloc(min, max, numBins);
+    cumulative = psHistogramAlloc(min, max, numBins);
+
+    if (!psVectorHistogram(histogram, myVector, NULL, NULL, 0)) {
+	// if psVectorHistogram returns false, we have a programming error
+	psAbort ("Unable to generate histogram");
+    }
+    if (psTraceGetLevel("psModules.astrom") >= 8) {
+	PS_VECTOR_PRINT_F32(histogram->bounds);
+	PS_VECTOR_PRINT_F32(histogram->nums);
+    }
+
+    // Convert the specific histogram to a cumulative histogram
+    // The cumulative histogram data points correspond to the UPPER bound value (N < Bin[i+1])
+    cumulative->nums->data.F32[0] = histogram->nums->data.F32[0];
+    for (long i = 1; i < histogram->nums->n; i++) {
+	cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i];
+	cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i];
+    }
+    if (psTraceGetLevel("psModules.astrom") >= 8) {
+	PS_VECTOR_PRINT_F32(cumulative->bounds);
+	PS_VECTOR_PRINT_F32(cumulative->nums);
+    }
+
+    // Find the bin which contains the first data point above the limit
+    long totalDataPoints = cumulative->nums->data.F32[numBins - 1];
+    psTrace("psModules.astrom", 6, "Total data points is %ld\n", totalDataPoints);
+
+    // find bin which is the lower bound of the limit value (value[bin] < f < value[bin+1]
+    long binMin;
+    PS_BIN_FOR_VALUE(binMin, cumulative->nums, minFrac * totalDataPoints, 0);
+    psTrace("psModules.astrom", 6, "The bin is %ld (%.4f to %.4f)\n", binMin, cumulative->bounds->data.F32[binMin], cumulative->bounds->data.F32[binMin+1]);
+
+    // Linear interpolation to the limit value in bin units
+    float valueMin;
+    PS_BIN_INTERPOLATE (valueMin, cumulative->nums, cumulative->bounds, binMin, totalDataPoints * minFrac);
+    psTrace("psModules.astrom", 6, "limit value is %f\n", valueMin);
+
+    // find bin which is the lower bound of the limit value (value[bin] < f < value[bin+1]
+    long binMax;
+    PS_BIN_FOR_VALUE(binMax, cumulative->nums, maxFrac * totalDataPoints, 0);
+    psTrace("psModules.astrom", 6, "The bin is %ld (%.4f to %.4f)\n", binMax, cumulative->bounds->data.F32[binMax], cumulative->bounds->data.F32[binMax+1]);
+
+    // Linear interpolation to the limit value in bin units
+    float valueMax;
+    PS_BIN_INTERPOLATE (valueMax, cumulative->nums, cumulative->bounds, binMax, totalDataPoints * maxFrac);
+    psTrace("psModules.astrom", 6, "limit value is %f\n", valueMax);
+
+    // Clean up
+    psFree(histogram);
+    psFree(cumulative);
+    psFree(stats);
+
+    return (valueMax - valueMin);
+}
 
 /******************************************************************************
@@ -1057,5 +1250,5 @@
     }
 
-    psLogMsg ("psastro", 3, "generate unique matches to reference stars: %ld matches -> %ld matches\n", matches->n, unique->n);
+    psLogMsg ("psModules.astrom", 3, "generate unique matches to reference stars: %ld matches -> %ld matches\n", matches->n, unique->n);
     psFree (refstarMatches);
 
Index: /branches/eam_branches/20091113/psModules/src/astrom/pmAstrometryObjects.h
===================================================================
--- /branches/eam_branches/20091113/psModules/src/astrom/pmAstrometryObjects.h	(revision 26229)
+++ /branches/eam_branches/20091113/psModules/src/astrom/pmAstrometryObjects.h	(revision 26230)
@@ -96,4 +96,8 @@
     int     nMatch;                     ///<
     double  nSigma;                     ///<
+    double  dXsys;			///< systematic error in X
+    double  dYsys;			///< systematic error in Y
+    double  dXrange;			///< 10% - 90% range X residuals (unmasked, high S/N)
+    double  dYrange;			///< 10% - 90% range Y residuals (unmasked, high S/N)
 }
 pmAstromFitResults;
@@ -355,4 +359,6 @@
 );
 
+float pmAstromVectorRange (psVector *myVector, float minFrac, float maxFrac, float stdevGuess);
+
 /// @}
 #endif // PM_ASTROMETRY_OBJECTS_H