Index: trunk/psLib/src/math/psPolynomial.c
===================================================================
--- trunk/psLib/src/math/psPolynomial.c	(revision 5074)
+++ trunk/psLib/src/math/psPolynomial.c	(revision 5090)
@@ -7,6 +7,6 @@
 *  polynomials.  It also contains a Gaussian functions.
 *
-*  @version $Revision: 1.122 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-09-20 20:48:35 $
+*  @version $Revision: 1.123 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-09-22 02:47:16 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -95,5 +95,5 @@
     unsigned int x = 0;
 
-    for (x = 0; x < poly->nX; x++) {
+    for (x = 0; x < poly->COOL_2D_nX; x++) {
         psFree(poly->coeff[x]);
         psFree(poly->coeffErr[x]);
@@ -110,6 +110,6 @@
     unsigned int y = 0;
 
-    for (x = 0; x < poly->nX; x++) {
-        for (y = 0; y < poly->nY; y++) {
+    for (x = 0; x < poly->COOL_3D_nX; x++) {
+        for (y = 0; y < poly->COOL_3D_nY; y++) {
             psFree(poly->coeff[x][y]);
             psFree(poly->coeffErr[x][y]);
@@ -132,7 +132,7 @@
     unsigned int z = 0;
 
-    for (x = 0; x < poly->nX; x++) {
-        for (y = 0; y < poly->nY; y++) {
-            for (z = 0; z < poly->nZ; z++) {
+    for (x = 0; x < poly->COOL_4D_nX; x++) {
+        for (y = 0; y < poly->COOL_4D_nY; y++) {
+            for (z = 0; z < poly->COOL_4D_nZ; z++) {
                 psFree(poly->coeff[x][y][z]);
                 psFree(poly->coeffErr[x][y][z]);
@@ -170,5 +170,5 @@
     chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *));
     for (unsigned int i = 0; i < maxChebyPoly; i++) {
-        chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD);
+        chebPolys[i] = psPolynomial1DAlloc(i, PS_POLYNOMIAL_ORD);
     }
 
@@ -182,8 +182,8 @@
 
         for (psS32 i = 2; i < maxChebyPoly; i++) {
-            for (psS32 j = 0; j < chebPolys[i - 1]->n; j++) {
+            for (psS32 j = 0; j < chebPolys[i - 1]->COOL_1D_n; j++) {
                 chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j];
             }
-            for (psS32 j = 0; j < chebPolys[i - 2]->n; j++) {
+            for (psS32 j = 0; j < chebPolys[i - 2]->COOL_1D_n; j++) {
                 chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j];
             }
@@ -210,11 +210,11 @@
             "---- Calling ordPolynomial1DEval(%lf)\n", x);
     psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 4,
-            "Polynomial order is %u\n", poly->n);
-    for (loop_x = 0; loop_x < poly->n; loop_x++) {
+            "Polynomial order is %u\n", poly->COOL_1D_n);
+    for (loop_x = 0; loop_x < poly->COOL_1D_n+1; loop_x++) {
         psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 4,
                 "Polynomial coeff[%u] is %lf\n", loop_x, poly->coeff[loop_x]);
     }
 
-    for (loop_x = 0; loop_x < poly->n; loop_x++) {
+    for (loop_x = 0; loop_x < poly->COOL_1D_n+1; loop_x++) {
         if (poly->mask[loop_x] == 0) {
             // XXX: If you set the tracelevel to 10 here, and later set the tracelevel to
@@ -238,10 +238,11 @@
     PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
     // XXX: Create a macro for this in psConstants.h
-    if (poly->n < 1) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error: Chebyshev polynomial as order %u.", poly->n);
+    if (poly->COOL_1D_n < 1) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error: Chebyshev polynomial is order %u.", poly->COOL_1D_n);
         return(NAN);
     }
     psVector *d;
-    unsigned int n = poly->n;
+    // XXX: n should be nTerms here (for clarity).  Or get rid of the variable.
+    unsigned int n = 1 + poly->COOL_1D_n;
     unsigned int i;
     psF64 tmp = 0.0;
@@ -302,9 +303,9 @@
     psPolynomial1D **chebPolys = NULL;
 
-    n = poly->n;
+    n = 1 + poly->COOL_1D_n;
     chebPolys = createChebyshevPolys(n);
 
     tmp = 0.0;
-    for (i=0;i<poly->n;i++) {
+    for (i=0;i<(1 + poly->COOL_1D_n);i++) {
         tmp+= (poly->coeff[i] * psPolynomial1DEval(x, chebPolys[i]));
     }
@@ -328,7 +329,7 @@
     psF64 ySum = 1.0;
 
-    for (loop_x = 0; loop_x < poly->nX; loop_x++) {
+    for (loop_x = 0; loop_x < (1 + poly->COOL_2D_nX); loop_x++) {
         ySum = xSum;
-        for (loop_y = 0; loop_y < poly->nY; loop_y++) {
+        for (loop_y = 0; loop_y < (1 + poly->COOL_2D_nY); loop_y++) {
             if (poly->mask[loop_x][loop_y] == 0) {
                 polySum += ySum * poly->coeff[loop_x][loop_y];
@@ -359,12 +360,12 @@
     // Determine how many Chebyshev polynomials
     // are needed, then create them.
-    maxChebyPoly = poly->nX;
-    if (poly->nY > maxChebyPoly) {
-        maxChebyPoly = poly->nY;
+    maxChebyPoly = 1 + poly->COOL_2D_nX;
+    if ((1 + poly->COOL_2D_nY) > maxChebyPoly) {
+        maxChebyPoly = (1 + poly->COOL_2D_nY);
     }
     chebPolys = createChebyshevPolys(maxChebyPoly);
 
-    for (loop_x = 0; loop_x < poly->nX; loop_x++) {
-        for (loop_y = 0; loop_y < poly->nY; loop_y++) {
+    for (loop_x = 0; loop_x < (1 + poly->COOL_2D_nX); loop_x++) {
+        for (loop_y = 0; loop_y < (1 + poly->COOL_2D_nY); loop_y++) {
             if (poly->mask[loop_x][loop_y] == 0) {
                 polySum += poly->coeff[loop_x][loop_y] *
@@ -394,9 +395,9 @@
     psF64 zSum = 1.0;
 
-    for (loop_x = 0; loop_x < poly->nX; loop_x++) {
+    for (loop_x = 0; loop_x < (1 + poly->COOL_3D_nX); loop_x++) {
         ySum = xSum;
-        for (loop_y = 0; loop_y < poly->nY; loop_y++) {
+        for (loop_y = 0; loop_y < (1 + poly->COOL_3D_nY); loop_y++) {
             zSum = ySum;
-            for (loop_z = 0; loop_z < poly->nZ; loop_z++) {
+            for (loop_z = 0; loop_z < (1 + poly->COOL_3D_nZ); loop_z++) {
                 if (poly->mask[loop_x][loop_y][loop_z] == 0) {
                     polySum += zSum * poly->coeff[loop_x][loop_y][loop_z];
@@ -430,16 +431,16 @@
     // Determine how many Chebyshev polynomials
     // are needed, then create them.
-    maxChebyPoly = poly->nX;
-    if (poly->nY > maxChebyPoly) {
-        maxChebyPoly = poly->nY;
-    }
-    if (poly->nZ > maxChebyPoly) {
-        maxChebyPoly = poly->nZ;
-    }
-    chebPolys = createChebyshevPolys(maxChebyPoly);
-
-    for (loop_x = 0; loop_x < poly->nX; loop_x++) {
-        for (loop_y = 0; loop_y < poly->nY; loop_y++) {
-            for (loop_z = 0; loop_z < poly->nZ; loop_z++) {
+    maxChebyPoly = poly->COOL_3D_nX;
+    if (poly->COOL_3D_nY > maxChebyPoly) {
+        maxChebyPoly = poly->COOL_3D_nY;
+    }
+    if (poly->COOL_3D_nZ > maxChebyPoly) {
+        maxChebyPoly = poly->COOL_3D_nZ;
+    }
+    chebPolys = createChebyshevPolys(maxChebyPoly + 1);
+
+    for (loop_x = 0; loop_x < (1 + poly->COOL_3D_nX); loop_x++) {
+        for (loop_y = 0; loop_y < (1 + poly->COOL_3D_nY); loop_y++) {
+            for (loop_z = 0; loop_z < (1 + poly->COOL_3D_nZ); loop_z++) {
                 if (poly->mask[loop_x][loop_y][loop_z] == 0) {
                     polySum += poly->coeff[loop_x][loop_y][loop_z] *
@@ -475,11 +476,11 @@
     psF64 tSum = 1.0;
 
-    for (loop_x = 0; loop_x < poly->nX; loop_x++) {
+    for (loop_x = 0; loop_x < (1 + poly->COOL_4D_nX); loop_x++) {
         ySum = xSum;
-        for (loop_y = 0; loop_y < poly->nY; loop_y++) {
+        for (loop_y = 0; loop_y < (1 + poly->COOL_4D_nY); loop_y++) {
             zSum = ySum;
-            for (loop_z = 0; loop_z < poly->nZ; loop_z++) {
+            for (loop_z = 0; loop_z < (1 + poly->COOL_4D_nZ); loop_z++) {
                 tSum = zSum;
-                for (loop_t = 0; loop_t < poly->nT; loop_t++) {
+                for (loop_t = 0; loop_t < (1 + poly->COOL_4D_nT); loop_t++) {
                     if (poly->mask[loop_x][loop_y][loop_z][loop_t] == 0) {
                         polySum += tSum * poly->coeff[loop_x][loop_y][loop_z][loop_t];
@@ -518,20 +519,20 @@
     // Determine how many Chebyshev polynomials
     // are needed, then create them.
-    maxChebyPoly = poly->nX;
-    if (poly->nY > maxChebyPoly) {
-        maxChebyPoly = poly->nY;
-    }
-    if (poly->nZ > maxChebyPoly) {
-        maxChebyPoly = poly->nZ;
-    }
-    if (poly->nT > maxChebyPoly) {
-        maxChebyPoly = poly->nT;
+    maxChebyPoly = poly->COOL_4D_nX;
+    if (poly->COOL_4D_nY > maxChebyPoly) {
+        maxChebyPoly = poly->COOL_4D_nY;
+    }
+    if (poly->COOL_4D_nZ > maxChebyPoly) {
+        maxChebyPoly = poly->COOL_4D_nZ;
+    }
+    if (poly->COOL_4D_nT > maxChebyPoly) {
+        maxChebyPoly = poly->COOL_4D_nT;
     }
     chebPolys = createChebyshevPolys(maxChebyPoly);
 
-    for (loop_x = 0; loop_x < poly->nX; loop_x++) {
-        for (loop_y = 0; loop_y < poly->nY; loop_y++) {
-            for (loop_z = 0; loop_z < poly->nZ; loop_z++) {
-                for (loop_t = 0; loop_t < poly->nT; loop_t++) {
+    for (loop_x = 0; loop_x < (1 + poly->COOL_4D_nX); loop_x++) {
+        for (loop_y = 0; loop_y < (1 + poly->COOL_4D_nY); loop_y++) {
+            for (loop_z = 0; loop_z < (1 + poly->COOL_4D_nZ); loop_z++) {
+                for (loop_t = 0; loop_t < (1 + poly->COOL_4D_nZ); loop_t++) {
                     if (poly->mask[loop_x][loop_y][loop_z][loop_t] == 0) {
                         polySum += poly->coeff[loop_x][loop_y][loop_z][loop_t] *
@@ -610,4 +611,7 @@
 /*****************************************************************************
     This routine must allocate memory for the polynomial structures.
+    XXX: replaces nterms variables with nOrder.  Lots of potential for bugs
+         here.  Probably should create separately named private variables
+         in these functions.
  *****************************************************************************/
 psPolynomial1D* psPolynomial1DAlloc(unsigned int n,
@@ -623,9 +627,9 @@
 
     newPoly->type = type;
-    newPoly->n = n;
+    newPoly->COOL_1D_n = n;
     newPoly->coeff = psAlloc(n * sizeof(psF64));
     newPoly->coeffErr = psAlloc(n * sizeof(psF64));
-    newPoly->mask = (psMaskType *)psAlloc(n * sizeof(psMaskType));
-    for (i = 0; i < n; i++) {
+    newPoly->mask = (psMaskType *)psAlloc((n + 1) * sizeof(psMaskType));
+    for (i = 0; i < (n + 1); i++) {
         newPoly->coeff[i] = 0.0;
         newPoly->coeffErr[i] = 0.0;
@@ -651,17 +655,17 @@
 
     newPoly->type = type;
-    newPoly->nX = nX;
-    newPoly->nY = nY;
-
-    newPoly->coeff = psAlloc(nX * sizeof(psF64 *));
-    newPoly->coeffErr = psAlloc(nX * sizeof(psF64 *));
-    newPoly->mask = (psMaskType **)psAlloc(nX * sizeof(psMaskType *));
-    for (x = 0; x < nX; x++) {
-        newPoly->coeff[x] = psAlloc(nY * sizeof(psF64));
-        newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64));
-        newPoly->mask[x] = (psMaskType *)psAlloc(nY * sizeof(psMaskType));
-    }
-    for (x = 0; x < nX; x++) {
-        for (y = 0; y < nY; y++) {
+    newPoly->COOL_2D_nX = nX;
+    newPoly->COOL_2D_nY = nY;
+
+    newPoly->coeff = psAlloc((1 + nX) * sizeof(psF64 *));
+    newPoly->coeffErr = psAlloc((1 + nX) * sizeof(psF64 *));
+    newPoly->mask = (psMaskType **)psAlloc((1 + nX) * sizeof(psMaskType *));
+    for (x = 0; x < (1 + nX); x++) {
+        newPoly->coeff[x] = psAlloc((1 + nY) * sizeof(psF64));
+        newPoly->coeffErr[x] = psAlloc((1 + nY) * sizeof(psF64));
+        newPoly->mask[x] = (psMaskType *)psAlloc((1 + nY) * sizeof(psMaskType));
+    }
+    for (x = 0; x < (1 + nX); x++) {
+        for (y = 0; y < (1 + nY); y++) {
             newPoly->coeff[x][y] = 0.0;
             newPoly->coeffErr[x][y] = 0.0;
@@ -691,24 +695,24 @@
 
     newPoly->type = type;
-    newPoly->nX = nX;
-    newPoly->nY = nY;
-    newPoly->nZ = nZ;
-
-    newPoly->coeff = psAlloc(nX * sizeof(psF64 **));
-    newPoly->coeffErr = psAlloc(nX * sizeof(psF64 **));
-    newPoly->mask = (psMaskType ***)psAlloc(nX * sizeof(psMaskType **));
-    for (x = 0; x < nX; x++) {
-        newPoly->coeff[x] = psAlloc(nY * sizeof(psF64 *));
-        newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64 *));
-        newPoly->mask[x] = (psMaskType **)psAlloc(nY * sizeof(psMaskType *));
-        for (y = 0; y < nY; y++) {
-            newPoly->coeff[x][y] = psAlloc(nZ * sizeof(psF64));
-            newPoly->coeffErr[x][y] = psAlloc(nZ * sizeof(psF64));
-            newPoly->mask[x][y] = (psMaskType *)psAlloc(nZ * sizeof(psMaskType));
-        }
-    }
-    for (x = 0; x < nX; x++) {
-        for (y = 0; y < nY; y++) {
-            for (z = 0; z < nZ; z++) {
+    newPoly->COOL_3D_nX = nX;
+    newPoly->COOL_3D_nY = nY;
+    newPoly->COOL_3D_nZ = nZ;
+
+    newPoly->coeff = psAlloc((nX + 1) * sizeof(psF64 **));
+    newPoly->coeffErr = psAlloc((nX + 1) * sizeof(psF64 **));
+    newPoly->mask = (psMaskType ***)psAlloc((nX + 1) * sizeof(psMaskType **));
+    for (x = 0; x < (1 + nX); x++) {
+        newPoly->coeff[x] = psAlloc((nY + 1) * sizeof(psF64 *));
+        newPoly->coeffErr[x] = psAlloc((nY + 1) * sizeof(psF64 *));
+        newPoly->mask[x] = (psMaskType **)psAlloc((nY + 1) * sizeof(psMaskType *));
+        for (y = 0; y < (nY + 1); y++) {
+            newPoly->coeff[x][y] = psAlloc((nZ + 1) * sizeof(psF64));
+            newPoly->coeffErr[x][y] = psAlloc((nZ + 1) * sizeof(psF64));
+            newPoly->mask[x][y] = (psMaskType *)psAlloc((nZ + 1) * sizeof(psMaskType));
+        }
+    }
+    for (x = 0; x < (nX + 1); x++) {
+        for (y = 0; y < (nY + 1); y++) {
+            for (z = 0; z < (nZ + 1); z++) {
                 newPoly->coeff[x][y][z] = 0.0;
                 newPoly->coeffErr[x][y][z] = 0.0;
@@ -742,31 +746,31 @@
 
     newPoly->type = type;
-    newPoly->nX = nX;
-    newPoly->nY = nY;
-    newPoly->nZ = nZ;
-    newPoly->nT = nT;
-
-    newPoly->coeff = psAlloc(nX * sizeof(psF64 ***));
-    newPoly->coeffErr = psAlloc(nX * sizeof(psF64 ***));
-    newPoly->mask = (psMaskType ****)psAlloc(nX * sizeof(psMaskType ***));
-    for (x = 0; x < nX; x++) {
-        newPoly->coeff[x] = psAlloc(nY * sizeof(psF64 **));
-        newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64 **));
-        newPoly->mask[x] = (psMaskType ***)psAlloc(nY * sizeof(psMaskType **));
-        for (y = 0; y < nY; y++) {
-            newPoly->coeff[x][y] = psAlloc(nZ * sizeof(psF64 *));
-            newPoly->coeffErr[x][y] = psAlloc(nZ * sizeof(psF64 *));
-            newPoly->mask[x][y] = (psMaskType **)psAlloc(nZ * sizeof(psMaskType *));
-            for (z = 0; z < nZ; z++) {
-                newPoly->coeff[x][y][z] = psAlloc(nT * sizeof(psF64));
-                newPoly->coeffErr[x][y][z] = psAlloc(nT * sizeof(psF64));
-                newPoly->mask[x][y][z] = (psMaskType *)psAlloc(nT * sizeof(psMaskType));
-            }
-        }
-    }
-    for (x = 0; x < nX; x++) {
-        for (y = 0; y < nY; y++) {
-            for (z = 0; z < nZ; z++) {
-                for (t = 0; t < nT; t++) {
+    newPoly->COOL_4D_nX = nX;
+    newPoly->COOL_4D_nY = nY;
+    newPoly->COOL_4D_nZ = nZ;
+    newPoly->COOL_4D_nT = nT;
+
+    newPoly->coeff = psAlloc((nX + 1) * sizeof(psF64 ***));
+    newPoly->coeffErr = psAlloc((nX + 1) * sizeof(psF64 ***));
+    newPoly->mask = (psMaskType ****)psAlloc((nX + 1) * sizeof(psMaskType ***));
+    for (x = 0; x < (nX + 1); x++) {
+        newPoly->coeff[x] = psAlloc((nY + 1) * sizeof(psF64 **));
+        newPoly->coeffErr[x] = psAlloc((nY + 1) * sizeof(psF64 **));
+        newPoly->mask[x] = (psMaskType ***)psAlloc((nY + 1) * sizeof(psMaskType **));
+        for (y = 0; y < (nY + 1); y++) {
+            newPoly->coeff[x][y] = psAlloc((nZ + 1) * sizeof(psF64 *));
+            newPoly->coeffErr[x][y] = psAlloc((nZ + 1) * sizeof(psF64 *));
+            newPoly->mask[x][y] = (psMaskType **)psAlloc((nZ + 1) * sizeof(psMaskType *));
+            for (z = 0; z < (nZ + 1); z++) {
+                newPoly->coeff[x][y][z] = psAlloc((nT + 1) * sizeof(psF64));
+                newPoly->coeffErr[x][y][z] = psAlloc((nT + 1) * sizeof(psF64));
+                newPoly->mask[x][y][z] = (psMaskType *)psAlloc((nT + 1) * sizeof(psMaskType));
+            }
+        }
+    }
+    for (x = 0; x < (nX + 1); x++) {
+        for (y = 0; y < (nY + 1); y++) {
+            for (z = 0; z < (nZ + 1); z++) {
+                for (t = 0; t < (nT + 1); t++) {
                     newPoly->coeff[x][y][z][t] = 0.0;
                     newPoly->coeffErr[x][y][z][t] = 0.0;
@@ -972,5 +976,5 @@
     }
 
-    // Allocate output vector
+    // Allocoutput vector
     tmp = psVectorAlloc(vecLen, PS_TYPE_F64);