Changeset 10775 for trunk/psModules/src/astrom/pmAstrometryWCS.c

Timestamp:

Dec 15, 2006, 7:34:54 PM (19 years ago)

Author:

magnier

Message:

polynomial fits

File:

• trunk/psModules/src/astrom/pmAstrometryWCS.c (modified) (9 diffs)

trunk/psModules/src/astrom/pmAstrometryWCS.c

@@ -7 +7 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-12-14 06:26:15 $
+ *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-12-16 05:34:54 $
  *
  *  Copyright 2006 Institute for Astronomy, University of Hawaii
@@ -190 +190 @@
             wcs->trans->x->coeff[0][1] = -wcs->cdelt2 * sin(rotate*PM_RAD_DEG); // == PC1_2
             wcs->trans->y->coeff[1][0] = +wcs->cdelt1 * sin(rotate*PM_RAD_DEG); // == PC2_1
-            wcs->trans->y->coeff[1][0] = +wcs->cdelt2 * cos(rotate*PM_RAD_DEG); // == PC2_2
+            wcs->trans->y->coeff[0][1] = +wcs->cdelt2 * cos(rotate*PM_RAD_DEG); // == PC2_2
             return wcs;
         }
@@ -207 +207 @@
                     if (i + j < 2)
                         continue;
-                    if (i + j > fitOrder)
+                    if (i + j > fitOrder) {
+                        wcs->trans->x->mask[i][j] = 1;
+                        wcs->trans->y->mask[i][j] = 1;
                         continue;
+                    }
                     sprintf (name, "PCA1X%1dY%1d", i, j);
                     wcs->trans->x->coeff[i][j] = pow(wcs->cdelt1, i) * pow(wcs->cdelt2, j) * psMetadataLookupF32 (&status, header, name);
@@ -487 +490 @@
 {
     // techinically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
-    if (!psPlaneDistortIsIdentity (fpa->toTPA))
+    if (!psPlaneDistortIsDiagonal (fpa->toTPA))
         psAbort ("psastro", "invalid TPA transformation");
 
@@ -505 +508 @@
     // crpix1,2 = X,Y(crval1,2)
     // start with linear solution for Xo,Yo:
-    double R  = (chip->toFPA->x->coeff[1][0]*chip->toFPA->x->coeff[0][1] - chip->toFPA->x->coeff[0][1]*chip->toFPA->x->coeff[1][0]);
-    double Xo = (chip->toFPA->x->coeff[0][0]*chip->toFPA->x->coeff[0][1] - chip->toFPA->x->coeff[0][0]*chip->toFPA->x->coeff[0][1])/R;
-    double Yo = (chip->toFPA->x->coeff[0][0]*chip->toFPA->x->coeff[1][0] - chip->toFPA->x->coeff[0][0]*chip->toFPA->x->coeff[1][0])/R;
+    double R  = (chip->toFPA->x->coeff[1][0]*chip->toFPA->y->coeff[0][1] - chip->toFPA->x->coeff[0][1]*chip->toFPA->y->coeff[0][1]);
+    double Xo = (chip->toFPA->y->coeff[0][0]*chip->toFPA->x->coeff[0][1] - chip->toFPA->x->coeff[0][0]*chip->toFPA->y->coeff[0][1])/R;
+    double Yo = (chip->toFPA->x->coeff[0][0]*chip->toFPA->y->coeff[1][0] - chip->toFPA->y->coeff[0][0]*chip->toFPA->x->coeff[1][0])/R;
 
     // iterate to actual solution: requires small non-linear terms
@@ -582 +585 @@
     double cdelt1 = fpa->toTPA->x->coeff[1][0][0][0]*fpa->toSky->Xs*PM_DEG_RAD;
     double cdelt2 = fpa->toTPA->y->coeff[0][1][0][0]*fpa->toSky->Ys*PM_DEG_RAD;
+    wcs->cdelt1 = cdelt1;
+    wcs->cdelt2 = cdelt2;
 
     // convert wcs->trans to a matrix which yields L,M in pixels
@@ -741 +746 @@
 }
 
-// check that the given psPlaneDistort is the identity
-bool psPlaneDistortIsIdentity (psPlaneDistort *distort)
+// check that the given psPlaneDistort is the identity * (Xs,Ys)
+bool psPlaneDistortIsDiagonal (psPlaneDistort *distort)
 {
 
@@ -787 +792 @@
             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 {
                 status &= !distort->x->mask[i][j][0][0];
-            } else {
-                if (i + j != 1) {
+                if ((i == 1) && (i + j == 1)) {
+                    // linear, diagonal terms must be 1.0
+                    status &= (fabs(distort->x->coeff[i][j][0][0]) > FLT_EPSILON);
+                } else {
                     // 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);
+                    status &= (fabs(distort->x->coeff[i][j][0][0]) < FLT_EPSILON);
                 }
             }
@@ -805 +811 @@
             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 {
                 status &= !distort->y->mask[i][j][0][0];
-            } else {
-                if (i + j != 1) {
+                if ((j == 1) && (i + j == 1)) {
+                    // linear, diagonal terms must be 1.0
+                    status &= (fabs(distort->y->coeff[i][j][0][0]) > FLT_EPSILON);
+                } else {
                     // 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);
+                    status &= (fabs(distort->y->coeff[i][j][0][0]) < FLT_EPSILON);
                 }
             }