Index: trunk/psModules/src/astrom/pmAstrometryObjects.c
===================================================================
--- trunk/psModules/src/astrom/pmAstrometryObjects.c	(revision 41493)
+++ trunk/psModules/src/astrom/pmAstrometryObjects.c	(revision 41892)
@@ -40,4 +40,5 @@
 // XXX this is defined in pmPSFtry.h, which makes no sense
 float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction);
+float pmAstrom2DSystematics (psVector *xPos, psVector *yPos, psVector *value);
 
 #define PM_ASTROMETRYOBJECTS_DEBUG 1
@@ -339,4 +340,11 @@
     psVector *yResGood = psVectorAllocEmpty (match->n, PS_TYPE_F32);
 
+    // we measure the stdev of the median residual in NxN bins.
+    // use only valid (not NAN) measurements
+    psVector *xPosValid = psVectorAllocEmpty (match->n, PS_TYPE_F32);
+    psVector *yPosValid = psVectorAllocEmpty (match->n, PS_TYPE_F32);
+    psVector *xResValid = psVectorAllocEmpty (match->n, PS_TYPE_F32);
+    psVector *yResValid = psVectorAllocEmpty (match->n, PS_TYPE_F32);
+
     for (int i = 0; i < match->n; i++) {
         if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue;
@@ -344,4 +352,18 @@
         pmAstromObj *rawStar = raw->data[pair->raw];
         if (!isfinite(rawStar->dMag)) continue;
+
+	bool isValid = true;
+	isValid = isValid & isfinite (x->data.F32[i]);
+	isValid = isValid & isfinite (y->data.F32[i]);
+	isValid = isValid & isfinite (xRes->data.F32[i]);
+	isValid = isValid & isfinite (yRes->data.F32[i]);
+	if (isValid) {
+	  psVectorAppend (xPosValid, x->data.F32[i]);
+	  psVectorAppend (yPosValid, y->data.F32[i]);
+	  psVectorAppend (xResValid, xRes->data.F32[i]);
+	  psVectorAppend (yResValid, yRes->data.F32[i]);
+	}
+
+	// for the systematic error, use only high S/N stars
         if (rawStar->dMag > 0.02) continue;
 
@@ -369,4 +391,5 @@
       results->dXsys   = results->dYsys   = NAN;
       results->dXrange = results->dYrange = NAN;
+      results->dXstdev = results->dYstdev = NAN;
     } else {
       results->dXsys = psVectorSystematicError (xResGood, xErr, 0.05);
@@ -375,8 +398,11 @@
       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);
+
+      results->dXstdev = pmAstrom2DSystematics (xPosValid, yPosValid, xResValid);
+      results->dYstdev = pmAstrom2DSystematics (xPosValid, yPosValid, yResValid);
+    }
+
+    psTrace ("psModules.astrom", 3, "dXsys: %f, dXrange: %f, dXstdev: %f\n", results->dXsys, results->dXrange, results->dXstdev);
+    psTrace ("psModules.astrom", 3, "dYsys: %f, dYrange: %f, dYstdev: %f\n", results->dYsys, results->dYrange, results->dYstdev);
 
     psFree (xErr);
@@ -386,4 +412,8 @@
     psFree (xResGood);
     psFree (yResGood);
+    psFree (xPosValid);
+    psFree (yPosValid);
+    psFree (xResValid);
+    psFree (yResValid);
 
     psFree (x);
@@ -395,4 +425,97 @@
 
     return (results);
+}
+
+# define VAL_COUNT 9.0
+# define MIN_COUNT 5.0
+
+float pmAstrom2DSystematics (psVector *xPos, psVector *yPos, psVector *value) {
+
+    // pre-filter the values to ensure no NANs, other invalid
+    if (xPos->n < VAL_COUNT*2*2) {
+      return NAN;
+    }
+
+    // generate a grid covering the full range of x,y and measure the median value in each
+    // grid cell.  calculate the stdev of those median values.
+
+    // VAL_COUNT (9) is the (min) goal density, but a single cell may have only MIN_COUNT (7)
+    // we have N points.  require a min of 9 pts per cell (configurable?).  grid is square.
+    // Ncell*Ncell*9 = Npts, Ncell = MIN(sqrt(Npts/9), 5)
+
+    int Ncell = PS_MIN(sqrt((float)(xPos->n / VAL_COUNT)), 2);  // have to at least have a 2x2 grid
+
+    // find the range of x,y values
+    float xMin = 1e9, xMax = -1e9, yMin = 1e9, yMax = 1e9;
+    for (int i = 0; i < xPos->n; i++) {
+	xMin = PS_MIN(xPos->data.F32[i],xMin);
+	xMax = PS_MAX(xPos->data.F32[i],xMax);
+	yMin = PS_MIN(yPos->data.F32[i],yMin);
+	yMax = PS_MAX(yPos->data.F32[i],yMax);
+    }
+
+    float xStep = Ncell / (xMax - xMin);
+    float yStep = Ncell / (yMax - yMin);
+  
+    if (isnan(xStep)) {
+      return NAN;
+    }
+    if (isnan(yStep)) {
+      return NAN;
+    }
+
+    psArray *xBin = psArrayAlloc (Ncell);
+    for (int ix = 0; ix < Ncell; ix++) {
+	psArray *yBin = psArrayAlloc (Ncell); 
+	xBin->data[ix] = yBin;
+	for (int iy = 0; iy < Ncell; iy++) {
+	    yBin->data[iy] = psVectorAllocEmpty(128, PS_TYPE_F32);
+	}
+    }    
+
+    // xValue = ix/xStep + xMin -> ix = (xValue - xMin) * xStep;
+
+    for (int i = 0; i < xPos->n; i++) {
+	int ix = PS_MIN(PS_MAX(0, (xPos->data.F32[i] - xMin) * xStep), Ncell - 1);
+	int iy = PS_MIN(PS_MAX(0, (yPos->data.F32[i] - yMin) * yStep), Ncell - 1);
+    
+	psArray *yBin = xBin->data[ix];
+	psVector *Bin = yBin->data[iy];
+
+	psVectorAppend(Bin, value->data.F32[i]);
+    }
+
+    psVector *sample = psVectorAllocEmpty (Ncell*Ncell, PS_TYPE_F32);
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+    // calculate the median for each vector and save on a vector
+    for (int ix = 0; ix < Ncell; ix++) {
+	psArray *yBin = xBin->data[ix];
+	for (int iy = 0; iy < Ncell; iy++) {
+	    psVector *Bin = yBin->data[iy];
+	    if (Bin->n < MIN_COUNT) { continue; }
+
+	    // psVectorStats resets stats so we can call this repeatedly
+	    psVectorStats (stats, Bin, NULL, NULL, 0); 
+	    if (isfinite(stats->sampleMedian)) {
+		psVectorAppend (sample, stats->sampleMedian);
+	    }
+	}
+    }    
+  
+    stats->options = PS_STAT_SAMPLE_STDEV;
+    psVectorStats (stats, sample, NULL, NULL, 0); 
+    float result = stats->sampleStdev;
+
+    if (!isfinite(stats->sampleStdev)) {
+      fprintf (stderr, "*** bad solution ***\n");
+    }
+
+    // NOTE: the elements of xBin are freed automatically, which extends down to the vectors
+    psFree (xBin); 
+    psFree (stats);
+    psFree (sample);
+  
+    return result;
 }