Changeset 15891 for trunk/psastro/src/psastroModelAnalysis.c

Timestamp:

Dec 22, 2007, 7:42:46 AM (18 years ago)

Author:

eugene

Message:

finished psastroModel for boresite, created gpcModel for rough theoretical model, added useModel and fixChips options

File:

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

trunk/psastro/src/psastroModelAnalysis.c

@@ -12 +12 @@
     }
 
+    pmFPAfile *output = psMetadataLookupPtr (&status, config->files, "PSASTRO.OUT.MODEL");
+    if (!status) psAbort ("Can't find output pmFPAfile PSASTRO.OUT.MODEL");
+
     // physical pixel scale in microns per pixel
     char *refChip = psMetadataLookupStr (&status, recipe, "PSASTRO.MODEL.REF.CHIP");
@@ -19 +22 @@
     } 
 
-    /* model analysis steps:
-
-    1) determine POS_ZERO:
-    * choose a reference chip
-    * what is posangle of chip
-    * posangle = FPA.POSANGLE - POS_ZERO
-    * POS_ZERO = FPA.POSANGLE - posangle
-
-    ** average over input images
-
-    2) determine boresite model
-
-    * find CRPIX1,2 for a series of images
-    * fit ellipse 
-
-    */
-
-    /***** find the POSANGLE offset *****/
+    /* model analysis:
+     *
+     * determine POS_ZERO via comparison of measured and reported posangles
+     * POS_ZERO = FPA.POSANGLE - posangle
+     *
+     * determine boresite model:
+     * X = Xo + R_X cos(FPA.POSANGLE - T_0) cos(P_0) + R_Y sin(FPA.POSANGLE - T_0) sin(P_0) 
+     * Y = Yo + R_Y sin(FPA.POSANGLE - T_0) cos(P_0) - R_X cos(FPA.POSANGLE - T_0) sin(P_0) 
+     * position of reported boresite in reference chip pixels
+     * Xo, Yo : true coordinate of boresite (rotator center) in reference chip pixels
+     * R_X, R_Y : amplitude of boresite offset
+     * T_0 : reference angle for rotator
+     * P_0 : orientation of boresite ellipse
+     */
 
     // select the input pmFPAfile pointers
@@ -44 +43 @@
     psArray *files = psListToArray (item->data.list);
 
-    // temp data vectors
+    // data storage vectors for measurements
     psVector *posZero  = psVectorAlloc (files->n, PS_TYPE_F32);
     psVector *Po       = psVectorAlloc (files->n, PS_TYPE_F32);
-    psVector *Lo       = psVectorAlloc (files->n, PS_TYPE_F32);
-    psVector *Mo       = psVectorAlloc (files->n, PS_TYPE_F32);
+    psVector *Xo       = psVectorAlloc (files->n, PS_TYPE_F32);
+    psVector *Yo       = psVectorAlloc (files->n, PS_TYPE_F32);
 
+    // counter for accepted measured values
     int n = 0;
 
-    // convert the headers for the input file into fpa astrometry terms
+    // Re-select the output fpa.  The output fpa needs to have valid astrometry for the
+    // refChip.  output->fpa is a copy of the pointer to one of the input->fpa, but the choice
+    // is arbitrary.  select a new one that has an existing ref chip
+    psFree (output->fpa);
+    output->fpa = NULL;
+
+    // extract the relevant measured and reported values from the reference chip
     for (int i = 0; i < files->n; i++) {
         psMetadataItem *file = files->data[i];
         pmFPAfile *input = file->data.V;
 
+        // reported rotator position angle (this should perhaps be ROT, not POS)
         double POSANGLE = psMetadataLookupF64 (&status, input->fpa->concepts, "FPA.POSANGLE"); 
         if (!status) psAbort ("missing FPA.POSANGLE");
 
-        // get chip from name
+        // get reference chip from name
         pmChip *chip = pmConceptsChipFromName (input->fpa, refChip);
         if (!chip) psAbort ("invalid chip name for reference");
 
+        fprintf (stderr, "input %d : %zx : %zx : %zx\n", i, (size_t) input, (size_t) chip, (size_t) chip->toFPA);
         if (!chip->toFPA) continue;
+
+        if (!output->fpa) {
+            // this one matches
+            output->fpa = psMemIncrRefCounter(input->fpa);
+        }
 
         // we have two measurements of the posangle (may be parity flipped to different quadrants)
@@ -72 +85 @@
         while (posZero->data.F32[n] <   0.0) posZero->data.F32[n] += 360.0;
 
-        Po->data.F32[n] = POSANGLE * PM_RAD_DEG;
-        Lo->data.F32[n] = chip->toFPA->x->coeff[0][0];
-        Mo->data.F32[n] = chip->toFPA->y->coeff[0][0];
-        fprintf (stderr, "%d : %f %f : %f = %f - %f\n", i, Lo->data.F32[n], Mo->data.F32[n], posZero->data.F32[n], POSANGLE, posangle);
+        Po->data.F32[n] = POSANGLE * PM_RAD_DEG; // reported position angle
+        Xo->data.F32[n] = chip->fromFPA->x->coeff[0][0]; // reported boresite x position in ref chip coordinates
+        Yo->data.F32[n] = chip->fromFPA->y->coeff[0][0]; // reported boresite y position in ref chip coordinates
+        fprintf (stderr, "%d : %f %f : %f = %f - %f\n", i, Xo->data.F32[n], Yo->data.F32[n], posZero->data.F32[n], POSANGLE, posangle);
         n ++;
     }
       
-    Lo->n = n;
-    Mo->n = n;
+    Xo->n = n;
+    Yo->n = n;
     Po->n = n;
     posZero->n = n;
@@ -88 +101 @@
 
     fprintf (stderr, "pos zero %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
+    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.POS_ZERO", PS_META_REPLACE, "offset between obs and meas posangle", stats->sampleMedian); 
 
-    psastroModelFitBoresite (Lo, Mo, Po);
+    psVector *params = psastroModelFitBoresite (Xo, Yo, Po);
+    if (params->n != 6) psAbort ("error");
+
+    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.X0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_X0]); 
+    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.Y0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_Y0]); 
+    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.RX", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_RX]); 
+    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.RY", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_RY]); 
+    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.T0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_T0]); 
+    psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.P0", PS_META_REPLACE, "boresite parameter", params->data.F32[PAR_P0]); 
+    psMetadataAddStr (output->fpa->concepts, PS_LIST_TAIL, "FPA.REF.CHIP", PS_META_REPLACE, "boresite parameter", refChip);
 
     return true;