Changeset 7014 for trunk/psastro/src/psastroWCS.c

Timestamp:

Apr 30, 2006, 12:15:03 PM (20 years ago)

Author:

eugene

Message:

first complete working version (chip only)

File:

• trunk/psastro/src/psastroWCS.c (modified) (7 diffs)

trunk/psastro/src/psastroWCS.c

@@ -2 +2 @@
 
 // interpret header WCS (only handles traditional WCS for the moment)
-bool pmAstromReadWCS (psPlaneTransform **toFPAOut, psProjection **toSkyOut, psMetadata *header) { 
+// plateScale is nominal physical scale on tangent plane (microns / arcsecond)
+bool pmAstromReadWCS (psProjection **toSkyOut, psPlaneDistort **toTPAout, psPlaneTransform **toFPAout, psMetadata *header, double plateScale) { 
 
     psPlaneTransform *toFPA;
+    psPlaneDistort *toTPA;
     psProjection *toSky;
     psProjectionType type;
@@ -19 +21 @@
 
     // determine projection type
+    // XXX add the Elixir DIS / WRP two-layer projection here
     type = PS_PROJ_NTYPE;
     if (!strcmp (&ctype[4], "-SIN")) type = PS_PROJ_SIN;
@@ -72 +75 @@
         
         // renormalize to cdelt1, cdelt2, etc
-        float scale = hypot (pc1_1, pc1_2);
+        double scale = hypot (pc1_1, pc1_2);
         cdelt1 = cdelt2 = scale;
         pc1_1 /= scale;
@@ -89 +92 @@
     //           scale = scale(i)/scale(0) (i == chip #)
     //           project crval1(0),crval2(0 using projection
-
-    // XXX EAM : psPlaneTransformAlloc uses nTerm not nOrder (bug 581)
-    // XXX EAM : I've fixed this in pslib eam_rel8_b2
     toFPA = psPlaneTransformAlloc (1, 1);
     
-    toFPA->x->coeff[0][0] = crpix1;
-    toFPA->x->coeff[1][0] = pc1_1;
-    toFPA->x->coeff[0][1] = pc1_2;
+    double cdelt = hypot (cdelt1, cdelt2) / plateScale;  // degrees / micron (eg, in fact, whatever unit user chooses for focal plane)
+    cdelt1 /= cdelt;
+    cdelt2 /= cdelt;
+
+    toFPA->x->coeff[0][0] = -(pc1_1*cdelt1*crpix1 + pc1_2*cdelt2*crpix2);
+    toFPA->x->coeff[1][0] = +(pc1_1*cdelt1);
+    toFPA->x->coeff[0][1] = +(pc1_2*cdelt2);
     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->coeff[0][0] = -(pc2_1*cdelt1*crpix1 + pc2_2*cdelt2*crpix2);
+    toFPA->y->coeff[1][0] = +(pc2_1*cdelt1);
+    toFPA->y->coeff[0][1] = +(pc2_2*cdelt2);
     toFPA->y->mask[1][1]  = 1;
-
-    // center of projection is (0,0) coordinate of TPA
-    toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, cdelt1*RAD_DEG, cdelt2*RAD_DEG, type);
-
-    *toFPAOut = toFPA;
-    *toSkyOut = toSky;
-
+    *toFPAout = toFPA;
+
+    if (*toSkyOut != NULL) {
+        if (*toTPAout == NULL) psAbort ("wcs", "projection defined, tangent-plane not defined");
+
+        psSphere sky;
+        psPlane tpa;
+
+        sky.r = crval1*RAD_DEG;
+        sky.d = crval2*RAD_DEG;
+        p_psProject (&tpa, &sky, *toSkyOut);
+        
+        // XXX for the moment, assume toTPA is the identity transformation
+        toFPA->x->coeff[0][0] = tpa.x;
+        toFPA->y->coeff[0][0] = tpa.y;
+    } else {
+        // XXX for now, use the identity for TPA <--> FPA
+        toTPA = psPlaneDistortIdentity ();
+
+        // center of projection is (0,0) coordinate of TPA
+        toSky = psProjectionAlloc (crval1*RAD_DEG, crval2*RAD_DEG, RAD_DEG*cdelt, RAD_DEG*cdelt, type);
+
+        *toTPAout = toTPA;
+        *toSkyOut = toSky;
+    }
     return true;
 }
@@ -115 +137 @@
 
 // convert toFPA / toSky components to traditional WCS
-bool pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header) { 
+// plateScale is nominal physical scale on tangent plane (microns / arcsecond)
+bool pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header, double plateScale) { 
 
     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");
+        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---SIN");
+        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", 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");
+        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---TAN");
+        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", 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");
+        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---AIT");
+        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", 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");
+        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE1", PS_META_REPLACE, "", "RA---PAR");
+        psMetadataAddStr (header, PS_LIST_TAIL, "CTYPE2", PS_META_REPLACE, "", "DEC--PAR");
         break;
       default:
@@ -139 +162 @@
     }
 
-    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]);
+    // XXX not really right: needs to deal with non-identity toTP coeffs
+    // XXX actually, totally wrong.  fix the conversions
+    // XXX need to handle the plateScale
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1",   PS_META_REPLACE, "", toSky->R*DEG_RAD);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL2",   PS_META_REPLACE, "", toSky->D*DEG_RAD);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1",   PS_META_REPLACE, "", toFPA->x->coeff[0][0]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2",   PS_META_REPLACE, "", toFPA->y->coeff[0][0]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1",   PS_META_REPLACE, "", toSky->Xs*DEG_RAD);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2",   PS_META_REPLACE, "", toSky->Ys*DEG_RAD);
+
+    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", toFPA->x->coeff[1][0]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", toFPA->x->coeff[0][1]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]);
 
     // alternative representations use
@@ -157 +183 @@
     return true;
 }
+
+psPlaneDistort *psPlaneDistortIdentity () {
+
+    psPlaneDistort *distort;
+
+    distort = psPlaneDistortAlloc (1, 1, 0, 0);
+    distort->x->coeff[0][0][0][0] = 0;
+    distort->x->coeff[1][0][0][0] = 1;
+    distort->x->coeff[0][1][0][0] = 0;
+    distort->x->mask [1][1][0][0] = 1;
+
+    distort->y->coeff[0][0][0][0] = 0;
+    distort->y->coeff[1][0][0][0] = 0;
+    distort->y->coeff[0][1][0][0] = 1;
+    distort->y->mask [1][1][0][0] = 1;
+
+    return distort;
+}