Index: /trunk/psastro/src/psastro.h
===================================================================
--- /trunk/psastro/src/psastro.h	(revision 10512)
+++ /trunk/psastro/src/psastro.h	(revision 10513)
@@ -43,6 +43,7 @@
 
 bool              pmAstromReadWCS (pmFPA *fpa, pmChip *chip, psMetadata *header, double plateScale, bool isMosaic);
-bool              pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header, double plateScale);
+bool              pmAstromWriteWCS (psPlaneTransform *toFPA, psPlaneDistort *toTPA, psProjection *toSky, psMetadata *header, double plateScale);
 psPlaneDistort   *psPlaneDistortIdentity ();
+bool psPlaneDistortIsIdentity (psPlaneDistort *distort);
 
 psArray          *psastroLoadRefstars (pmConfig *config);
Index: /trunk/psastro/src/psastroChipAstrom.c
===================================================================
--- /trunk/psastro/src/psastroChipAstrom.c	(revision 10512)
+++ /trunk/psastro/src/psastroChipAstrom.c	(revision 10513)
@@ -55,5 +55,5 @@
 		    return false;
 		}
-		pmAstromWriteWCS (chip->toFPA, fpa->toSky, updates, plateScale);
+		pmAstromWriteWCS (chip->toFPA, fpa->toTangentPlane, fpa->toSky, updates, plateScale);
 		psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_DATA_METADATA, "psastro header stats", updates);
 		psFree (updates);
Index: /trunk/psastro/src/psastroWCS.c
===================================================================
--- /trunk/psastro/src/psastroWCS.c	(revision 10512)
+++ /trunk/psastro/src/psastroWCS.c	(revision 10513)
@@ -208,5 +208,9 @@
 // convert toFPA / toSky components to traditional WCS
 // plateScale is nominal physical scale on tangent plane (microns / arcsecond)
-bool pmAstromWriteWCS (psPlaneTransform *toFPA, psProjection *toSky, psMetadata *header, double plateScale) { 
+// this requires toTP to be the identity transformation
+bool pmAstromWriteWCS (psPlaneTransform *toFPA, psPlaneDistort *toTPA, psProjection *toSky, psMetadata *header, double plateScale) { 
+
+    // techinically, we can have a plate scale here (toTPA:dx,dy != 1)
+    if (!psPlaneDistortIsIdentity (toTPA)) psAbort ("psastro", "invalid TPA transformation");
 
     switch (toSky->type) {
@@ -232,7 +236,23 @@
     }
 
-    // XXX not really right: needs to deal with non-identity toTP coeffs
+# if (0)
+    // XXX not really right: needs to deal with non-identity to coeffs
     // XXX actually, totally wrong.  fix the conversions
     // XXX need to handle the plateScale
+    
+    // solve for CDELT1,2:
+    cdelt1 = toSky->Xs*DEG_RAD*toTPA->x->coeff[1][0][0][0];
+    cdelt2 = toSky->Ys*DEG_RAD*toTPA->y->coeff[0][1][0][0];
+
+    // L,M = toFPA(X,Y)
+    // solve for CRPIX1,2 (Xo,Yo) : L,M(Xo,Yo) = 0,0
+    // linear solution for Xo,Yo:
+    R  = (xsum[1][0]*ysum[0][1] - xsum[0][1]*ysum[1][0]);
+    Xo = det*(ysum[0][0]*xsum[0][1] - xsum[0][0]*ysum[0][1]);
+    Yo = det*(xsum[0][0]*ysum[1][0] - ysum[0][0]*xsum[1][0]);
+
+    if (
+# endif
+
     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);
@@ -362,5 +382,5 @@
 }
 
-// XXX for the moment, we mimic the limited Elixir mosastro orders
+// construct a psPlaneDistort which is the identify transformation
 psPlaneDistort *psPlaneDistortIdentity (int order) {
 
@@ -376,6 +396,6 @@
 	for (int j = 0; j <= order; j++) {
 	    if (i + j > order) {
-		distort->x->mask [1][1][0][0] = 1;
-		distort->y->mask [1][1][0][0] = 1;
+		distort->x->mask [i][j][0][0] = 1;
+		distort->y->mask [i][j][0][0] = 1;
 	    } 
 	}
@@ -386,2 +406,66 @@
     return distort;
 }
+
+// check that the given psPlaneDistort is the identity
+bool psPlaneDistortIsIdentity (psPlaneDistort *distort) {
+
+    int order;
+    bool status;
+
+    // we currently only support up to 3rd order polynomials
+    if (distort->x->nX < 1) return false;
+    if (distort->x->nY < 1) return false;
+    if (distort->y->nX < 1) return false;
+    if (distort->y->nY < 1) return false;
+
+    if (distort->x->nX > 3) return false;
+    if (distort->x->nY > 3) return false;
+    if (distort->y->nX > 3) return false;
+    if (distort->y->nY > 3) return false;
+    
+    if (distort->x->nZ > 0) return false;
+    if (distort->x->nT > 0) return false;
+    if (distort->y->nZ > 0) return false;
+    if (distort->y->nT > 0) return false;
+    
+    if (distort->x->nX != distort->x->nY) return false;
+    if (distort->y->nX != distort->y->nY) return false;
+
+    status = true;
+    order = distort->x->nX;
+    for (int i = 0; i <= order; i++) {
+	for (int j = 0; j <= order; j++) {
+	    if (i + j > order) {
+		// high-order cross terms must be masked (eg, x^3 y^2)
+		status &= !distort->x->mask[i][j][0][0];
+	    } else {
+		if (i + j != 1) {
+		    // non-linear and off-diagonal terms must be 0 (eg, x^2, x y)
+		    status &= (distort->x->coeff[i][j][0][0] == 0.0);
+		} else {
+		    // linear, diagonal terms must be 1.0
+		    status &= (distort->x->coeff[i][j][0][0] == 1.0);
+		}
+	    }
+	}
+    }
+
+    order = distort->y->nX;
+    for (int i = 0; i <= order; i++) {
+	for (int j = 0; j <= order; j++) {
+	    if (i + j > order) {
+		// high-order cross terms must be masked (eg, x^3 y^2)
+		status &= !distort->y->mask[i][j][0][0];
+	    } else {
+		if (i + j != 1) {
+		    // non-linear and off-diagonal terms must be 0 (eg, x^2, x y)
+		    status &= (distort->y->coeff[i][j][0][0] == 0.0);
+		} else {
+		    // linear, diagonal terms must be 1.0
+		    status &= (distort->y->coeff[i][j][0][0] == 1.0);
+		}
+	    }
+	}
+    }
+    return status;
+}