Index: branches/eam_branches/ipp-dev-20210817/psModules/src/astrom/pmAstrometryObjects.c
===================================================================
--- branches/eam_branches/ipp-dev-20210817/psModules/src/astrom/pmAstrometryObjects.c	(revision 41813)
+++ branches/eam_branches/ipp-dev-20210817/psModules/src/astrom/pmAstrometryObjects.c	(revision 41814)
@@ -340,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;
@@ -345,4 +352,18 @@
         pmAstromObj *rawStar = raw->data[pair->raw];
         if (!isfinite(rawStar->dMag)) continue;
+
+	bool isValid = true;
+	isValid = isValid & isfinite (xFit->data.F32[i]);
+	isValid = isValid & isfinite (yFit->data.F32[i]);
+	isValid = isValid & isfinite (xRes->data.F32[i]);
+	isValid = isValid & isfinite (yRes->data.F32[i]);
+	if (isValid) {
+	  psVectorAppend (xPosValid, xFit->data.F32[i]);
+	  psVectorAppend (yPosValid, yFit->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;
 
@@ -377,10 +398,10 @@
       results->dYrange = pmAstromVectorRange (yResGood, 0.1, 0.9, results->yStats->clippedStdev);
 
-      results->dXstdev = pmAstrom2DSystematics (xFit, yFit, xRes);
-      results->dYstdev = pmAstrom2DSystematics (xFit, yFit, yRes);
-    }
-
-    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);
@@ -390,4 +411,8 @@
     psFree (xResGood);
     psFree (yResGood);
+    psFree (xPosValid);
+    psFree (yPosValid);
+    psFree (xResValid);
+    psFree (yResValid);
 
     psFree (x);
Index: branches/eam_branches/ipp-dev-20210817/psModules/src/astrom/pmAstrometryObjects.h
===================================================================
--- branches/eam_branches/ipp-dev-20210817/psModules/src/astrom/pmAstrometryObjects.h	(revision 41813)
+++ branches/eam_branches/ipp-dev-20210817/psModules/src/astrom/pmAstrometryObjects.h	(revision 41814)
@@ -102,4 +102,6 @@
     double  dXrange;			///< 10% - 90% range X residuals (unmasked, high S/N)
     double  dYrange;			///< 10% - 90% range Y residuals (unmasked, high S/N)
+    double  dXstdev;			///< stdev of median residual in NxN bins
+    double  dYstdev;			///< stdev of median residual in NxN bins
 }
 pmAstromFitResults;
Index: branches/eam_branches/ipp-dev-20210817/psastro/src/psastroMosaicOneChip.c
===================================================================
--- branches/eam_branches/ipp-dev-20210817/psastro/src/psastroMosaicOneChip.c	(revision 41813)
+++ branches/eam_branches/ipp-dev-20210817/psastro/src/psastroMosaicOneChip.c	(revision 41814)
@@ -134,4 +134,8 @@
     int astNstar = results->yStats->clippedNvalues;
 
+    // astStdev is the average 1D stdev of median residuals in arcsec ('results' are in FPA units = microns)
+    // the median residuals are calculated in a grid of N x N bins 
+    float astStdev = 0.5*(results->xStdev + results->yStdev) * plateScale;
+
     // if we clip away too many stars, the order may be invalid
     if (order == 3) { minNstar = PS_MAX (15, minNstar); }
@@ -142,5 +146,6 @@
 
     // XXX should these result in errors or be handled another way?
-    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d", astError, astNstar);
+    // psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d", astError, astNstar);
+    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d, stdev: %f arcsec", astError, astNstar, astStdev);
     if ((maxError > 0) && (astError > maxError)) {
         psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
@@ -155,4 +160,5 @@
     psMetadataAddF32 (updates, PS_LIST_TAIL, "PERROR",   PS_META_REPLACE, "astrometry error (pixels)", pixError);
     psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "astrometry error (arcsec)", astError);
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERSTD",   PS_META_REPLACE, "astrometry stdev (arcsec)", astStdev);
     if (validSolution) {
         psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
Index: branches/eam_branches/ipp-dev-20210817/psastro/src/psastroOneChipFit.c
===================================================================
--- branches/eam_branches/ipp-dev-20210817/psastro/src/psastroOneChipFit.c	(revision 41813)
+++ branches/eam_branches/ipp-dev-20210817/psastro/src/psastroOneChipFit.c	(revision 41814)
@@ -192,4 +192,8 @@
     float astError = 0.5*(rawXstdev + rawYstdev) * plateScale;
 
+    // astStdev is the average 1D stdev of median residuals in arcsec ('results' are in FPA units = microns)
+    // the median residuals are calculated in a grid of N x N bins 
+    float astStdev = 0.5*(results->xStdev + results->yStdev) * plateScale;
+
     // x and y are forced to use the same subset of values:
     int astNstar = results->yStats->clippedNvalues;
@@ -198,5 +202,5 @@
 
     // XXX should these result in errors or be handled another way?
-    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d", astError, astNstar);
+    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d, stdev: %f arcsec", astError, astNstar, astStdev);
     if (astError > maxError) {
         psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
@@ -211,4 +215,5 @@
     psMetadataAddF32 (updates, PS_LIST_TAIL, "PERROR",   PS_META_REPLACE, "astrometry error (pixels)", pixError);
     psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "astrometry error (arcsec)", astError);
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERSTD",   PS_META_REPLACE, "astrometry stdev (arcsec)", astStdev);
     if (validSolution) {
         psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));