Changeset 26227

Timestamp:

Nov 20, 2009, 3:24:18 PM (17 years ago)

Author:

eugene

Message:

add additional statistical measures to astrometry calibration

File:

• branches/eam_branches/20091113/psastro/src/psastroOneChipFit.c (modified) (3 diffs)

branches/eam_branches/20091113/psastro/src/psastroOneChipFit.c

@@ -135 +135 @@
         // toSky converts from FPA & TPA units (microns) to sky units (radians)
         float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
-        // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
-        float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
+
+        float rawXstdev = psStatsGetValue (results->xStats, psStatsStdevOption(results->xStats->options));
+        float rawYstdev = psStatsGetValue (results->yStats, psStatsStdevOption(results->yStats->options));
+
+        float astError = 0.5*(rawXstdev + rawYstdev) * plateScale;
         int astNstar = results->yStats->clippedNvalues;
         psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar);
@@ -150 +153 @@
     float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
 
+    float rawXstdev = psStatsGetValue (results->xStats, psStatsStdevOption(results->xStats->options));
+    float rawYstdev = psStatsGetValue (results->yStats, psStatsStdevOption(results->yStats->options));
+
     // pixError is the average 1D scatter in pixels ('results' are in FPA units = microns)
-    // float pixError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale;
-    float pixError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) / pixelScale;
+    float pixError = 0.5*(rawXstdev + rawYstdev) / pixelScale;
 
     // astError is the average 1D scatter in arcsec ('results' are in FPA units = microns)
-    // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
-    float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
+    float astError = 0.5*(rawXstdev + rawYstdev) * plateScale;
+
+    // x and y are forced to use the same subset of values:
     int astNstar = results->yStats->clippedNvalues;
 
@@ -184 +190 @@
     psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars
 
+    // additional error measurements:
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CDX",   PS_META_REPLACE, "chip astrometry X stdev (arcsec)", rawXstdev * plateScale);
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CSX",   PS_META_REPLACE, "chip astrometry X systematic err (arcsec)", results->dXsys * plateScale);
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CRX",   PS_META_REPLACE, "chip astrometry X 10-90 percentile (arcsec)", results->dXrange * plateScale);
+
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CDY",   PS_META_REPLACE, "chip astrometry Y stdev (arcsec)", rawYstdev * plateScale);
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CSY",   PS_META_REPLACE, "chip astrometry Y systematic err (arcsec)", results->dYsys * plateScale);
+    psMetadataAddF32 (updates, PS_LIST_TAIL, "AST_CRY",   PS_META_REPLACE, "chip astrometry Y 10-90 percentile (arcsec)", results->dYrange * plateScale);
+
     // XXX check if we correctly applied the new transformation:
     if (psTraceGetLevel("psastro.dump") > 0) {