Changeset 9732 for trunk/psastro/src/psastroUtils.c

Timestamp:

Oct 24, 2006, 1:58:13 PM (20 years ago)

Author:

Paul Price

Message:

Changed definition of psArrayAlloc.

File:

• trunk/psastro/src/psastroUtils.c (modified) (10 diffs)

trunk/psastro/src/psastroUtils.c

@@ -10 +10 @@
 
     for (int i = 0; i < rawstars->n; i++) {
-        pmAstromObj *raw = rawstars->data[i];
-        
-        psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip);
-        psPlaneDistortApply (raw->TP, fpa->toTangentPlane, raw->FP, 0.0, 0.0);
-        p_psDeproject (raw->sky, raw->TP, fpa->projection);
+        pmAstromObj *raw = rawstars->data[i];
+
+        psPlaneTransformApply (raw->FP, chip->toFPA, raw->chip);
+        psPlaneDistortApply (raw->TP, fpa->toTangentPlane, raw->FP, 0.0, 0.0);
+        p_psDeproject (raw->sky, raw->TP, fpa->projection);
     }
 
     for (int i = 0; i < refstars->n; i++) {
-        pmAstromObj *ref = refstars->data[i];
-        psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
+        pmAstromObj *ref = refstars->data[i];
+        psPlaneTransformApply (ref->chip, chip->fromFPA, ref->FP);
     }
 
@@ -35 +35 @@
 
     for (int i = 0; i < fpa->chips->n; i++) {
-        pmChip *chip = fpa->chips->data[i];
-        for (int j = 0; j < chip->cells->n; j++) {
-            pmCell *cell = chip->cells->data[j];
-            for (int k = 0; k < cell->readouts->n; k++) {
-                pmReadout *readout = cell->readouts->data[k];
-
-                psArray *stars = psMetadataLookupPtr (&status, readout->analysis, "STARS.FULLSET");
-                stars = psArraySort (stars, pmAstromObjSortByMag);
-
-                int nSubset = PS_MIN (MAX_NSTARS, stars->n);
-                psArray *subset = psArrayAlloc (nSubset);
-
-                for (int i = 0; i < nSubset; i++) {
-                    subset->data[i] = stars->data[i];
-                }
-                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "STARS.SUBSET", PS_DATA_ARRAY, "stars from analysis", subset);
-            }
-        }
+        pmChip *chip = fpa->chips->data[i];
+        for (int j = 0; j < chip->cells->n; j++) {
+            pmCell *cell = chip->cells->data[j];
+            for (int k = 0; k < cell->readouts->n; k++) {
+                pmReadout *readout = cell->readouts->data[k];
+
+                psArray *stars = psMetadataLookupPtr (&status, readout->analysis, "STARS.FULLSET");
+                stars = psArraySort (stars, pmAstromObjSortByMag);
+
+                int nSubset = PS_MIN (MAX_NSTARS, stars->n);
+                psArray *subset = psArrayAlloc (nSubset);
+
+                for (int i = 0; i < nSubset; i++) {
+                    subset->data[i] = stars->data[i];
+                }
+                psMetadataAdd (readout->analysis, PS_LIST_TAIL, "STARS.SUBSET", PS_DATA_ARRAY, "stars from analysis", subset);
+            }
+        }
     }
     return true;
@@ -64 +64 @@
     // 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];
-        psastroChipAstrom (chip, config);
-    }
-
-    // second stage: re-normalize the chip terms, passing the 
+        pmChip *chip = fpa->chips->data[i];
+        psastroChipAstrom (chip, config);
+    }
+
+    // second stage: re-normalize the chip terms, passing the
     // average rotation and offset values to the fpa.toSky
     if (RENORM) {
 
-        // this code is needed for the mosastro stage, with multiple chip solutions
-
-        double dX, dY, dT, dN;
-        dX = dY = dT = dN = 0;
-
-        psPlane origin, P1, P2;
-        origin.x = 0;
-        origin.y = 0;
-
-        // 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];
-
-            // offset of chip
-            psCoordChipToFPA (&P1, &origin, iChip);
-            psCoordChipToFPA (&P2, &origin, oChip);
-            dX += (P2.x - P1.x);
-            dY += (P2.y - P1.y);
-
-            // get parity-independent rotations for old and new solutions
-            double T1 = psPlaneTransformGetRotation (iChip->toFPA);
-            double T2 = psPlaneTransformGetRotation (oChip->toFPA);
-            dT += T2 - T1;
-            dN ++;
-        }
-
-        dT /= dN;
-        dX /= dN;
-        dY /= dN;
-
-        // R(T)
-        double PC1_1 = fpa->toTPA->x->coeff[1][0][0][0];
-        double PC1_2 = fpa->toTPA->x->coeff[0][1][0][0];
-        double PC2_1 = fpa->toTPA->y->coeff[1][0][0][0];
-        double PC2_2 = fpa->toTPA->y->coeff[0][1][0][0];
-
-        // R(dT)
-        double dPC1_1 = +cos (dT);
-        double dPC1_2 = +sin (dT);
-        double dPC2_1 = -sin (dT);
-        double dPC2_2 = +cos (dT);
-
-        // R'(T) = R(T) * R(dT)
-        double pc1_1 = PC1_1*dPC1_1 + PC1_2*dPC2_1;
-        double pc1_2 = PC1_1*dPC1_2 + PC1_2*dPC2_2;
-        double pc2_1 = PC2_1*dPC1_1 + PC2_2*dPC2_1;
-        double pc2_2 = PC2_1*dPC1_2 + PC2_2*dPC2_2;
-
-        double det = 1.0 / (pc1_1*pc2_2 - pc1_2*pc2_1);
-
-        // R'(-T)  (matrix inverse, not just rotation inverse -- keeps parity)
-        double pi1_1 = +pc2_2 * det;
-        double pi1_2 = -pc1_2 * det;
-        double pi2_1 = -pc2_1 * det;
-        double pi2_2 = +pc1_1 * det;
-    
-        // apply the new modifcations in rotation and boresite
-        for (int i = 0; i < fpa->chips->n; i++) {
-            pmChip *oChip = fpa->chips->data[i];
-
-            // r(T)
-            double pr1_1 = oChip->toFPA->x->coeff[1][0];
-            double pr1_2 = oChip->toFPA->x->coeff[0][1];
-            double pr2_1 = oChip->toFPA->y->coeff[1][0];
-            double pr2_2 = oChip->toFPA->y->coeff[0][1];
-
-            // ri'(T) = R(T) r(t)
-            double ri1_1 = PC1_1*pr1_1 + PC1_2*pr2_1;
-            double ri1_2 = PC1_1*pr1_2 + PC1_2*pr2_2;
-            double ri2_1 = PC2_1*pr1_1 + PC2_2*pr2_1;
-            double ri2_2 = PC2_1*pr1_2 + PC2_2*pr2_2;
-
-            // r'(T) = R'(-T) ri'(T)
-            oChip->toFPA->x->coeff[1][0] = pi1_1*ri1_1 + pi1_2*ri2_1;
-            oChip->toFPA->x->coeff[0][1] = pi1_1*ri1_2 + pi1_2*ri2_2;
-            oChip->toFPA->y->coeff[1][0] = pi2_1*ri1_1 + pi2_2*ri2_1;
-            oChip->toFPA->y->coeff[0][1] = pi2_1*ri1_2 + pi2_2*ri2_2;
-
-            double dx = PC1_1*oChip->toFPA->x->coeff[0][0] + PC1_2*oChip->toFPA->y->coeff[0][0] + dX;
-            double dy = PC2_1*oChip->toFPA->x->coeff[0][0] + PC2_2*oChip->toFPA->y->coeff[0][0] + dY;
-
-            oChip->toFPA->x->coeff[0][0] = pi1_1*dx + pi1_2*dy;
-            oChip->toFPA->y->coeff[0][0] = pi2_1*dx + pi2_2*dy;
-        }
-
-        fpa->toTPA->x->coeff[0][0][0][0] -= dX;
-        fpa->toTPA->y->coeff[0][0][0][0] -= dY;
-
-        fpa->toTPA->x->coeff[1][0][0][0] = pc1_1;
-        fpa->toTPA->x->coeff[0][1][0][0] = pc1_2;
-        fpa->toTPA->y->coeff[1][0][0][0] = pc2_1;
-        fpa->toTPA->y->coeff[0][1][0][0] = pc2_2;
-    }
-    return true; 
+        // this code is needed for the mosastro stage, with multiple chip solutions
+
+        double dX, dY, dT, dN;
+        dX = dY = dT = dN = 0;
+
+        psPlane origin, P1, P2;
+        origin.x = 0;
+        origin.y = 0;
+
+        // 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];
+
+            // offset of chip
+            psCoordChipToFPA (&P1, &origin, iChip);
+            psCoordChipToFPA (&P2, &origin, oChip);
+            dX += (P2.x - P1.x);
+            dY += (P2.y - P1.y);
+
+            // get parity-independent rotations for old and new solutions
+            double T1 = psPlaneTransformGetRotation (iChip->toFPA);
+            double T2 = psPlaneTransformGetRotation (oChip->toFPA);
+            dT += T2 - T1;
+            dN ++;
+        }
+
+        dT /= dN;
+        dX /= dN;
+        dY /= dN;
+
+        // R(T)
+        double PC1_1 = fpa->toTPA->x->coeff[1][0][0][0];
+        double PC1_2 = fpa->toTPA->x->coeff[0][1][0][0];
+        double PC2_1 = fpa->toTPA->y->coeff[1][0][0][0];
+        double PC2_2 = fpa->toTPA->y->coeff[0][1][0][0];
+
+        // R(dT)
+        double dPC1_1 = +cos (dT);
+        double dPC1_2 = +sin (dT);
+        double dPC2_1 = -sin (dT);
+        double dPC2_2 = +cos (dT);
+
+        // R'(T) = R(T) * R(dT)
+        double pc1_1 = PC1_1*dPC1_1 + PC1_2*dPC2_1;
+        double pc1_2 = PC1_1*dPC1_2 + PC1_2*dPC2_2;
+        double pc2_1 = PC2_1*dPC1_1 + PC2_2*dPC2_1;
+        double pc2_2 = PC2_1*dPC1_2 + PC2_2*dPC2_2;
+
+        double det = 1.0 / (pc1_1*pc2_2 - pc1_2*pc2_1);
+
+        // R'(-T)  (matrix inverse, not just rotation inverse -- keeps parity)
+        double pi1_1 = +pc2_2 * det;
+        double pi1_2 = -pc1_2 * det;
+        double pi2_1 = -pc2_1 * det;
+        double pi2_2 = +pc1_1 * det;
+
+        // apply the new modifcations in rotation and boresite
+        for (int i = 0; i < fpa->chips->n; i++) {
+            pmChip *oChip = fpa->chips->data[i];
+
+            // r(T)
+            double pr1_1 = oChip->toFPA->x->coeff[1][0];
+            double pr1_2 = oChip->toFPA->x->coeff[0][1];
+            double pr2_1 = oChip->toFPA->y->coeff[1][0];
+            double pr2_2 = oChip->toFPA->y->coeff[0][1];
+
+            // ri'(T) = R(T) r(t)
+            double ri1_1 = PC1_1*pr1_1 + PC1_2*pr2_1;
+            double ri1_2 = PC1_1*pr1_2 + PC1_2*pr2_2;
+            double ri2_1 = PC2_1*pr1_1 + PC2_2*pr2_1;
+            double ri2_2 = PC2_1*pr1_2 + PC2_2*pr2_2;
+
+            // r'(T) = R'(-T) ri'(T)
+            oChip->toFPA->x->coeff[1][0] = pi1_1*ri1_1 + pi1_2*ri2_1;
+            oChip->toFPA->x->coeff[0][1] = pi1_1*ri1_2 + pi1_2*ri2_2;
+            oChip->toFPA->y->coeff[1][0] = pi2_1*ri1_1 + pi2_2*ri2_1;
+            oChip->toFPA->y->coeff[0][1] = pi2_1*ri1_2 + pi2_2*ri2_2;
+
+            double dx = PC1_1*oChip->toFPA->x->coeff[0][0] + PC1_2*oChip->toFPA->y->coeff[0][0] + dX;
+            double dy = PC2_1*oChip->toFPA->x->coeff[0][0] + PC2_2*oChip->toFPA->y->coeff[0][0] + dY;
+
+            oChip->toFPA->x->coeff[0][0] = pi1_1*dx + pi1_2*dy;
+            oChip->toFPA->y->coeff[0][0] = pi2_1*dx + pi2_2*dy;
+        }
+
+        fpa->toTPA->x->coeff[0][0][0][0] -= dX;
+        fpa->toTPA->y->coeff[0][0][0][0] -= dY;
+
+        fpa->toTPA->x->coeff[1][0][0][0] = pc1_1;
+        fpa->toTPA->x->coeff[0][1][0][0] = pc1_2;
+        fpa->toTPA->y->coeff[1][0][0][0] = pc2_1;
+        fpa->toTPA->y->coeff[0][1][0][0] = pc2_2;
+    }
+    return true;
 }
 
@@ -174 +174 @@
 
     for (int i = 0; i < input->nX; i++) {
-        for (int j = 0; j < input->nY; j++) {
-            output->mask[i][j]     = input->mask[i][j];
-            output->coeff[i][j]    = input->coeff[i][j];
-            output->coeffErr[i][j] = input->coeffErr[i][j];
-        }
+        for (int j = 0; j < input->nY; j++) {
+            output->mask[i][j]     = input->mask[i][j];
+            output->coeff[i][j]    = input->coeff[i][j];
+            output->coeffErr[i][j] = input->coeffErr[i][j];
+        }
     }
     return (output);
@@ -188 +188 @@
 
     for (int i = 0; i < input->nX; i++) {
-        for (int j = 0; j < input->nY; j++) {
-            for (int k = 0; k < input->nZ; k++) {
-                for (int m = 0; m < input->nT; m++) {
-                    output->mask[i][j][k][m]     = input->mask[i][j][k][m];
-                    output->coeff[i][j][k][m]    = input->coeff[i][j][k][m];
-                    output->coeffErr[i][j][k][m] = input->coeffErr[i][j][k][m];
-                }
-            }
-        }
+        for (int j = 0; j < input->nY; j++) {
+            for (int k = 0; k < input->nZ; k++) {
+                for (int m = 0; m < input->nT; m++) {
+                    output->mask[i][j][k][m]     = input->mask[i][j][k][m];
+                    output->coeff[i][j][k][m]    = input->coeff[i][j][k][m];
+                    output->coeffErr[i][j][k][m] = input->coeffErr[i][j][k][m];
+                }
+            }
+        }
     }
     return (output);
@@ -206 +206 @@
 
     for (int i = 0; i < input->x->nX; i++) {
-        for (int j = 0; j < input->x->nY; j++) {
-            for (int k = 0; k < input->x->nZ; k++) {
-                for (int m = 0; m < input->x->nT; m++) {
-                    // x-terms
-                    output->x->mask[i][j][k][m]     = input->x->mask[i][j][k][m];
-                    output->x->coeff[i][j][k][m]    = input->x->coeff[i][j][k][m];
-                    output->x->coeffErr[i][j][k][m] = input->x->coeffErr[i][j][k][m];
-                    // y-terms
-                    output->y->mask[i][j][k][m]     = input->y->mask[i][j][k][m];
-                    output->y->coeff[i][j][k][m]    = input->y->coeff[i][j][k][m];
-                    output->y->coeffErr[i][j][k][m] = input->y->coeffErr[i][j][k][m];
-                }
-            }
-        }
+        for (int j = 0; j < input->x->nY; j++) {
+            for (int k = 0; k < input->x->nZ; k++) {
+                for (int m = 0; m < input->x->nT; m++) {
+                    // x-terms
+                    output->x->mask[i][j][k][m]     = input->x->mask[i][j][k][m];
+                    output->x->coeff[i][j][k][m]    = input->x->coeff[i][j][k][m];
+                    output->x->coeffErr[i][j][k][m] = input->x->coeffErr[i][j][k][m];
+                    // y-terms
+                    output->y->mask[i][j][k][m]     = input->y->mask[i][j][k][m];
+                    output->y->coeff[i][j][k][m]    = input->y->coeff[i][j][k][m];
+                    output->y->coeffErr[i][j][k][m] = input->y->coeffErr[i][j][k][m];
+                }
+            }
+        }
     }
     return (output);
@@ -229 +229 @@
 
     for (int i = 0; i < input->x->nX; i++) {
-        for (int j = 0; j < input->x->nY; j++) {
-            // x-terms
-            output->x->mask[i][j]     = input->x->mask[i][j];
-            output->x->coeff[i][j]    = input->x->coeff[i][j];
-            output->x->coeffErr[i][j] = input->x->coeffErr[i][j];
-            // y-terms
-            output->y->mask[i][j]     = input->y->mask[i][j];
-            output->y->coeff[i][j]    = input->y->coeff[i][j];
-            output->y->coeffErr[i][j] = input->y->coeffErr[i][j];
-        }
+        for (int j = 0; j < input->x->nY; j++) {
+            // x-terms
+            output->x->mask[i][j]     = input->x->mask[i][j];
+            output->x->coeff[i][j]    = input->x->coeff[i][j];
+            output->x->coeffErr[i][j] = input->x->coeffErr[i][j];
+            // y-terms
+            output->y->mask[i][j]     = input->y->mask[i][j];
+            output->y->coeff[i][j]    = input->y->coeff[i][j];
+            output->y->coeffErr[i][j] = input->y->coeffErr[i][j];
+        }
     }
     return (output);
@@ -287 +287 @@
     if (map->y->nX < 1) return 0;
     if (map->y->nY < 1) return 0;
-    
+
     double pc1_1 = map->x->coeff[1][0];
     double pc1_2 = map->x->coeff[0][1];
@@ -299 +299 @@
     double t1 = -atan2 (px*pc1_2, px*pc1_1);
     double t2 = +atan2 (py*pc2_1, py*pc2_2);
-    
+
     // careful near -pi,+pi boundary...
     if (t1 - t2 > M_PI/2) t2 += 2*M_PI;
@@ -307 +307 @@
     while (theta < M_PI) theta += 2*M_PI;
     while (theta > M_PI) theta -= 2*M_PI;
-    
+
     return (theta);
 }