Index: trunk/psLib/src/astro/psCoord.c
===================================================================
--- trunk/psLib/src/astro/psCoord.c	(revision 6383)
+++ trunk/psLib/src/astro/psCoord.c	(revision 6384)
@@ -10,6 +10,6 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.108 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-02-08 21:34:53 $
+*  @version $Revision: 1.109 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-02-08 21:52:05 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -58,12 +58,9 @@
 "transform" is linear.
  
-XXX: Since there is now a general psPlaneTransformInvert() function, we
-should rename this.
+XXX: This code no longer makes sense.  The merge must be reviewed.
  
-XXX: This code no longer makes sense.  The merge must be reviewed.
+XXX: below is the code using the standard matrix representation.  note that
+this inversion requires x->nX == 1, y->nY == 1 and x->nY <= 1, y->nX <= 1
  *****************************************************************************/
-// XXX EAM : below is the code using the standard matrix representation.
-//           note that this inversion requires x->nX == 1, y->nY == 1 and
-//           x->nY <= 1, y->nX <= 1
 psPlaneTransform *p_psPlaneTransformLinearInvert(psPlaneTransform *transform)
 {
@@ -177,5 +174,4 @@
 }
 
-// XXX: Must rewrite code and tests to use these functions.
 psPlane* psPlaneAlloc(void)
 {
@@ -326,23 +322,4 @@
 }
 
-
-/******************************************************************************
-XXX: Private Function.
- 
-XXX: Optimize this.  We don't need a while loop here.
- 
-piNormalize(): take an input angle in radians and convert it to the range 0:2*PI.
- *****************************************************************************/
-psF32 piNormalize(psF32 angle)
-{
-    while (angle < FLT_EPSILON) {
-        angle+=M_PI*2;
-    }
-
-    while (angle >= (M_PI*2)) {
-        angle-=M_PI*2;
-    }
-    return(angle);
-}
 
 void projectionFree(psProjection *p)
@@ -455,6 +432,7 @@
 }
 
-psPlane* psProject(const psSphere* coord,
-                   const psProjection* projection)
+psPlane* psProject(
+    const psSphere* coord,
+    const psProjection* projection)
 {
     return(p_psProject(NULL, coord, projection));
@@ -648,4 +626,5 @@
     //
     // XXX: rename, or verify, or recode, after the poly norder/nterm change.
+
     psPlaneTransform *myPT = NULL;
     if (out == NULL) {
@@ -791,5 +770,4 @@
  
 XXX: What about nRejIter?  Iterations?
-XXX: Use static vectors for internal data.
 XXX: This code has problems with data that corresponds to a non-linear fit.
  *****************************************************************************/
@@ -805,4 +783,5 @@
     PS_ASSERT_PTR_NON_NULL(dest, NULL);
 
+    psBool rc = true;
     psS32 numCoords = PS_MIN(source->n, dest->n);
     psS32 order = PS_MAX(trans->x->nX, trans->x->nY);
@@ -874,18 +853,28 @@
     // Solution via LU Decomposition
     //
-    // XXX: Check return codes
-    //
+    psVector *xSolution;
+    psVector *ySolution;
     psVector *permutation = psVectorAlloc(nCoeff, PS_TYPE_F64); // Permutation vector for LU Decomposition
     psImage *luMatrix = psMatrixLUD(NULL, &permutation, matrix); // LU decomposed matrix
-    psVector *xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
-    psVector *ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
-
-    //
-    // Stuff coefficients into transformation
-    //
-    for (psS32 i = 0, ijIndex = 0; i < order; i++) {
-        for (psS32 j = 0; j < order - i; j++, ijIndex++) {
-            trans->x->coeff[i][j] = xSolution->data.F64[ijIndex];
-            trans->y->coeff[i][j] = ySolution->data.F64[ijIndex];
+    if (luMatrix == NULL) {
+        rc = false;
+        psError(PS_ERR_UNKNOWN, true, "psMatrixLUD() returned NULL.  Returning FALSE.\n");
+    } else {
+        xSolution = psMatrixLUSolve(NULL, luMatrix, xVector, permutation); // Solution in x
+        ySolution = psMatrixLUSolve(NULL, luMatrix, yVector, permutation); // Solution in y
+
+        if ((xSolution == NULL) || (ySolution == NULL)) {
+            rc = false;
+            psError(PS_ERR_UNKNOWN, true, "psMatrixLUSolve() returned NULL.  Returning FALSE.\n");
+        } else {
+            //
+            // Stuff coefficients into transformation
+            //
+            for (psS32 i = 0, ijIndex = 0; i < order; i++) {
+                for (psS32 j = 0; j < order - i; j++, ijIndex++) {
+                    trans->x->coeff[i][j] = xSolution->data.F64[ijIndex];
+                    trans->y->coeff[i][j] = ySolution->data.F64[ijIndex];
+                }
+            }
         }
     }
@@ -900,5 +889,5 @@
     psFree(yVector);
 
-    return(true);
+    return(rc);
 }
 
@@ -907,10 +896,10 @@
 psPlaneTransformInvert(out, in, region, nSamples)
  
-// XXX: Use static data structures.
  *****************************************************************************/
-psPlaneTransform *psPlaneTransformInvert(psPlaneTransform *out,
-        const psPlaneTransform *in,
-        psRegion region,
-        int nSamples)
+psPlaneTransform *psPlaneTransformInvert(
+    psPlaneTransform *out,
+    const psPlaneTransform *in,
+    psRegion region,
+    int nSamples)
 {
     PS_ASSERT_PTR_NON_NULL(in, NULL);
@@ -1015,8 +1004,8 @@
     if (rc == true) {
         return(myPT);
-    }
-
-    // XXX: Generate an error message, or warning message.
-    return(NULL);
+    } else {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: could not fit a transform to coordinates.");
+        return(NULL);
+    }
 }