Index: trunk/psastro/src/psastroWCS.c
===================================================================
--- trunk/psastro/src/psastroWCS.c	(revision 6176)
+++ trunk/psastro/src/psastroWCS.c	(revision 7014)
@@ -2,7 +2,9 @@
 
 // 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,4 +21,5 @@
 
     // 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,5 +75,5 @@
 	
 	// 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,25 +92,44 @@
     //           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,22 +137,23 @@
 
 // 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,15 +162,18 @@
     }
 
-    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,2 +183,20 @@
     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;
+}