Index: trunk/psModules/src/astrom/pmAstrometryUtils.c
===================================================================
--- trunk/psModules/src/astrom/pmAstrometryUtils.c	(revision 24096)
+++ trunk/psModules/src/astrom/pmAstrometryUtils.c	(revision 35561)
@@ -23,4 +23,7 @@
 #include "pmAstrometryUtils.h"
 
+// this is used by the test output block
+static int Nout = 0;
+
 // given a 2D transformation -- L(x,y),M(x,y) -- find the coordinates x,y
 // for which L,M = 0,0. tol is the allowed error on x,y.
@@ -33,4 +36,20 @@
     double Xo = (trans->y->coeff[0][0]*trans->x->coeff[0][1] - trans->x->coeff[0][0]*trans->y->coeff[0][1])/R;
     double Yo = (trans->x->coeff[0][0]*trans->y->coeff[1][0] - trans->y->coeff[0][0]*trans->x->coeff[1][0])/R;
+
+    if (0) {
+      // this is a test output block, not used in normal operations
+      char filename[64];
+      snprintf (filename, 64, "trans.%03d.md", Nout);
+      FILE *f = fopen (filename, "w");
+
+      psMetadata *md = psMetadataAlloc();
+
+      psPolynomial2DtoMetadata (md, trans->x, "X");
+      psPolynomial2DtoMetadata (md, trans->y, "Y");
+      psMetadataConfigPrint (f, md);
+      psFree (md);
+      fclose (f);
+      Nout ++;
+    }
 
     // iterate to actual solution: requires small non-linear terms
@@ -52,11 +71,13 @@
         for (int i = 0; (dPos > tol) && (i < 20); i++) {
             // NOTE: order for Alpha is: [y][x]
+            // original Alpha->data.F32[1][0] = psPolynomial2DEval (XdY, Xo, Yo);
+            // original Alpha->data.F32[0][1] = psPolynomial2DEval (YdX, Xo, Yo);
             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[0][1] = psPolynomial2DEval (XdY, Xo, Yo);
+            Alpha->data.F32[1][0] = psPolynomial2DEval (YdX, Xo, Yo);
             Alpha->data.F32[1][1] = psPolynomial2DEval (YdY, Xo, Yo);
 
-            Beta->data.F32[0] = psPolynomial2DEval (trans->x, Xo, Yo);
-            Beta->data.F32[1] = psPolynomial2DEval (trans->y, Xo, Yo);
+            Beta->data.F32[0] = -1.0 * psPolynomial2DEval (trans->x, Xo, Yo);
+            Beta->data.F32[1] = -1.0 * psPolynomial2DEval (trans->y, Xo, Yo);
 
             if (!psMatrixGJSolve (Alpha, Beta)) {
@@ -71,8 +92,8 @@
 	    }
 
-            Xo -= Beta->data.F32[0];
-            Yo -= Beta->data.F32[1];
+            Xo += Beta->data.F32[0];
+            Yo += Beta->data.F32[1];
             dPos = hypot(Beta->data.F32[0], Beta->data.F32[1]);
-
+	    psTrace ("psastro", 5, "i: %d, Xo,Yo : %f %f  dX, dY: %f %f dPos: %f\n", i, Xo, Yo, Beta->data.F32[0], Beta->data.F32[1], dPos);
         }
         psFree (Alpha);