Index: trunk/psModules/src/astrom/pmAstrometryObjects.c
===================================================================
--- trunk/psModules/src/astrom/pmAstrometryObjects.c	(revision 24034)
+++ trunk/psModules/src/astrom/pmAstrometryObjects.c	(revision 26260)
@@ -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,8 +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;
+
+    // *** 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 = (psVector *) psBinaryOp (NULL, x, "-", xFit);
+    if (!xRes) abort();
+    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);
@@ -311,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);
+}
 
 /******************************************************************************
@@ -634,17 +837,17 @@
         }
 
-# if 0
-        char line[16];
-        psFits *fits = psFitsOpen ("grid.image.fits", "w");
-        psFitsWriteImage (fits, NULL, gridNP, 0, NULL);
-        psFitsClose (fits);
-        fprintf (stderr, "wrote grid image, press return to continue\n");
-        fgets (line, 15, stdin);
-# endif
+	if (psTraceGetLevel("psModules.astrom") >= 5) {
+	    char line[16];
+	    psFits *fits = psFitsOpen ("grid.image.fits", "w");
+	    psFitsWriteImage (fits, NULL, gridNP, 0, NULL);
+	    psFitsClose (fits);
+	    fprintf (stderr, "wrote grid image, press return to continue\n");
+	    fgets (line, 15, stdin);
+	}
 
         // only check bins with at least 1/2 of max bin
         // XXX requiring at least 3 matches in bin
         int minNpts = PS_MAX (0.5*imStats->max, 5);
-        psTrace("psModule.astrom", 5, "minNpts: %d, min: %d, max: %d, median: %f, stdev: %f", minNpts, (int)(imStats->min), (int)(imStats->max), imStats->sampleMedian, imStats->sampleStdev);
+        psTrace("psModule.astrom", 4, "minNpts: %d, min: %d, max: %d, median: %f, stdev: %f", minNpts, (int)(imStats->min), (int)(imStats->max), imStats->sampleMedian, imStats->sampleStdev);
 
         // find the 'best' bin
@@ -687,5 +890,6 @@
 
         // XXX this function is crashing
-        // pmAstromVisualPlotGridMatch(raw, ref, gridNP, stats->offset.x, stats->offset.y, maxOffpix, Scale, Offset);
+        pmAstromVisualPlotGridMatch(raw, ref, gridNP, stats->offset.x, stats->offset.y, maxOffpix, Scale, Offset);
+        pmAstromVisualPlotGridMatchOverlay(raw, ref);
 
         psFree (imStats);
@@ -962,2 +1166,91 @@
 */
 
+/*****************************************************************************/
+static void pmAstromMatchInfoFree (pmAstromMatchInfo *info)
+{
+    if (info == NULL) return;
+    return;
+}
+
+
+/*****************************************************************************/
+pmAstromMatchInfo *pmAstromMatchInfoAlloc()
+{
+    pmAstromMatchInfo *info = psAlloc (sizeof(pmAstromMatchInfo));
+    psMemSetDeallocator(info, (psFreeFunc) pmAstromMatchInfoFree);
+
+    info->match = NULL;
+    info->radius = NAN;
+
+    return (info);
+}
+
+// generate a unique set of matches (choose closest match)
+psArray *pmAstromRadiusMatchUniq (psArray *rawstars, psArray *refstars, psArray *matches) {
+
+    // I have the matches between the refstars and the rawstars.  
+    // For each refstar, find the single match which has the smallest radius and reject
+    // all others.  
+
+    // create an array of refstars->n arrays, each containing all of the matches for the
+    // given refstar.  
+
+    psArray *refstarMatches = psArrayAlloc (refstars->n);
+
+    for (int i = 0; i < matches->n; i++) {
+
+	pmAstromMatch *match = matches->data[i];
+	
+	// accumulate this refstar match on the array for this refstar (create if needed)
+	psArray *refSet = refstarMatches->data[match->ref];
+	if (!refSet) {
+	    refstarMatches->data[match->ref] = psArrayAllocEmpty (8);
+	    refSet = refstarMatches->data[match->ref];
+	}
+
+	pmAstromMatchInfo *matchInfo = pmAstromMatchInfoAlloc();
+
+	pmAstromObj *refStar = refstars->data[match->ref];
+	pmAstromObj *rawStar = rawstars->data[match->raw];
+	
+	matchInfo->match = match; // reference to the match of interest
+	matchInfo->radius = hypot (refStar->FP->x - rawStar->FP->x, refStar->FP->y - rawStar->FP->y);
+
+	psArrayAdd (refSet, 8, matchInfo); // matchInfo->match is just a reference
+	psFree (matchInfo);
+    }
+
+    // we now have a set of matches for each refstar and their distances; create a new set
+    // keeping only the closest entry for each match
+
+    psArray *unique = psArrayAllocEmpty (PS_MAX(16, matches->n / 2));
+    for (int i = 0; i < refstars->n; i++) {
+
+	psArray *refSet = refstarMatches->data[i];
+	if (!refSet) continue;
+	if (refSet->n == 0) continue; // not certain how this can happen...
+
+	if (refSet->n == 1) {
+	    pmAstromMatchInfo *matchInfo = refSet->data[0];
+	    psArrayAdd (unique, 32, matchInfo->match);
+	    continue;
+	}
+
+	pmAstromMatchInfo *matchInfo = refSet->data[0];
+	float minRadius = matchInfo->radius;
+	pmAstromMatch *minMatch = matchInfo->match;
+	for (int j = 1; j < refSet->n; j++) {
+	    pmAstromMatchInfo *matchInfo = refSet->data[j];
+	    if (minRadius < matchInfo->radius) continue;
+	    minMatch = matchInfo->match;
+	    minRadius = matchInfo->radius;
+	}
+	
+	psArrayAdd (unique, 32, minMatch); // minMatch is just a reference to a match on matches,
+    }
+
+    psLogMsg ("psModules.astrom", 3, "generate unique matches to reference stars: %ld matches -> %ld matches\n", matches->n, unique->n);
+    psFree (refstarMatches);
+
+    return unique;
+}