Index: trunk/psLib/src/astro/psCoord.c
===================================================================
--- trunk/psLib/src/astro/psCoord.c	(revision 41531)
+++ trunk/psLib/src/astro/psCoord.c	(revision 41896)
@@ -40,4 +40,10 @@
 #include "psMinimizePolyFit.h"
 
+# define TEST_SAVE_INVERSE_TRANSFORM 0
+
+# if (TEST_SAVE_INVERSE_TRANSFORM)
+# include "psBinaryOp.h"
+# endif
+
 # define ELIXIR_CODE 1
 
@@ -86,5 +92,6 @@
     }
 
-    psPlaneTransform *out = psPlaneTransformAlloc(1, 1);
+    // since the output polynomial is 1st order, a Chebyshev is not really useful
+    psPlaneTransform *out = psPlaneTransformAlloc(1, 1, PS_POLYNOMIAL_ORD); 
 
     psF64 r12 = 0.0;
@@ -191,5 +198,5 @@
 }
 
-psPlaneTransform* psPlaneTransformAlloc(int order1, int order2)
+psPlaneTransform* psPlaneTransformAlloc(int order1, int order2, psPolynomialType type)
 {
     PS_ASSERT_INT_NONNEGATIVE(order1, NULL);
@@ -198,6 +205,6 @@
     psPlaneTransform *pt = psAlloc(sizeof(psPlaneTransform));
 
-    pt->x = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, order1, order2);
-    pt->y = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, order1, order2);
+    pt->x = psPolynomial2DAlloc(type, order1, order2);
+    pt->y = psPolynomial2DAlloc(type, order1, order2);
 
     psMemSetDeallocator(pt, (psFreeFunc) planeTransformFree);
@@ -797,4 +804,10 @@
     }
 
+    // both polynomials in both input transforms must match type -- and for now be ORD
+    psAssert (trans1->x->type == PS_POLYNOMIAL_ORD, "fix for CHEB");
+    psAssert (trans1->y->type == PS_POLYNOMIAL_ORD, "fix for CHEB");
+    psAssert (trans2->x->type == PS_POLYNOMIAL_ORD, "fix for CHEB");
+    psAssert (trans2->y->type == PS_POLYNOMIAL_ORD, "fix for CHEB");
+
     //
     // Determine the size of the new psPlaneTransform.
@@ -814,5 +827,5 @@
     psPlaneTransform *myPT = NULL;
     if (out == NULL) {
-        myPT = psPlaneTransformAlloc(orderX, orderY);
+	myPT = psPlaneTransformAlloc(orderX, orderY, PS_POLYNOMIAL_ORD);
     } else {
         if ((out->x->nX == orderX) &&
@@ -834,5 +847,5 @@
         } else {
             psFree(out);
-            myPT = psPlaneTransformAlloc(orderX, orderY);
+            myPT = psPlaneTransformAlloc(orderX, orderY, PS_POLYNOMIAL_ORD);
         }
     }
@@ -1037,5 +1050,5 @@
     psPlaneTransform *myPT = NULL;
     if (out == NULL) {
-        myPT = psPlaneTransformAlloc(order, order);
+      myPT = psPlaneTransformAlloc(order, order, PS_POLYNOMIAL_CHEB);
     } else {
       // the user has supplied a model with a specific order : fit that order
@@ -1071,4 +1084,59 @@
     result &= psVectorFitPolynomial2D(myPT->x, NULL, 0, xOut, NULL, xIn, yIn);
     result &= psVectorFitPolynomial2D(myPT->y, NULL, 0, yOut, NULL, xIn, yIn);
+
+# if (TEST_SAVE_INVERSE_TRANSFORM)
+
+    psVector *xFit = psPolynomial2DEvalVector (myPT->x, xIn, yIn);
+    psVector *xRes = (psVector *) psBinaryOp (NULL, xOut, "-", xFit);
+
+    psVector *yFit = psPolynomial2DEvalVector (myPT->y, xIn, yIn);
+    psVector *yRes = (psVector *) psBinaryOp (NULL, yOut, "-", yFit);
+
+    static int nOut = 0;
+    char filename[1024];
+    snprintf (filename, 1024, "test.fit.%03d.ply", nOut); 
+    FILE *fp = fopen (filename, "w");
+    
+    fprintf (fp, " ---- xFit ---- \n");
+
+    for (int ix = 0; ix < myPT->x->nX + 1; ix++) {
+      for (int iy = 0; iy < myPT->x->nY + 1; iy++) {
+	fprintf (fp, "%18.12e ", myPT->x->coeff[ix][iy]);
+      }
+      fprintf (fp, "\n");
+    }
+    fprintf (fp, " ---- yFit ---- \n");
+
+    for (int ix = 0; ix < myPT->y->nX + 1; ix++) {
+      for (int iy = 0; iy < myPT->y->nY + 1; iy++) {
+	fprintf (fp, "%18.12e ", myPT->y->coeff[ix][iy]);
+      }
+      fprintf (fp, "\n");
+    }
+    fclose (fp);
+    
+    snprintf (filename, 1024, "test.fit.%03d.dat", nOut); nOut ++;
+    FILE *f1 = fopen (filename, "w");
+    for (int i = 0; i < xFit->n; i++) {
+      fprintf (f1, "%d : %f %f : %f %f : %f %f : %f %f\n", i,
+	       xIn->data.F64[i], yIn->data.F64[i],
+	       xOut->data.F64[i], yOut->data.F64[i],
+	       xFit->data.F64[i], yFit->data.F64[i],
+	       xRes->data.F64[i], yRes->data.F64[i]);
+    }
+    fclose (f1);
+
+    psStats *myStats = psStatsAlloc (PS_STAT_SAMPLE_STDEV);
+    psVectorStats (myStats, xRes, NULL, NULL, 0); float dX = myStats->sampleStdev;
+    psVectorStats (myStats, yRes, NULL, NULL, 0); float dY = myStats->sampleStdev;
+    fprintf (stderr, "xRes Sigma: %f  --  yRes Sigma %f\n", dX, dY);
+
+    psFree (myStats);
+    psFree (xFit);
+    psFree (yFit);
+    psFree (xRes);
+    psFree (yRes);
+
+# endif
 
     psFree(inCoord);
Index: trunk/psLib/src/astro/psCoord.h
===================================================================
--- trunk/psLib/src/astro/psCoord.h	(revision 41531)
+++ trunk/psLib/src/astro/psCoord.h	(revision 41896)
@@ -178,6 +178,7 @@
 
 psPlaneTransform* psPlaneTransformAlloc(
-    int order1,                        ///< The order of the x term in the transform.
-    int order2                         ///< The order of the y term in the transform.
+    int order1,			      ///< The order of the x term in the transform.
+    int order2,			      ///< The order of the y term in the transform.
+    psPolynomialType type	      ///< The polynomial type (ORD or CHEB) for this transform
 ) PS_ATTR_MALLOC;