Changeset 5575 for trunk/psastro/src/psastroBuildFPA.c

Timestamp:

Nov 22, 2005, 10:38:21 AM (21 years ago)

Author:

eugene

Message:

fixed up WCS

File:

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

trunk/psastro/src/psastroBuildFPA.c

@@ -20 +20 @@
     pmFPA *fpa = pmFPAAlloc (NULL, NULL);
 
+    // identify matrix for fpa to tpa
+    fpa->toTPA   = psPlaneDistortAlloc (1, 1, 0, 0);
+    fpa->toTPA->x->coeff[1][0][0][0] = 1;
+    fpa->toTPA->x->mask[1][1][0][0]  = 1;
+
+    fpa->toTPA->y->coeff[0][1][0][0] = 1;
+    fpa->toTPA->y->mask[1][1][0][0]  = 1;
+    fpa->fromTangentPlane = psPlaneDistortInvert (fpa->toTangentPlane);
+
     psArrayRealloc (fpa->chips, Nchips);
     for (int i = 0; i < Nchips; i++) {
@@ -30 +39 @@
             // XXX EAM : extend this function to take more than one header
             cell->header = header;      
-            pmCellInterpretWCS (cell, header);
+            pmAstromReadWCS (&chip->toFPA, &fpa->toSky, header);
+            chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA);
+
+            // allocate identity matrix for cell to chip
+            cell->toChip   = psPlaneTransformAlloc (1, 1);
+            cell->toChip->x->coeff[1][0] = 1;
+            cell->toChip->x->mask[1][1]  = 1;
+            cell->toChip->y->coeff[0][1] = 1;
+            cell->toChip->y->mask[1][1]  = 1;
 
             psArrayRealloc (cell->readouts, Nreadouts);
@@ -52 +69 @@
 }
 
-// interpret header WCS
-bool pmCellInterpretWCS (pmCell *cell, psMetadata *header) { 
-
-    pmChip *chip;
-    pmFPA  *fpa;
+// interpret header WCS (only handles traditional WCS for the moment)
+bool pmAstromReadWCS (psPlaneTransform **toFPAOut, psProjection **toSkyOut, psMetadata *header) { 
+
+    psPlaneTransform *toFPA;
+    psProjection *toSky;
+    psProjectionType type;
     bool status;
-    psProjectionType type;
     float crval1, crval2, crpix1, crpix2, cdelt1, cdelt2;
     float pc1_1, pc1_2, pc2_1, pc2_2;
@@ -136 +153 @@
 got_matrix:
 
-    // XXX EAM : these must already be set
-    chip = cell->parent;
-    fpa  = chip->parent;
-
-    // set toChip to identity as default
-    // XXX EAM : can we actually use higher order?
-    int nX = psMetadataLookupS32 (&status, header, "PSASTRO.CHIP.NX");
-    if (!status) nX = 1;
-    int nY = psMetadataLookupS32 (&status, header, "PSASTRO.CHIP.NY");
-    if (!status) nY = 1;
-
-    cell->toChip   = psPlaneTransformAlloc (1, 1);
-    cell->toChip->x->coeff[1][0] = 1;
-    cell->toChip->x->mask[1][1]  = 1;
-
-    cell->toChip->y->coeff[0][1] = 1;
-    cell->toChip->y->mask[1][1]  = 1;
-
     // XXX EAM : if fpa->toSky and fpa->toTPA are already defined, then the
     //           toFPA must be modified to match the crval(i), scale(i) and crpix(i)
@@ -161 +160 @@
     // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581)
     // XXX EAM : I've fixed this in pslib eam_rel8_b2
-    chip->toFPA   = psPlaneTransformAlloc (1, 1);
+    toFPA = psPlaneTransformAlloc (1, 1);
     
-    chip->toFPA->x->coeff[0][0] = crpix1;
-    chip->toFPA->x->coeff[1][0] = pc1_1;
-    chip->toFPA->x->coeff[0][1] = pc1_2;
-    chip->toFPA->x->mask[1][1]  = 1;
-
-    chip->toFPA->y->coeff[0][0] = crpix2;
-    chip->toFPA->y->coeff[1][0] = pc2_1;
-    chip->toFPA->y->coeff[0][1] = pc2_2;
-    chip->toFPA->y->mask[1][1]  = 1;
-
-    chip->fromFPA = p_psPlaneTransformLinearInvert (chip->toFPA);
-
-    // set toTPA to identity as default
-    // XXX EAM : psPlaneDistortAlloc uses nTerm not nOrder (bug 581)
-    if (fpa->toTPA == NULL) {
-        fpa->toTPA   = psPlaneDistortAlloc (1, 1, 0, 0);
-        fpa->toTPA->x->coeff[1][0][0][0] = 1;
-        fpa->toTPA->x->mask[1][1][0][0]  = 1;
-
-        fpa->toTPA->y->coeff[0][1][0][0] = 1;
-        fpa->toTPA->y->mask[1][1][0][0]  = 1;
-        fpa->fromTangentPlane = psPlaneDistortInvert (fpa->toTangentPlane);
-    } else {
-        psLogMsg ("psastro", 2, "warning: fpa distortion already defined\n");
-    }
+    toFPA->x->coeff[0][0] = crpix1;
+    toFPA->x->coeff[1][0] = pc1_1;
+    toFPA->x->coeff[0][1] = pc1_2;
+    toFPA->x->mask[1][1]  = 1;
+
+    toFPA->y->coeff[0][0] = crpix2;
+    toFPA->y->coeff[1][0] = pc2_1;
+    toFPA->y->coeff[0][1] = pc2_2;
+    toFPA->y->mask[1][1]  = 1;
 
     // center of projection is (0,0) coordinate of TPA
-    fpa->toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type);
+    toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type);
+
+    *toFPAOut = toFPA;
+    *toSkyOut = toSky;
+
     return true;
 }
+
+
+// convert toFPA / toSky components to traditional WCS
+bool pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header) { 
+
+    switch (toSky->type) {
+      case PS_PROJ_SIN:
+        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---SIN");
+        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--SIN");
+        break;
+      case PS_PROJ_TAN:
+        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---TAN");
+        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--TAN");
+        break;
+      case PS_PROJ_AIT:
+        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---AIT");
+        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--AIT");
+        break;
+      case PS_PROJ_PAR:
+        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE1", PS_DATA_STRING | PS_META_REPLACE, "", "RA---PAR");
+        psMetadataAdd (header, PS_LIST_TAIL, "CTYPE2", PS_DATA_STRING | PS_META_REPLACE, "", "DEC--PAR");
+        break;
+      default:
+        psLogMsg ("psastro", 2, "warning: unknown projection type %s\n", toSky->type);
+        return false;
+    }
+
+    psMetadataAdd (header, PS_LIST_TAIL, "CRVAL1", PS_DATA_F32 | PS_META_REPLACE, "", toSky->R*DEG_RAD);
+    psMetadataAdd (header, PS_LIST_TAIL, "CRVAL2", PS_DATA_F32 | PS_META_REPLACE, "", toSky->D*DEG_RAD);
+    psMetadataAdd (header, PS_LIST_TAIL, "CRPIX1", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[0][0]);
+    psMetadataAdd (header, PS_LIST_TAIL, "CRPIX2", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[0][0]);
+    psMetadataAdd (header, PS_LIST_TAIL, "CDELT1", PS_DATA_F32 | PS_META_REPLACE, "", toSky->Xs*DEG_RAD);
+    psMetadataAdd (header, PS_LIST_TAIL, "CDELT2", PS_DATA_F32 | PS_META_REPLACE, "", toSky->Ys*DEG_RAD);
+
+    psMetadataAdd (header, PS_LIST_TAIL, "PC001001", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[1][0]);
+    psMetadataAdd (header, PS_LIST_TAIL, "PC001002", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->x->coeff[0][1]);
+    psMetadataAdd (header, PS_LIST_TAIL, "PC002001", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[1][0]);
+    psMetadataAdd (header, PS_LIST_TAIL, "PC002002", PS_DATA_F32 | PS_META_REPLACE, "", toFPA->y->coeff[0][1]);
+
+    // alternative representations use
+    // CD1_1 = PC001001*CDELT1, etc
+    // make these representations optional
+
+    return true;
+}