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

Timestamp:

Dec 9, 2006, 11:11:55 AM (19 years ago)

Author:

eugene

Message:

higher-order terms

File:

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

trunk/psastro/src/psastroWCS.c

@@ -284 +284 @@
 
         // 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 +301 @@
         }
     }
+
     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 +329 @@
     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