Index: /trunk/psLib/src/astro/psCoord.c
===================================================================
--- /trunk/psLib/src/astro/psCoord.c	(revision 3714)
+++ /trunk/psLib/src/astro/psCoord.c	(revision 3715)
@@ -10,6 +10,6 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.64 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-04-19 04:16:02 $
+*  @version $Revision: 1.65 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-04-19 05:46:24 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -91,27 +91,4 @@
 XXX: Since thre is now a general psPlaneTransformInvert() function, we
 should rename this.
- 
-X1, Y1 = (2, 3)
- 
-    X2 = 1 + 4 + 9   = 14
-    Y2 = 4 + 10 + 18 = 32
- 
-in inverse
- 
-Y1 = (32 - ((5/2) * 14) - 4 + ((5)/2)) / (6 - ((15)/2));
-Y1 = (32 - 35 - 4 + (5/2)) / (6 - 15/2);
-Y1 = (-7 + (5/2)) / (12/2 - 15/2);
-Y1 = (-7 + (5/2)) / (-3/2);
-Y1 = (-14/2 + 5/2) / (-3/2);
-Y1 = (-9/2) / (-3/2);
-Y1 = (9/2) / (3/2);
-Y1 = 3
- 
-X1 = (Y2 - ((F/C) * X2) - D + ((F*A)/C)) / (E - ((F*B)/C));
-X1 = (32 - ((6/3) * 14) - 4 + ((6)/3)) / (5 - ((12)/3));
-X1 = (32 - (2 * 14) - 4 + 2 / (5 - 4);
-X1 = (32 - 28 - 4 + 2) / 1;
-X1 = 2;
- 
  *****************************************************************************/
 psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform)
@@ -128,5 +105,4 @@
     psF64 F = 0.0;
 
-    // XXX: Test this for correctness.
     A = transform->x->coeff[0][0];
     if (transform->x->nX >= 2) {
@@ -144,5 +120,4 @@
     }
 
-    // XXX: Use the constructor here.
     psPlaneTransform *out = psPlaneTransformAlloc(2, 2);
 
@@ -877,9 +852,9 @@
                 for (psS32 t2x = 0 ; t2x < trans2->nX ; t2x++) {
                     for (psS32 t2y = 0 ; t2y < trans2->nY ; t2y++) {
-                        /* VERIFY MACRO
+                        /* Possible debug-only macro which checks these coords?
                         if ((t1x+t2x) >= orderX)
-                            printf("CRAP 1\n");
+                            printf("BAD 1\n");
                         if ((t1y+t2y) >= orderY)
-                            printf("CRAP 2\n");
+                            printf("BAD 2\n");
                         */
                         out->coeff[t1x+t2x][t1y+t2y]+= (trans1->coeff[t1x][t1y] * trans2->coeff[t2x][t2y]);
@@ -1051,5 +1026,4 @@
 
     psS32 numCoords = PS_MIN(source->n, dest->n);
-    // This is not really necessary because of above conditionals.
     psS32 order = PS_MAX(trans->x->nX, trans->x->nY);
 
Index: /trunk/psLib/src/astronomy/psCoord.c
===================================================================
--- /trunk/psLib/src/astronomy/psCoord.c	(revision 3714)
+++ /trunk/psLib/src/astronomy/psCoord.c	(revision 3715)
@@ -10,6 +10,6 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.64 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-04-19 04:16:02 $
+*  @version $Revision: 1.65 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-04-19 05:46:24 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -91,27 +91,4 @@
 XXX: Since thre is now a general psPlaneTransformInvert() function, we
 should rename this.
- 
-X1, Y1 = (2, 3)
- 
-    X2 = 1 + 4 + 9   = 14
-    Y2 = 4 + 10 + 18 = 32
- 
-in inverse
- 
-Y1 = (32 - ((5/2) * 14) - 4 + ((5)/2)) / (6 - ((15)/2));
-Y1 = (32 - 35 - 4 + (5/2)) / (6 - 15/2);
-Y1 = (-7 + (5/2)) / (12/2 - 15/2);
-Y1 = (-7 + (5/2)) / (-3/2);
-Y1 = (-14/2 + 5/2) / (-3/2);
-Y1 = (-9/2) / (-3/2);
-Y1 = (9/2) / (3/2);
-Y1 = 3
- 
-X1 = (Y2 - ((F/C) * X2) - D + ((F*A)/C)) / (E - ((F*B)/C));
-X1 = (32 - ((6/3) * 14) - 4 + ((6)/3)) / (5 - ((12)/3));
-X1 = (32 - (2 * 14) - 4 + 2 / (5 - 4);
-X1 = (32 - 28 - 4 + 2) / 1;
-X1 = 2;
- 
  *****************************************************************************/
 psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform)
@@ -128,5 +105,4 @@
     psF64 F = 0.0;
 
-    // XXX: Test this for correctness.
     A = transform->x->coeff[0][0];
     if (transform->x->nX >= 2) {
@@ -144,5 +120,4 @@
     }
 
-    // XXX: Use the constructor here.
     psPlaneTransform *out = psPlaneTransformAlloc(2, 2);
 
@@ -877,9 +852,9 @@
                 for (psS32 t2x = 0 ; t2x < trans2->nX ; t2x++) {
                     for (psS32 t2y = 0 ; t2y < trans2->nY ; t2y++) {
-                        /* VERIFY MACRO
+                        /* Possible debug-only macro which checks these coords?
                         if ((t1x+t2x) >= orderX)
-                            printf("CRAP 1\n");
+                            printf("BAD 1\n");
                         if ((t1y+t2y) >= orderY)
-                            printf("CRAP 2\n");
+                            printf("BAD 2\n");
                         */
                         out->coeff[t1x+t2x][t1y+t2y]+= (trans1->coeff[t1x][t1y] * trans2->coeff[t2x][t2y]);
@@ -1051,5 +1026,4 @@
 
     psS32 numCoords = PS_MIN(source->n, dest->n);
-    // This is not really necessary because of above conditionals.
     psS32 order = PS_MAX(trans->x->nX, trans->x->nY);
 
Index: /trunk/psLib/test/astronomy/tst_psAstrometry01.c
===================================================================
--- /trunk/psLib/test/astronomy/tst_psAstrometry01.c	(revision 3714)
+++ /trunk/psLib/test/astronomy/tst_psAstrometry01.c	(revision 3715)
@@ -5,6 +5,6 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-04-19 04:16:02 $
+*  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-04-19 05:46:24 $
 *
 * XXX: Must test
@@ -41,4 +41,8 @@
                               {NULL}
                           };
+
+#define PS_PERCENT_COMPARE(X, Y, PERCENT_FRACTION) \
+(fabs((Y)-(X))/fabs(X) < (PERCENT_FRACTION))
+
 #define RA  1.0
 #define DEC  2.0
@@ -764,6 +768,14 @@
 }
 
-#define NUM_DATA 100
-
+
+#define TST05_X_MIN 0.0
+#define TST05_X_MAX 10.0
+#define TST05_Y_MIN 0.0
+#define TST05_Y_MAX 10.0
+#define TST05_NUM_X 10.0
+#define TST05_NUM_Y 10.0
+#define TST05_X_POLY_ORDER 2
+#define TST05_Y_POLY_ORDER 2
+#define TST05_NUM_DATA (TST05_NUM_X * TST05_NUM_Y)
 /******************************************************************************
 XXX: This is only a rudimentary test of the psPlaneTransformFit() function.
@@ -775,17 +787,65 @@
 {
     bool testStatus = true;
-    psPlaneTransform *trans = psPlaneTransformAlloc(2, 2);
-    psArray *src = psArrayAlloc(NUM_DATA);
-    psArray *dst = psArrayAlloc(NUM_DATA);
+    psPlaneTransform *trans = psPlaneTransformAlloc(TST05_X_POLY_ORDER, TST05_Y_POLY_ORDER);
+    psPlaneTransform *transInit = psPlaneTransformAlloc(TST05_X_POLY_ORDER, TST05_Y_POLY_ORDER);
+    psArray *src = psArrayAlloc((int) TST05_NUM_DATA);
+    psArray *dst = psArrayAlloc((int) TST05_NUM_DATA);
     bool rc;
 
-    for (psS32 i = 0 ; i < NUM_DATA ; i++) {
-        src->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
-        dst->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
-        ((psVector *) src->data[i])->data.F32[0] = (psF32) i;
-        ((psVector *) src->data[i])->data.F32[1] = (psF32) i;
-        ((psVector *) dst->data[i])->data.F32[0] = (psF32) 2+i;
-        ((psVector *) dst->data[i])->data.F32[1] = (psF32) 3+i;
-    }
+
+    //
+    // We set an arbitrary non-linear transformation.
+    //
+
+    for (psS32 x = 0 ; x < TST05_X_POLY_ORDER ; x++) {
+        for (psS32 y = 0 ; y < TST05_Y_POLY_ORDER ; y++) {
+            transInit->x->coeff[x][y] = ((psF32) (random() % 1000000)) / 100000.0;
+            transInit->y->coeff[x][y] = ((psF32) (random() % 1000000)) / 100000.0;
+        }
+    }
+
+    // okay, for now, it's linear.
+    transInit->x->coeff[1][1] = 0.0;
+    transInit->y->coeff[1][1] = 0.0;
+
+    /*
+        //
+        // We set an arbitrary linear transformation.
+        // 
+        transInit->x->coeff[0][0] = 1.0;
+        transInit->x->coeff[0][1] = 3.0;
+        transInit->x->coeff[1][0] = 2.0;
+        transInit->x->coeff[1][1] = 7.0;
+        transInit->y->coeff[0][0] = 4.0;
+        transInit->y->coeff[0][1] = 6.0;
+        transInit->y->coeff[1][0] = 5.0;
+        transInit->y->coeff[1][1] = 3.0;
+    */
+    //
+    // We generate a grid of input data points in the x1,y1 plane, calculate the
+    // corresponding valuesin the x2,y2 in the transformed plane, and set these to
+    // the src and dst psVectors.
+    //
+    psPlane *inData = psPlaneAlloc();
+    psPlane *outData = NULL;
+    psS32 i = 0;
+    for (inData->x = TST05_X_MIN ; inData->x < TST05_X_MAX ; inData->x+= (TST05_X_MAX-TST05_X_MIN)/TST05_NUM_X) {
+        for (inData->y = TST05_Y_MIN ; inData->y < TST05_Y_MAX ; inData->y+= (TST05_Y_MAX-TST05_Y_MIN)/TST05_NUM_Y) {
+            src->data[i] = (psPtr *) psPlaneAlloc();
+            dst->data[i] = (psPtr *) psPlaneAlloc();
+
+            ((psPlane *) src->data[i])->x = inData->x;
+            ((psPlane *) src->data[i])->y = inData->y;
+
+            outData = psPlaneTransformApply(NULL, transInit, inData);
+
+            ((psPlane *) dst->data[i])->x = outData->x;
+            ((psPlane *) dst->data[i])->y = outData->y;
+            psFree(outData);
+
+            i++;
+        }
+    }
+    psFree(inData);
 
     printf("----------------------------------------------------------------------------------\n");
@@ -817,8 +877,33 @@
     rc = psPlaneTransformFit(trans, src, dst, 100, 100.0);
     if (rc != true) {
-        printf("TEST ERROR: psPlaneTransformFit() returned TRUE.\n");
+        printf("TEST ERROR: psPlaneTransformFit() returned FALSE.\n");
         testStatus = false;
     }
 
+    printf("The initial transformation was:\n");
+    printTrans(transInit);
+    printf("The derived transformation is:\n");
+    printTrans(trans);
+
+    //
+    // For the initial grid of input points, we transform them to output points with
+    // the derived transformation, and verify that they are within 10%.
+    //
+
+    for (psS32 i = 0 ; i < src->n ; i++) {
+        psPlane *inData = (psPlane *) src->data[i];
+        psPlane *outData = (psPlane *) dst->data[i];
+        psPlane *outDataDeriv = psPlaneTransformApply(NULL, trans, inData);
+
+        if (!PS_PERCENT_COMPARE(outDataDeriv->x, outData->x, 0.10) ||
+                !PS_PERCENT_COMPARE(outDataDeriv->y, outData->y, 0.10)) {
+            printf("TEST ERROR: the derived output coords were (%f, %f) should have been (%f, %f).\n",
+                   outDataDeriv->x, outDataDeriv->y, outData->x, outData->y);
+            testStatus = false;
+        }
+        psFree(outDataDeriv);
+    }
+
+    psFree(transInit);
     psFree(trans);
     psFree(src);
@@ -886,9 +971,9 @@
         trans->x->coeff[0][1] = ((psF32) (random() % 1000000)) / 100000.0;
         trans->x->coeff[1][0] = ((psF32) (random() % 1000000)) / 100000.0;
-        trans->x->coeff[1][1] = 0.1;
+        trans->x->coeff[1][1] = 0.0;
         trans->y->coeff[0][0] = ((psF32) (random() % 1000000)) / 100000.0;
         trans->y->coeff[0][1] = ((psF32) (random() % 1000000)) / 100000.0;
         trans->y->coeff[1][0] = ((psF32) (random() % 1000000)) / 100000.0;
-        trans->y->coeff[1][1] = 0.1;
+        trans->y->coeff[1][1] = 0.0;
 
         transOut = psPlaneTransformInvert(NULL, trans, myRegion, 10);