Changeset 20033 for branches/cnb_branch_20080830/psastro/src/psastroModelAnalysis.c

Timestamp:

Oct 9, 2008, 2:39:01 PM (18 years ago)

Author:

beaumont

Message:

merged mainline into branch. resolved conflicts. added plots.

Location:

branches/cnb_branch_20080830/psastro/src

Files:

• . (modified) (1 prop)
• psastroModelAnalysis.c (modified) (9 diffs)

branches/cnb_branch_20080830/psastro/src

@@ -15 +15 @@
 psastroModel
 gpcModel
+psastroModelFit

@@ -15 +15 @@
 psastroModel
 gpcModel
+psastroModelFit

branches/cnb_branch_20080830/psastro/src/psastroModelAnalysis.c

@@ -1 +1 @@
 # include "psastroStandAlone.h"
+# define NONLIN_TOL 0.001
 
 bool psastroModelAnalysis (pmConfig *config) {
@@ -12 +13 @@
     }
 
-    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
+    // reference chip for boresite parameters
     char *refChip = psMetadataLookupStr (&status, recipe, "PSASTRO.MODEL.REF.CHIP");
     if (!refChip) {
@@ -21 +19 @@
         return false; 
     } 
+
+    pmFPAfile *output = psMetadataLookupPtr (&status, config->files, "PSASTRO.OUT.MODEL");
+    if (!status) psAbort ("Can't find output pmFPAfile PSASTRO.OUT.MODEL");
+
+    // measure the boresite position from a rotation sequence?
+    bool fitBoresite = psMetadataLookupBool (&status, recipe, "PSASTRO.MODEL.FIT.BORESITE");
+    if (!fitBoresite) {
+        psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.X0", PS_META_REPLACE, "boresite parameter", 0.0); 
+        psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.Y0", PS_META_REPLACE, "boresite parameter", 0.0); 
+        psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.RX", PS_META_REPLACE, "boresite parameter", 0.0); 
+        psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.RY", PS_META_REPLACE, "boresite parameter", 0.0); 
+        psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.T0", PS_META_REPLACE, "boresite parameter", 0.0); 
+        psMetadataAddF32 (output->fpa->concepts, PS_LIST_TAIL, "FPA.BORE.P0", PS_META_REPLACE, "boresite parameter", 0.0); 
+        psMetadataAddStr (output->fpa->concepts, PS_LIST_TAIL, "FPA.REF.CHIP", PS_META_REPLACE, "boresite parameter", refChip);
+        return true; 
+    } 
+
+    char *outroot = psMetadataLookupStr (&status, config->arguments, "OUTPUT");
+    if (!status || !outroot) psAbort ("Can't find outroot on config->arguments");
 
     /* model analysis:
@@ -35 +52 @@
      * T_0 : reference angle for rotator
      * P_0 : orientation of boresite ellipse
+     *
      */
 
@@ -58 +76 @@
     output->fpa = NULL;
 
+    char filename[256];
+    snprintf (filename, 256, "%s.bore", outroot);
+    FILE *outfile = fopen (filename, "w");
+    if (!outfile) {
+        psError(PSASTRO_ERR_ARGUMENTS, true, "cannot open %s for output", filename);
+        return false;
+    }
+
+    float posBoundary = 0.0;
+
     // extract the relevant measured and reported values from the reference chip
     for (int i = 0; i < files->n; i++) {
@@ -63 +91 @@
         pmFPAfile *input = file->data.V;
 
-        // reported rotator position angle (this should perhaps be ROT, not POS)
+        // reported rotator position angle
         double POSANGLE = psMetadataLookupF64 (&status, input->fpa->concepts, "FPA.POSANGLE"); 
         if (!status) psAbort ("missing FPA.POSANGLE");
@@ -71 +99 @@
         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;
 
@@ -80 +107 @@
 
         // we have two measurements of the posangle (may be parity flipped to different quadrants)
-        float posangle = PM_DEG_RAD * atan2 (chip->toFPA->y->coeff[1][0], chip->toFPA->x->coeff[1][0]);
-        posZero->data.F32[n] = POSANGLE - posangle;
-        while (posZero->data.F32[n] > 360.0) posZero->data.F32[n] -= 360.0;
-        while (posZero->data.F32[n] <   0.0) posZero->data.F32[n] += 360.0;
+        // atan2 returns values in the range 0-2pi
+        // all posZero values should be clustered in some region, but we need to flip over the 0,360 boundary correctly.
+        // push all to one side or the other
+        float chipAngle = PM_DEG_RAD * atan2 (chip->toFPA->y->coeff[1][0], chip->toFPA->x->coeff[1][0]);
+        float fpaAngle = PM_DEG_RAD * atan2 (input->fpa->toTPA->y->coeff[1][0], input->fpa->toTPA->x->coeff[1][0]);
+
+        posZero->data.F32[n] = POSANGLE - chipAngle - fpaAngle;
+        if (n == 0) {
+            posBoundary = posZero->data.F32[n] + 180.0;
+        } else {
+            while (posZero->data.F32[n] > posBoundary) posZero->data.F32[n] -= 360.0;
+            while (posZero->data.F32[n] < posBoundary - 360.0) posZero->data.F32[n] += 360.0;
+        }
 
         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);
+        float xc = chip->fromFPA->x->coeff[0][0]; // reported boresite x position in ref chip coordinates
+        float yc = chip->fromFPA->y->coeff[0][0]; // reported boresite y position in ref chip coordinates
+        // XXX this can also be derived from toFPA via GetCenter....
+        
+        psPlane *PT = psPlaneTransformGetCenter (chip->toFPA, NONLIN_TOL);
+        Xo->data.F32[n] = PT->x; // reported boresite x position in ref chip coordinates
+        Yo->data.F32[n] = PT->y; // reported boresite y position in ref chip coordinates
+        psFree (PT);
+
+        fprintf (outfile, "%d : %f %f : %f = %f - %f - %f | %f %f\n", i, Xo->data.F32[n], Yo->data.F32[n], posZero->data.F32[n], POSANGLE, chipAngle, fpaAngle, xc, yc);
         n ++;
     }
@@ -100 +143 @@
     psVectorStats (stats, posZero, NULL, NULL, 0);
 
-    fprintf (stderr, "pos zero %f +/- %f\n", stats->sampleMedian, stats->sampleStdev);
+    fprintf (outfile, "# 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); 
+    fclose (outfile);
 
-    psVector *params = psastroModelFitBoresite (Xo, Yo, Po);
+    psVector *params = psastroModelFitBoresite (Xo, Yo, Po, outroot);
     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]); 

@@ -15 +15 @@
 psastroModel
 gpcModel
+psastroModelFit