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);