Index: trunk/psLib/src/astro/psCoord.c
===================================================================
--- trunk/psLib/src/astro/psCoord.c	(revision 3977)
+++ trunk/psLib/src/astro/psCoord.c	(revision 3990)
@@ -10,6 +10,6 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.67 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-05-19 05:18:20 $
+*  @version $Revision: 1.68 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-05-19 23:57:36 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -89,5 +89,5 @@
       (Y2 * (1.0 / (F - ((C*E)/B))));
  
-XXX: Since thre is now a general psPlaneTransformInvert() function, we
+XXX: Since there is now a general psPlaneTransformInvert() function, we
 should rename this.
  *****************************************************************************/
@@ -122,4 +122,25 @@
     psPlaneTransform *out = psPlaneTransformAlloc(2, 2);
 
+    /* This is sample code from IfA.  It didn't work initially, and I did not
+       spend any time debugging it.
+     
+        psF64 a = transform->x->coeff[1][0];
+        psF64 b = transform->x->coeff[0][1];
+        psF64 c = transform->y->coeff[1][0];
+        psF64 d = transform->y->coeff[0][1];
+        psF64 e = transform->x->coeff[0][0];
+        psF64 f = transform->y->coeff[0][0];
+     
+        psF64 invDet = 1.0 / (a * d - b * c); // Inverse of the determinant
+     
+        // Not entirely sure why this works, but it appears to do so....................................!
+        out->x->coeff[1][0] = invDet * a;
+        out->x->coeff[0][1] = - invDet * b;
+        out->y->coeff[1][0] = - invDet * c;
+        out->y->coeff[0][1] = invDet * d;
+     
+        out->x->coeff[0][0] = - invDet * (d * e + c * f);
+        out->y->coeff[0][0] = - invDet * (b * e + a * f);
+    */
     out->x->coeff[0][0] = (-D + ((F*A)/C)) / (E - ((F*B)/C));
     out->x->coeff[1][0] = -(F/C) / (E - ((F*B)/C));
@@ -1032,4 +1053,6 @@
     psS32 numCoords = PS_MIN(source->n, dest->n);
     psS32 order = PS_MAX(trans->x->nX, trans->x->nY);
+    order = PS_MAX(order, trans->y->nX);
+    order = PS_MAX(order, trans->y->nY);
 
     //
@@ -1133,5 +1156,5 @@
 psPlaneTransform *psPlaneTransformInvert(psPlaneTransform *out,
         const psPlaneTransform *in,
-        psRegion region,
+        psRegion *region,
         int nSamples)
 {
@@ -1151,5 +1174,5 @@
         psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Input transformation must have same nX==nY.");
     }
-    psS32 order = PS_MAX(in->x->nX, in->x->nY);
+    psS32 order = in->x->nX;
 
     psPlaneTransform *myPT = NULL;
@@ -1197,5 +1220,7 @@
 
     //
-    // Initialize the grid.
+    // Initialize the grid.  Since we want the inverse of the transformation, the
+    // inCoords are written to the outData vector, and the outCoords are written
+    // to the inData vector.
     //
     psS32 cnt = 0;
@@ -1205,5 +1230,4 @@
             inCoord->x = region.x0 + ((psF32) xint) * ((region.x1 - region.x0) / ((psF32) nSamples));
             (void)psPlaneTransformApply(outCoord, in, inCoord);
-
             ((psPlane *) outData->data[cnt])->x = inCoord->x;
             ((psPlane *) outData->data[cnt])->y = inCoord->y;
@@ -1214,4 +1238,5 @@
         }
     }
+    // XXX: what values should be used here?
     bool rc = psPlaneTransformFit(myPT, inData, outData, 10, 100.0);