Changeset 5510

Timestamp:

Nov 13, 2005, 1:16:22 PM (21 years ago)

Author:

eugene

Message:

development work

Location:

trunk/psastro

Files:

• doc (added)
• doc/notes.txt (added)
• src/pmAstrom.c (modified) (4 diffs)
• src/pmAstrom.h (modified) (2 diffs)
• src/pmAstromGrid.c (modified) (4 diffs)
• src/psastro.c (modified) (1 diff)
• src/psastroBuildFPA.c (modified) (1 diff)
• src/psastroUtils.c (modified) (2 diffs)

trunk/psastro/src/pmAstrom.c

@@ -54 +54 @@
 }
 
-pmAstromMatchedListFit (pmFPA *fpa, psArray *st1, psArray *st2, psArray *match, psMetadata *config) {
+// take two matched star lists and fit a psPlaneTransform between them
+// 
+psPlaneTransform *pmAstromMatchedListFit (psPlaneTransform *map, psArray *st1, psArray *st2, psArray *match, psMetadata *config) {
+
+    if (map == NULL) {
+        // int nX = psMetadataLookupS32 (&status, myHeader, "CHIP.NX");
+        // int nY = psMetadataLookupS32 (&status, myHeader, "CHIP.NY");
+        map = psPlaneTransform (2, 2);
+    }
+
+    psStats *stats = psStatsAlloc (CLIPPED_MEAN);
+    stats->nSigma = psMetadataLookupS32 (&status, myHeader, "CHIP.NITER");
+    stats->nSigma = psMetadataLookupS32 (&status, myHeader, "CHIP.NSIGMA");
 
     psVector *X = psVectorAlloc (match->n);
@@ -62 +74 @@
   
     // take the matched stars, first fit 
-    for (i = 0; i < match->n; i++) {
+    for (int i = 0; i < match->n; i++) {
 
         pair = match->data[i];
@@ -68 +80 @@
         ob2 = st2->data[pair->i2];
 
-        X->data.F64[i] = ob1->P;
-        Y->data.F64[i] = ob1->Q;
+        X->data.F64[i] = ob1->chip.x;
+        Y->data.F64[i] = ob1->chip.y;
 
-        x->data.F64[i] = ob2->P;
-        y->data.F64[i] = ob2->Q;
-
-        // use one or the other...
-        psPolynomial2D *xt = psVectorFitPolynomial2D (NULL, stats, NULL, 0, X, NULL, x, y);
-        psPolynomial2D *xt = psVectorClipFitPolynomial2D (NULL, stats, NULL, 0, X, NULL, x, y);
-
-        psPolynomial2D *yt = psVectorFitPolynomial2D (NULL, stats, NULL, 0, Y, NULL, x, y);
-        psPolynomial2D *yt = psVectorClipFitPolynomial2D (NULL, stats, NULL, 0, Y, NULL, x, y);
+        x->data.F64[i] = ob2->FP.x;
+        y->data.F64[i] = ob2->FP.y;
     }
 
-    // apply fitted polynomials to fpa
+    // no masking, no errors
+    psVectorClipFitPolynomial2D (map->x, stats, NULL, 0, X, NULL, x, y);
+    psVectorClipFitPolynomial2D (map->y, stats, NULL, 0, Y, NULL, x, y);
+    psFree (x);
+    psFree (y);
+    psFree (X);
+    psFree (Y);
+    psFree (stats);
+
+    return (map);
 }
 
@@ -103 +117 @@
         newObj = pmAstromObjCopy (oldObj);
 
-        P = oldObj->FP.x - xCenter;
-        Q = oldObj->FP.y - yCenter;
+        X = oldObj->FP.x - xCenter;
+        Y = oldObj->FP.y - yCenter;
         
-        newObj->FP.x = P*cs + Q*sn;
-        newObj->FP.y = Q*cs - P*sn;
+        newObj->FP.x = X*cs + Y*sn;
+        newObj->FP.y = Y*cs - X*sn;
         
         new->data[i] = newObj;

trunk/psastro/src/pmAstrom.h

@@ -7 +7 @@
 
 typedef struct {
-    double X, Y;
-    double P, Q;
-    double L, M;
-    double R, D;
+    psPlane pix;
+    psPlane cell;
+    psPlane chip;
+    psPlane FP;
+    psPlane TP;
+    psSphere sky;
     double Mag, dMag;
 } pmAstromObj;
@@ -26 +28 @@
 
 typedef struct {
+    psPlane center;
+    psPlane offset;
     double angle;
-    double dP;
-    double dQ;
     double minMetric;
     double minVar;

trunk/psastro/src/pmAstromGrid.c

@@ -47 +47 @@
         st2r = pmAstromRotateObj (st2, center, angle);
         newStat = pmAstromGridMatchAngle (st1, st2r, config);
-        newStat.angle = angle;
+        newStat.angle  = angle;
+        newStat.center = center;
         if (newStat.minMetric < minStat.minMetric) {
             minStat = newStat;
@@ -147 +148 @@
                     minMetric = metric;
                     minVar    = var;
-                    minP      = i;
-                    minQ      = j;
+                    minX      = i;
+                    minY      = j;
                 }
             }
@@ -154 +155 @@
 
         // convert the bin to delta-delta
-        matchStats.dX        = DP[minY][minX] / NP[minY][minX];
-        matchStats.dY        = DQ[minY][minX] / NP[minY][minX];
+        matchStats.offset.x  = DP[minY][minX] / NP[minY][minX];
+        matchStats.offset.y  = DQ[minY][minX] / NP[minY][minX];
         matchStats.minMetric = minMetric;
         matchStats.minVar    = minVar;
@@ -164 +165 @@
 
 // apply the measured FPA offset and rotation (stat) to the fpa astrom structures
-psFPA *pmAstromApplyGridMatch (psFPA *fpa, pmAstromGridMatchStat stat) {
+psFPA *pmAstromGridApply (psPlaneTransform *map, pmAstromGridMatchStat stat) {
 
-    // stat.angle modifies fpa.toTangentPlane (effective angle)
-    // stat.dP, stat.dQ modifies fpa.projection (Ro, Do)
+    // stat.angle, stat.center, stat.offse
+    // I think i need to know the center reference....
 
     return (fpa);

trunk/psastro/src/psastro.c

@@ -22 +22 @@
     psastroProjectRawstars (subset);
 
-    // load the corresponding reference data (DVO command)
-    psArray *refstars = psastroLoadReference (argv[3]);
-
-    // use the header & config info to project refstars on the focal plane
-    psastroProjectRefstars (refstars, subset);
-
     // fpa and subset point to the same astrometry terms
-    psastroChipAstrom (refstars, subset);
+    psastroChipAstrom (subset, config);
 
     // write out data (cmp file)

trunk/psastro/src/psastroBuildFPA.c

@@ -136 +136 @@
     // set toChip to identity as default
     // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581)
+    // XXX EAM : define these in the config?
+    // int nX = psMetadataLookupS32 (&status, myHeader, "CHIP.NX");
+    // int nY = psMetadataLookupS32 (&status, myHeader, "CHIP.NY");
     cell->toChip   = psPlaneTransformAlloc (2, 2);
     cell->toChip->x->coeff[1][0] = 1;

trunk/psastro/src/psastroUtils.c

@@ -197 +197 @@
 }
 
-bool psastroProjectRawstars (pmFPA *fpa) {
+bool psastroProjectFPA (pmFPA *fpa, bool toSky) {
 
     for (int i = 0; i < fpa->chips->n; i++) {
@@ -205 +205 @@
             for (int k = 0; k < cell->readouts->n; k++) {
                 pmReadout *readout = cell->readouts->data[k];
-                for (int m = 0; m < readout->stars->n; m++) {
-                    pmAstromObj *star = readout->data[m];
-                    // apply readout offsets
-                    // apply cell PlaneTransform
-                    // apply chip PlaneTransform
-                    psPlaneTransformApply (star-> 
-
-            }
-            chip->cells->data[j] = cell;
-        }
-        fpa->chips->data[i] = chip;
-    }
-    return (fpa);
+                if (toSky) {
+                    psastroProjectRawstars (readout->stars, readout);
+                } else {
+                    psastroProjectRefstars (readout->stars, readout);
+                }
+            }
+        }
+    }
+    return true;
 }
  
-bool psastroChipAstrom (psArray *refstars, pmFPA *subset) {
-
-    // do this loop over the readouts 
-
-    // find initial offset / rotation
-    stat = pmAstromGridMatch (subset, refstars, config);
-
-    pmAstromModifyFPA (fpa, stat);
-
-    // use fit result to re-project rawstars
-    psastroProjectRawstars (fpa, subset);
-    psastroProjectRefstars (fpa, refstars);
+bool psastroProjectRawstars (psArray *stars, pmReadout *readout) {
+
+    pmCell *cell = readout->cell;
+    pmChip *chip = cell->chip;
+    pmFPA  *fpa  = chip->cell;
+
+    for (int i = 0; i < readout->stars->n; i++) {
+        pmAstromObj *star = readout->stars->data[i];
+        psCoordReadoutToCell (&star->cell, &star->pix, readout);
+        psCoordCellToChip (&star->chip, &star->cell, cell);
+        psCoordChipToFP (&star->FP, &star->chip, chip);
+        psCoordFPtoTP (&star->TP, &star->FP, fpa);
+        psCoordTPtoSky (&star->sky, &star->TP, fpa);
+    }
+    return true;
+}
+ 
+bool psastroProjectRefstars (psArray *stars, pmReadout *readout) {
+
+    pmCell *cell = readout->cell;
+    pmChip *chip = cell->chip;
+    pmFPA  *fpa  = chip->cell;
+
+    for (int i = 0; i < stars->n; i++) {
+        pmAstromObj *star = stars->data[i];
+
+        psCoordSkytoTP (&star->TP, &star->sky, fpa);
+        psCoordTPtoFP (&star->FP, &star->TP, fpa);
+        psCoordFPtoChip (&star->chip, &star->FP, chip);
+        psCoordChipToCell (&star->cell, &star->chip, cell);
+        psCoordReadoutToCell (&star->pix, &star->cell, readout);
+    }
+    return true;
+}
+ 
+// measure per-chip astrometry terms 
+bool psastroChipAstrom (pmFPA *fpa, psMetadata *config) {
+
+    pmFPA *raw = pmFPACopy (fpa);
+
+    // first pass: measure the per-chip solutions, modify the chip.toFPA terms
+    for (int i = 0; i < fpa->chips->n; i++) {
+        pmChip *chip = fpa->chips->data[i];
+
+        // cells->n > 1 is not yet well-defined
+        if (chip->cells->n > 1) {
+            psLogMsg ("pmSourcesReadCMP", 3, "undefined behavior for nCells > 1");
+            return false;
+        }
+
+        for (int j = 0; j < chip->cells->n; j++) {
+            pmCell *cell = chip->cells->data[j];
+
+            // readouts->n > 1 is not yet well-defined
+            if (cell->readouts->n > 1) {
+                psLogMsg ("pmSourcesReadCMP", 3, "undefined behavior for nReadouts > 1");
+                return false;
+            }
+
+            // load the corresponding reference data (DVO command)
+            psArray *refstars = psastroLoadReference (cell->header, config);
+
+            for (int k = 0; k < cell->readouts->n; k++) {
+
+                pmReadout *readout = cell->readouts->data[k];
+
+                // use the header & config info to project refstars on the focal plane
+                psastroProjectRefstars (refstars, readout);
+
+                // find initial offset / rotation
+                stat = pmAstromGridMatch (readout->stars, refstars, config);
+
+                // adjust the chip.toFPA terms only
+                pmAstromGridApply (chip->toFPA, stat);
+
+                // use fit result to re-project rawstars
+                psastroProjectRawstars (readout->stars, readout);
+                psastroProjectRefstars (refstars,       readout);
     
-    // use small radius to match stars
-    match = pmAstromRadiusMatch (rawstars, refstars, options);
-
-    // fit astrometric terms
-    output = pmAstromMatchedListFit (fpa, subset, refstars, match, options);
-}
-
+                // use small radius to match stars
+                match = pmAstromRadiusMatch (rawstars, refstars, options);
+
+                // fit astrometric terms
+                pmAstromMatchedListFit (chip->toFPA, readout->stars, refstars, match, options);
+            }
+        }
+    }
+
+    // second stage: re-normalize the chip terms, passing the 
+    // average rotation and offset values to the fpa.toSky
+            
+# if (0)
+    // code disabled for now 
+    // this is not as trivial as it seems at first:
+    // applying a new rotation to the FPA implies changes to the 
+    // (x,y) reference coordinates for the chips as well as a rotation
+    // I need to calculate both the x,y offset effect and the rotation
+    // in one step.  what about the effect of higher order terms in either
+    // the chip.toFPA or the fpa.toTPA? 
+    // this calculation also needs to account for scale effects between the 
+    // fpa and the chips.
+
+    // calculate the average rotation and boresite offset relative to raw
+    for (int i = 0; i < fpa->chips->n; i++) {
+        pmChip *iChip = raw->chips->data[i];
+        pmChip *oChip = fpa->chips->data[i];
+
+        T1 = psPlaneTransformGetRotation (iChip->toFPA);
+        T2 = psPlaneTransformGetRotation (oChip->toFPA);
+        dT += T1 - T2;
+        dN ++;
+
+        
+    }
+
+    dT /= dN;
+
+    psPlaneTransformSetRotation (fpa->toTPA, dT);
+        
+    // apply the new modifcations in rotation and boresite
+    for (int i = 0; i < fpa->chips->n; i++) {
+        pmChip *iChip = raw->chips->data[i];
+        pmChip *oChip = fpa->chips->data[i];
+
+        psPlaneTransformSetRotation (oChip->toFPA, -dT);
+    }
+# endif 
+}