Index: /trunk/psastro/src/psastroWCS.c
===================================================================
--- /trunk/psastro/src/psastroWCS.c	(revision 10593)
+++ /trunk/psastro/src/psastroWCS.c	(revision 10594)
@@ -213,4 +213,9 @@
     // techinically, we can have a plate scale here (toTPA:dx,dy != 1)
     if (!psPlaneDistortIsIdentity (toTPA)) psAbort ("psastro", "invalid TPA transformation");
+    
+    // XXX require toFPA->x->nX == toFPA->x->nY
+    // XXX require toFPA->y->nX == toFPA->y->nY
+    // XXX require toFPA->x->nX == toFPA->y->nX
+    // XXX require toFPA->nX == 1,2,3
 
     switch (toSky->type) {
@@ -241,17 +246,155 @@
     // 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];
+    /* discussion of the coord transformations:
+    X,Y: coord on a chip in pixels
+    L,M: coord on the focal plane (pixels)
+    P,Q: coord in the tangent plane (microns or mm?)
+    R,D: coord on the sky 
+    
+    this function creates WCS terms which convert directly from chip to sky.
+    this function requires a linear, unrotated toTPA distortion term
+    toTPA->x,y->coeff[1][0],[0][1] defines the detector scale (microns / pixel)
+    tpSky->Xs,Ys defines the plate scale (radians / micron)
+    */
+    
+    // solve for CDELT1,2 (degrees / pixel)
+    cdelt1 = DEG_RAD*toSky->Xs*toTPA->x->coeff[1][0][0][0];
+    cdelt2 = DEG_RAD*toSky->Ys*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
+    xcoeff = toFPA->x->coeff;
+    ycoeff = toFPA->y->coeff;
+    R  = (xcoeff[1][0]*ycoeff[0][1] - xcoeff[0][1]*ycoeff[1][0]);
+    Xo = det*(ycoeff[0][0]*xcoeff[0][1] - xcoeff[0][0]*ycoeff[0][1]);
+    Yo = det*(xcoeff[0][0]*ycoeff[1][0] - ycoeff[0][0]*xcoeff[1][0]);
+
+    if (toFPA->x->nX > 1) {
+
+	psPolynomial2D *XdX = psPolynomial2D_dX(toFPA->x);
+	psPolynomial2D *XdY = psPolynomial2D_dY(toFPA->x);
+
+	psPolynomial2D *YdX = psPolynomial2D_dX(toFPA->y);
+	psPolynomial2D *YdY = psPolynomial2D_dY(toFPA->y);
+
+	psImage *Alpha = psImageAlloc (2, 2, PS_DATA_F32);
+	psVector *Beta = psVectorAlloc (2, PS_DATA_F32);
+
+	// XXX this loop is rather arbitrary in length...
+	for (int i = 0; i < 10; i++) {
+	    // NOTE: order is: [y][x]
+	    Alpha->data.F32[0][0] = psPolynomial2DEval (XdX, Xo, Yo);
+	    Alpha->data.F32[1][0] = psPolynomial2DEval (XdY, Xo, Yo);
+	    Alpha->data.F32[0][1] = psPolynomial2DEval (YdX, Xo, Yo);
+	    Alpha->data.F32[1][1] = psPolynomial2DEval (YdY, Xo, Yo);
+
+	    Beta->data.F32[0] = psPolynomial2DEval (toFPA->x, Xo, Yo);
+	    Beta->data.F32[1] = psPolynomial2DEval (toFPA->y, Xo, Yo);
+
+	    psMatrixGJSolveF32 (Alpha, Beta);
+	
+	    Xo += Beta->data.F32[0];
+	    Yo += Beta->data.F32[1];
+	}
+    }
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX1", 	PS_META_REPLACE, "", Xo);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", 	PS_META_REPLACE, "", Yo);
+
+    switch (toFPA->x->nX) {
+
+      case 1:
+	/* the linear solution can be analytically inverted */
+	pc1_1 = xcoeff[1][0] / cdelt1;
+	pc1_2 = xcoeff[0][1] / cdelt2;
+	pc2_1 = ycoeff[1][0] / cdelt1;
+	pc2_2 = ycoeff[0][1] / cdelt2;
+
+      case 2:
+	a10 = xcoeff[1][0] + 2.0*xcoeff[2][0]*Xo + xcoeff[1][1]*Yo;
+	a01 = xcoeff[0][1] + 2.0*xcoeff[0][2]*Yo + xcoeff[1][1]*Xo;
+	a11 = xcoeff[1][1];
+	a20 = xcoeff[2][0];
+	a02 = xcoeff[0][2];
+
+	b10 = ycoeff[1][0] + 2.0*ycoeff[2][0]*Xo + ycoeff[1][1]*Yo;
+	b01 = ycoeff[0][1] + 2.0*ycoeff[0][2]*Yo + ycoeff[1][1]*Xo;
+	b20 = ycoeff[2][0];
+	b11 = ycoeff[1][1];
+	b02 = ycoeff[0][2];
+
+	coords[0].pc1_1 = a10 / coords[0].cdelt1;
+	coords[0].pc1_2 = a01 / coords[0].cdelt2;
+	coords[0].pc2_1 = b10 / coords[0].cdelt1;
+	coords[0].pc2_2 = b01 / coords[0].cdelt2;
+
+	coords[0].polyterms[0][0] = a20 / SQ(coords[0].cdelt1);
+	coords[0].polyterms[1][0] = a11 / (coords[0].cdelt1*coords[0].cdelt2);
+	coords[0].polyterms[2][0] = a02 / SQ(coords[0].cdelt2);
+
+	coords[0].polyterms[0][1] = b20 / SQ(coords[0].cdelt1);
+	coords[0].polyterms[1][1] = b11 / (coords[0].cdelt1*coords[0].cdelt2);
+	coords[0].polyterms[2][1] = b02 / SQ(coords[0].cdelt2);
+	for (i = 3; i < 7; i++) {
+	    coords[0].polyterms[i][0] = coords[0].polyterms[i][1] = 0.0;
+	}
+	break;
+      
+      case 3:
+	a10 = xcoeff[1][0] + 2*xcoeff[2][0]*Xo +   xcoeff[1][1]*Yo + 3*xcoeff[3][0]*Xo*Xo + 2*xcoeff[2][1]*Xo*Yo + xcoeff[1][2]*Yo*Yo;
+	a01 = xcoeff[0][1] + 2*xcoeff[0][2]*Yo +   xcoeff[1][1]*Xo + 3*xcoeff[0][3]*Yo*Yo + 2*xcoeff[1][2]*Xo*Yo + xcoeff[2][1]*Xo*Xo;
+	a20 = xcoeff[2][0] + 3*xcoeff[3][0]*Xo +   xcoeff[2][1]*Yo;
+	a11 = xcoeff[1][1] + 2*xcoeff[2][1]*Xo + 2*xcoeff[1][2]*Yo;
+	a02 = xcoeff[0][2] + 3*xcoeff[0][3]*Yo +   xcoeff[1][2]*Xo;
+	a30 = xcoeff[3][0];
+	a21 = xcoeff[2][1];
+	a12 = xcoeff[1][2];
+	a03 = xcoeff[0][3];
+
+	b10 = ycoeff[1][0] + 2*ycoeff[2][0]*Xo +   ycoeff[1][1]*Yo + 3*ycoeff[3][0]*Xo*Xo + 2*ycoeff[2][1]*Xo*Yo + ycoeff[1][2]*Yo*Yo;
+	b01 = ycoeff[0][1] + 2*ycoeff[0][2]*Yo +   ycoeff[1][1]*Xo + 3*ycoeff[0][3]*Yo*Yo + 2*ycoeff[1][2]*Xo*Yo + ycoeff[2][1]*Xo*Xo;
+	b20 = ycoeff[2][0] + 3*ycoeff[3][0]*Xo +   ycoeff[2][1]*Yo;
+	b11 = ycoeff[1][1] + 2*ycoeff[2][1]*Xo + 2*ycoeff[1][2]*Yo;
+	b02 = ycoeff[0][2] + 3*ycoeff[0][3]*Yo +   ycoeff[1][2]*Xo;
+	b30 = ycoeff[3][0];
+	b21 = ycoeff[2][1];
+	b12 = ycoeff[1][2];
+	b03 = ycoeff[0][3];
+
+	coords[0].pc1_1 = a10 / coords[0].cdelt1;
+	coords[0].pc1_2 = a01 / coords[0].cdelt2;
+	coords[0].pc2_1 = b10 / coords[0].cdelt1;
+	coords[0].pc2_2 = b01 / coords[0].cdelt2;
+
+	coords[0].polyterms[0][0] = a20 / SQ(coords[0].cdelt1);
+	coords[0].polyterms[1][0] = a11 / (coords[0].cdelt1*coords[0].cdelt2);
+	coords[0].polyterms[2][0] = a02 / SQ(coords[0].cdelt2);
+
+	coords[0].polyterms[3][0] = a30 / (SQ(coords[0].cdelt1)*coords[0].cdelt1);
+	coords[0].polyterms[4][0] = a21 / (SQ(coords[0].cdelt1)*coords[0].cdelt2);
+	coords[0].polyterms[5][0] = a12 / (SQ(coords[0].cdelt2)*coords[0].cdelt1);
+	coords[0].polyterms[6][0] = a03 / (SQ(coords[0].cdelt2)*coords[0].cdelt2);
+
+	coords[0].polyterms[0][1] = b20 / SQ(coords[0].cdelt1);
+	coords[0].polyterms[1][1] = b11 / (coords[0].cdelt1*coords[0].cdelt2);
+	coords[0].polyterms[2][1] = b02 / SQ(coords[0].cdelt2);
+
+	coords[0].polyterms[3][1] = b30 / (SQ(coords[0].cdelt1)*coords[0].cdelt1);
+	coords[0].polyterms[4][1] = b21 / (SQ(coords[0].cdelt1)*coords[0].cdelt2);
+	coords[0].polyterms[5][1] = b12 / (SQ(coords[0].cdelt2)*coords[0].cdelt1);
+	coords[0].polyterms[6][1] = b03 / (SQ(coords[0].cdelt2)*coords[0].cdelt2);
+	break;
+
+      default:
+	fprintf (stderr, "error: invalid order %d\n", coords[0].Npolyterms);
+	exit (2);
+    }
+
+    while (coords[0].crval1 < 0) coords[0].crval1 += 360.0;
+    while (coords[0].crval1 > 360.0) coords[0].crval1 -= 360.0;
+
+
+# else 
 
     psMetadataAddF32 (header, PS_LIST_TAIL, "CRVAL1", 	PS_META_REPLACE, "", toSky->R*DEG_RAD);
@@ -266,4 +409,6 @@
     psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", toFPA->y->coeff[1][0]/plateScale);
     psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", toFPA->y->coeff[0][1]/plateScale);
+
+# endif
 
     // alternative representations use