Changeset 10880 for trunk/psastro/src/psastroMosaicOneChip.c

Timestamp:

Jan 2, 2007, 3:26:52 PM (20 years ago)

Author:

eugene

Message:

fixed mosaic analysis; cleaned up test code (some in traces, some removed), general output cleanups

File:

• trunk/psastro/src/psastroMosaicOneChip.c (modified) (5 diffs)

trunk/psastro/src/psastroMosaicOneChip.c

@@ -1 +1 @@
 # include "psastro.h"
 
-bool psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, int iteration) {
+# define REQUIRED_RECIPE_VALUE(VALUE, NAME, TYPE, MESSAGE)\
+  VALUE = psMetadataLookup##TYPE (&status, recipe, NAME); \
+  if (!status) { \
+   psError(PSASTRO_ERR_CONFIG, false, MESSAGE); \
+   return false; } 
+
+bool psastroMosaicOneChip (pmChip *chip, pmReadout *readout, psMetadata *recipe, psMetadata *updates, bool nonlinear) {
 
     bool status;
@@ -20 +26 @@
     if (match == NULL) return false;
 
-    if ((iteration < 1) || (iteration > 2)) {
-        psError(PSASTRO_ERR_UNKNOWN, false, "invalid iteration number %d\n", iteration);
-        return false;
+    // correct radius to FP units (physical pixel scale in microns per pixel)
+    REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32, "Failed to lookup pixel scale"); 
+
+    // allowed limits for valid solutions
+    REQUIRED_RECIPE_VALUE (float maxError, "PSASTRO.MOSAIC.MAX.ERROR", F32, "failed to find single-chip max allowed error\n");
+    REQUIRED_RECIPE_VALUE (int minNstar, "PSASTRO.MOSAIC.MIN.NSTAR", S32, "failed to find single-chip min allowed stars\n");
+
+    // set the order of the per-chip fit (higher order only if nonlinear == true)
+    int order = 1;
+    if (nonlinear) {
+        // select the desired chip order
+        REQUIRED_RECIPE_VALUE (order, "PSASTRO.MOSAIC.CHIP.ORDER", S32, "failed to find mosaic chip-level fit order\n");
+
+        // modify the order to correspond to the actual number of matched stars:
+        if ((match->n < 11) && (order >= 3)) order = 2;
+        if ((match->n <  7) && (order >= 2)) order = 1;
+        if ((match->n <  4) && (order >= 1)) order = 0;
+        if (order < 1) {
+            psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n); 
+            return false; 
+        } 
     }
 
-    // correct to FP units (physical pixel scale in microns per pixel)
-    double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
-    if (!status) { 
-        psError(PS_ERR_IO, false, "Failed to lookup pixel scale"); 
-        return false; 
-    } 
-
-    // set the order of the per-chip fit (iteration 1: 1, iteration 2: from recipe)
-    int order = 1;
-    if (iteration == 2) {
-        order = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.CHIP.ORDER");
-        if (!status) {
-            psError(PSASTRO_ERR_UNKNOWN, false, "failed to find mosaic chip-level fit order\n");
-            return false;
-        }
-    }
+    // create output toFPA; set masks appropriate to the Elixir DVO astrometry format
     psFree (chip->toFPA);
     chip->toFPA = psPlaneTransformAlloc (order, order);
@@ -64 +73 @@
     }
 
-    // allowed limits for valid solutions
-    // XXX have these values depend on the iteration?
-    float maxError = psMetadataLookupF32 (&status, recipe, "PSASTRO.MOSAIC.MAX.ERROR");
-    int minNstar = psMetadataLookupS32 (&status, recipe, "PSASTRO.MOSAIC.MIN.NSTAR");
+    // toSky converts from FPA & TPA units (microns) to sky units (radians)
+    pmFPA *fpa = chip->parent;
+    float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
 
-    // astError is the average 1D scatter in pixels ('results' are in FPA units = microns)
-    // XXX currently using supplied pixel scale: use toFPA->x->coeff[1][0], etc?
-    float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) / pixelScale;
+    // 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;
+
+    // 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;
     int astNstar = results->yStats->clippedNvalues;
 
@@ -77 +87 @@
 
     // XXX should these result in errors or be handled another way?
-    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f, Nstars: %d", astError, astNstar);
+    psLogMsg ("psastro", PS_LOG_INFO, "astrometry solution: error: %f arcsec, Nstars: %d", astError, astNstar);
     if (astError > maxError) {
         psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
@@ -88 +98 @@
 
     // DVO expects NASTRO = 0 if we fail to find a solution
+    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);
     if (validSolution) {
-        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "", astNstar);
-        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", astError/sqrt(astNstar));
+        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
+        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
     } else {
-        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "", 0);
-        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "", 0.0);
+        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0);
+        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", 0);
     }
-    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "", astError);
     psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars
+
+    // determine fromFPA transformation and apply new transformation to raw & ref stars
+    psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
 
     psFree (fitStats);