Index: /trunk/psastro/src/psastroWCS.c
===================================================================
--- /trunk/psastro/src/psastroWCS.c	(revision 10599)
+++ /trunk/psastro/src/psastroWCS.c	(revision 10600)
@@ -284,4 +284,5 @@
 
 	// XXX this loop is rather arbitrary in length...
+	// XXX measure the error and use as criterion
 	for (int i = 0; i < 10; i++) {
 	    // NOTE: order is: [y][x]
@@ -300,94 +301,27 @@
 	}
     }
+
     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);
+    psPolynomial2D *xWCS = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, toFPA->x->nX, toFPA->x->nY);
+    psPolynomial2D *yWCS = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, toFPA->y->nX, toFPA->y->nY);
+
+    psPolynomial2D *xPx = psPolynomial2DCopy (toFPA->x);
+    psPolynomial2D *yPx = psPolynomial2DCopy (toFPA->y);
+
+    // skip the zero order terms
+    // XXX double check that these relationships are correct
+    for (int i = 0; i < toFPA->x->nX; i++) {
+	psPolynomial2D *xPy = psPolynomial2DCopy (xPx);
+	psPolynomial2D *yPy = psPolynomial2DCopy (yPx);
+	for (int j = 0; j < toFPA->x->nY; j++) {
+	    xWCS->coords[i][j] = psPolynomial2DEval (xPy, Xo, Yo) / (i*j) / pow(cdelt1, i) / pow(cdelt2, j);
+	    yWCS->coords[i][j] = psPolynomial2DEval (yPy, Xo, Yo) / (i*j) / pow(cdelt1, i) / pow(cdelt2, j);
+	    psPolynomial2D_dY(xPy, xPy);
+	    psPolynomial2D_dY(yPy, yPy);
+	}
+	psPolynomial2D_dX(xPx, xPx);
+	psPolynomial2D_dX(yPx, yPx);
     }
 
@@ -395,5 +329,32 @@
     while (coords[0].crval1 > 360.0) coords[0].crval1 -= 360.0;
 
-
+    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, "", Xo);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CRPIX2", 	PS_META_REPLACE, "", Yo);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT1", 	PS_META_REPLACE, "", cdelt1);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "CDELT2", 	PS_META_REPLACE, "", cdelt2);
+
+    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001001", PS_META_REPLACE, "", xWCS->coeff[1][0]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "PC001002", PS_META_REPLACE, "", xWCS->coeff[0][1]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002001", PS_META_REPLACE, "", yWCS->coeff[1][0]);
+    psMetadataAddF32 (header, PS_LIST_TAIL, "PC002002", PS_META_REPLACE, "", yWCS->coeff[0][1]);
+
+    // XXX respect the masks
+    for (int i = 0; i < xWCS->nX; i++) {
+	for (int j = 0; j < xWCS->nX; j++) {
+	    if (i + j < 2) continue;
+	    sprintf (name, "PCA1dX%1dY%1d", i, j);
+	    psMetadataAddF32 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", xWCS->coeff[i][j]);
+	}
+    }
+    for (int i = 0; i < yWCS->nX; i++) {
+	for (int j = 0; j < yWCS->nX; j++) {
+	    if (i + j < 2) continue;
+	    sprintf (name, "PCA2dX%1dY%1d", i, j);
+	    psMetadataAddF32 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", yWCS->coeff[i][j]);
+	}
+    }
+	    
 # else