Index: trunk/psModules/src/astrom/pmAstrometryWCS.c
===================================================================
--- trunk/psModules/src/astrom/pmAstrometryWCS.c	(revision 10775)
+++ trunk/psModules/src/astrom/pmAstrometryWCS.c	(revision 10781)
@@ -7,6 +7,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-12-16 05:34:54 $
+ *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-12-17 09:46:09 $
  *
  *  Copyright 2006 Institute for Astronomy, University of Hawaii
@@ -312,12 +312,7 @@
 // interpret header WCS (only handles traditional WCS for the moment)
 // plateScale is nominal physical scale on tangent plane (radians / TPA physical units)
-bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, const psMetadata *header, double plateScale, bool isMosaic)
+bool pmAstromWCStoFPA (pmFPA *fpa, pmChip *chip, const pmAstromWCS *wcs, double plateScale, bool isMosaic)
 {
     psPlaneTransform *toFPA;
-
-    pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
-    if (!wcs) {
-        return false;
-    }
 
     /* at this point, we have extracted from the header the WCS terms in the form of a polynomial,
@@ -479,15 +474,13 @@
              chip->fromFPA->y->coeff[1][0], chip->fromFPA->y->coeff[0][1]);
 
-    psFree (wcs);
-
     return true;
 }
 
-// convert toFPA / toSky components to pmAstromWCS, then write to the headers
-// plateScale is nominal physical scale on tangent plane (microns / arcsecond)
-// this requires toTP to be the identity transformation
-bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip)
-{
-    // techinically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
+// convert toFPA / toSky components to pmAstromWCS
+// tolerance is in pixels
+pmAstromWCS *pmAstromWCSfromFPA (const pmFPA *fpa, const pmChip *chip, float tol)
+{
+
+    // technically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
     if (!psPlaneDistortIsDiagonal (fpa->toTPA))
         psAbort ("psastro", "invalid TPA transformation");
@@ -508,5 +501,5 @@
     // crpix1,2 = X,Y(crval1,2)
     // start with linear solution for Xo,Yo:
-    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 R  = (chip->toFPA->x->coeff[1][0]*chip->toFPA->y->coeff[0][1] - chip->toFPA->x->coeff[0][1]*chip->toFPA->y->coeff[1][0]);
     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;
@@ -525,5 +518,8 @@
         // XXX this loop is rather arbitrary in length...
         // XXX measure the error and use as criterion
-        for (int i = 0; i < 10; i++) {
+        /* this is the Newton-Raphson method to solve for Xo,Yo - it needs the high order terms to be small */
+        // Xo,Yo are in pixels;
+        float dPos = tol + 1;
+        for (int i = 0; (dPos > tol) && (i < 10); i++) {
             // NOTE: order is: [y][x]
             Alpha->data.F32[0][0] = psPolynomial2DEval (XdX, Xo, Yo);
@@ -537,6 +533,7 @@
             psMatrixGJSolveF32 (Alpha, Beta);
 
-            Xo += Beta->data.F32[0];
-            Yo += Beta->data.F32[1];
+            Xo -= Beta->data.F32[0];
+            Yo -= Beta->data.F32[1];
+            dPos = hypot(Xo,Yo);
         }
         psFree (Alpha);
@@ -595,5 +592,36 @@
         }
     }
-
+    return wcs;
+}
+
+// interpret header WCS (only handles traditional WCS for the moment)
+// plateScale is nominal physical scale on tangent plane (radians / TPA physical units)
+bool pmAstromReadWCS (pmFPA *fpa, pmChip *chip, const psMetadata *header, double plateScale, bool isMosaic)
+{
+    pmAstromWCS *wcs = pmAstromWCSfromHeader (header);
+    if (!wcs) {
+        return false;
+    }
+
+    bool status = pmAstromWCStoFPA (fpa, chip, wcs, plateScale, isMosaic);
+
+    psFree (wcs);
+    return status;
+}
+
+// convert toFPA / toSky components to pmAstromWCS, then write to the headers
+// tolerance is in pixels
+bool pmAstromWriteWCS (psMetadata *header, const pmFPA *fpa, const pmChip *chip, float tol)
+{
+    // XXX require chip->toFPA->x->nX == chip->toFPA->x->nY
+    // XXX require chip->toFPA->y->nX == chip->toFPA->y->nY
+    // XXX require chip->toFPA->x->nX == chip->toFPA->y->nX
+    // XXX require chip->toFPA->nX == 1,2,3
+
+    // technically, we can have a plate scale here (fpa->toTPA:dx,dy != 1)
+    if (!psPlaneDistortIsDiagonal (fpa->toTPA))
+        psAbort ("psastro", "invalid TPA transformation");
+
+    pmAstromWCS *wcs = pmAstromWCSfromFPA(fpa, chip, tol);
     pmAstromWCStoHeader (header, wcs);