Changeset 19513

Timestamp:

Sep 11, 2008, 3:08:49 PM (18 years ago)

Author:

eugene

Message:

adding a test of the original guess

Location:

trunk/psastro/src

Files:

• psastroAnalysis.c (modified) (1 diff)
• psastroAstromGuess.c (modified) (7 diffs)

trunk/psastro/src/psastroAnalysis.c

@@ -81 +81 @@
     // psastroZeroPoint (config);
 
+    psastroAstromGuessCheck (config);
+
     // XXX how do we specify stack astrometry?
     // psastroStackAstrom (config, refs);

trunk/psastro/src/psastroAstromGuess.c

@@ -1 +1 @@
 # include "psastroInternal.h"
+# define DEBUG 0
 
 // this function loads the header WCS astrometry terms into the fpa terms and applies the
@@ -52 +53 @@
     } 
 
+    psVector *cornerL = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerM = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerP = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerQ = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerR = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerD = psVectorAllocEmpty (100, PS_TYPE_F32);
+
     pmFPA *fpa = input->fpa;
+
+    if (DEBUG) psastroDumpCorners ("corners.up.guess1.dat", "corners.dn.guess1.dat", fpa);
 
     // load mosaic-level astrometry?
@@ -77 +87 @@
         }
 
+        // report and save the current best guess for the chip 0,0 pixel coordinates
+        { 
+            psPlane ptCH, ptFP, ptTP;
+            psSphere ptSky;
+
+            ptCH.x = 0;
+            ptCH.y = 0;
+            psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
+            psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
+            psDeproject (&ptSky, &ptTP, fpa->toSky);
+            psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
+
+            psVectorAppend (cornerL, ptFP.x);
+            psVectorAppend (cornerM, ptFP.y);
+            psVectorAppend (cornerP, ptTP.x);
+            psVectorAppend (cornerQ, ptTP.y);
+            psVectorAppend (cornerR, ptSky.r);
+            psVectorAppend (cornerD, ptSky.d);
+        }
+
         // apply the new WCS guess data to all of the data in the readouts
         while ((cell = pmFPAviewNextCell (view, fpa, 1)) != NULL) {
@@ -85 +115 @@
             while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
                 if (! readout->data_exists) { continue; }
-
-                // report the current best guess for the cell 0,0 pixel coordinate
-                { 
-                  psPlane ptCH, ptFP, ptTP;
-                  psSphere ptSky;
-
-                  ptCH.x = 0;
-                  ptCH.y = 0;
-                  psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
-                  psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
-                  psDeproject (&ptSky, &ptTP, fpa->toSky);
-                  psLogMsg ("psastro", 2, "0,0 pix for chip,cell %d,%d = %f,%f\n", view->chip, view->cell, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
-                }
 
                 psArray *rawstars = psMetadataLookupPtr (NULL, readout->analysis, "PSASTRO.RAWSTARS");
@@ -132 +149 @@
         }
     }
+
+    if (DEBUG) psastroDumpCorners ("corners.up.guess2.dat", "corners.dn.guess2.dat", fpa);
 
     // how many total sources are available to us?
@@ -149 +168 @@
     psMetadataAddF32 (recipe, PS_LIST_TAIL, "DEC_MIN", PS_META_REPLACE, "", DECmin);
     psMetadataAddF32 (recipe, PS_LIST_TAIL, "DEC_MAX", PS_META_REPLACE, "", DECmax);
+
+    psMetadataAddVector (input->fpa->analysis, PS_LIST_TAIL, "CORNER.L", PS_META_REPLACE, "corner pixel", cornerL);
+    psMetadataAddVector (input->fpa->analysis, PS_LIST_TAIL, "CORNER.M", PS_META_REPLACE, "corner pixel", cornerM);
+    psMetadataAddVector (input->fpa->analysis, PS_LIST_TAIL, "CORNER.P", PS_META_REPLACE, "corner pixel", cornerP);
+    psMetadataAddVector (input->fpa->analysis, PS_LIST_TAIL, "CORNER.Q", PS_META_REPLACE, "corner pixel", cornerQ);
+    psMetadataAddVector (input->fpa->analysis, PS_LIST_TAIL, "CORNER.R", PS_META_REPLACE, "corner pixel", cornerR);
+    psMetadataAddVector (input->fpa->analysis, PS_LIST_TAIL, "CORNER.D", PS_META_REPLACE, "corner pixel", cornerD);
+
+    psFree (cornerL);
+    psFree (cornerM);
+    psFree (cornerP);
+    psFree (cornerQ);
+    psFree (cornerR);
+    psFree (cornerD);
 
     psFree (view);
@@ -201 +234 @@
     return true;
 }
+
+// we made a guess at the beginning; how does the guess compare with the result?
+bool psastroAstromGuessCheck (pmConfig *config) {
+
+    bool status;
+
+    // select the input data sources
+    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
+    if (!input) {
+        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
+        return false;
+    }
+
+    psVector *cornerLn = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerMn = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerPn = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerQn = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerRn = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerDn = psVectorAllocEmpty (100, PS_TYPE_F32);
+
+    pmFPA *fpa = input->fpa;
+
+    if (DEBUG) psastroDumpCorners ("corners.up.guess3.dat", "corners.dn.guess3.dat", fpa);
+
+    pmChip *chip = NULL;
+    pmFPAview *view = pmFPAviewAlloc (0);
+
+    while ((chip = pmFPAviewNextChip (view, fpa, 1)) != NULL) {
+        if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; }
+
+        psPlane ptCH, ptFP, ptTP;
+        psSphere ptSky;
+
+        ptCH.x = 0;
+        ptCH.y = 0;
+        psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
+        psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
+        psDeproject (&ptSky, &ptTP, fpa->toSky);
+        psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
+
+        // new corner locations based on the calibrated astrometry
+        psVectorAppend (cornerLn, ptFP.x);
+        psVectorAppend (cornerMn, ptFP.y);
+        psVectorAppend (cornerPn, ptTP.x);
+        psVectorAppend (cornerQn, ptTP.y);
+        psVectorAppend (cornerRn, ptSky.r);
+        psVectorAppend (cornerDn, ptSky.d);
+    }
+
+    psVector *cornerLo = psMetadataLookupPtr (&status, input->fpa->analysis, "CORNER.L");
+    psVector *cornerMo = psMetadataLookupPtr (&status, input->fpa->analysis, "CORNER.M");
+    psVector *cornerPo = psMetadataLookupPtr (&status, input->fpa->analysis, "CORNER.P");
+    psVector *cornerQo = psMetadataLookupPtr (&status, input->fpa->analysis, "CORNER.Q");
+    psVector *cornerRo = psMetadataLookupPtr (&status, input->fpa->analysis, "CORNER.R");
+    psVector *cornerDo = psMetadataLookupPtr (&status, input->fpa->analysis, "CORNER.D");
+
+    // compare the old R,D values projected to the same tangent plane as the new R,D values:
+
+    psVector *cornerPs = psVectorAllocEmpty (100, PS_TYPE_F32);
+    psVector *cornerQs = psVectorAllocEmpty (100, PS_TYPE_F32);
+
+    for (int i = 0; i < cornerRo->n; i++) {
+        
+        psPlane ptTP;
+        psSphere ptSky;
+
+        ptSky.r = cornerRo->data.F32[i];
+        ptSky.d = cornerDo->data.F32[i];
+
+        psProject (&ptTP, &ptSky, fpa->toSky);
+        psVectorAppend (cornerPs, ptTP.x);
+        psVectorAppend (cornerQs, ptTP.y);
+    }
+
+    psPlaneTransform *map = psPlaneTransformAlloc (1, 1);
+    map->x->coeffMask[1][1] = PS_POLY_MASK_SET;
+    map->y->coeffMask[1][1] = PS_POLY_MASK_SET;
+    
+    psVectorFitPolynomial2D (map->x, NULL, 0, cornerPn, NULL, cornerPs, cornerQs);
+    psVectorFitPolynomial2D (map->y, NULL, 0, cornerQn, NULL, cornerPs, cornerQs);
+    
+    // apply the linear fit...
+    psVector *cornerPf = psPolynomial2DEvalVector (map->x, cornerPs, cornerQs);
+    psVector *cornerQf = psPolynomial2DEvalVector (map->y, cornerPs, cornerQs);
+
+    // ...and calculate the residual between Pn,Qn and Pf,Qf
+    psVector *cornerPd = (psVector *) psBinaryOp (NULL, cornerPn, "-", cornerPf);
+    psVector *cornerQd = (psVector *) psBinaryOp (NULL, cornerQn, "-", cornerQf);
+
+    psStats *statsP = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+    psStats *statsQ = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+
+    psVectorStats (statsP, cornerPd, NULL, NULL, 0);
+    psVectorStats (statsQ, cornerQd, NULL, NULL, 0);
+
+    float angle = atan2 (map->y->coeff[1][0], map->x->coeff[1][0]);
+    float scale = hypot (map->y->coeff[1][0], map->x->coeff[1][0]);
+    
+    psLogMsg ("psastro", 3, "boresite offset  : %f,%f\n", map->x->coeff[0][0], map->y->coeff[0][0]);
+    psLogMsg ("psastro", 3, "boresite angle   : %f, scale: %f", angle*PS_DEG_RAD, scale);
+    psLogMsg ("psastro", 3, "boresite scatter : %f,%f\n", statsP->sampleStdev, statsQ->sampleStdev);
+
+    // write the elapsed time here; this will be updated in psastroMosaicAstrometry, if called
+    psMetadata *header = psMetadataLookupMetadata (&status, input->fpa->analysis, "PSASTRO.HEADER");
+    if (!header) {
+        header = psMetadataAlloc();
+        psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", header);
+        psFree (header);  // drop this reference
+    }
+
+    psMetadataAddF32 (header, PS_LIST_TAIL, "AST_R0", PS_META_REPLACE, "boresite offset in RA (TP units)", map->x->coeff[0][0]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "AST_D0", PS_META_REPLACE, "boresite offset in DEC (TP units)", map->y->coeff[0][0]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "AST_T0", PS_META_REPLACE, "boresite angle (degrees)", angle*PS_DEG_RAD);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "AST_S0", PS_META_REPLACE, "boresite scale correction", scale);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "AST_RS", PS_META_REPLACE, "boresite scatter in RA (TP units)", statsP->sampleStdev);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "AST_DS", PS_META_REPLACE, "boresite scatter in DEC (TP units)", statsQ->sampleStdev);
+
+    if (0) {
+        FILE *f = fopen ("corners.dat", "w");
+        for (int i = 0; i < cornerRo->n; i++) {
+            fprintf (f, "%10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f  |  %10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f\n",
+                     cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i], 
+                     cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]);
+        }
+        fclose (f);
+    }
+
+    psFree (cornerPf);
+    psFree (cornerQf);
+    psFree (cornerPd);
+    psFree (cornerQd);
+
+    psFree (statsP);
+    psFree (statsQ);
+
+    psFree (cornerLn);
+    psFree (cornerMn);
+    psFree (cornerPn);
+    psFree (cornerQn);
+    psFree (cornerRn);
+    psFree (cornerDn);
+    psFree (cornerPs);
+    psFree (cornerQs);
+    psFree (map);
+    psFree (view);
+    
+
+    return true;
+}