Index: /trunk/psLib/src/math/psPolynomial.c
===================================================================
--- /trunk/psLib/src/math/psPolynomial.c	(revision 5065)
+++ /trunk/psLib/src/math/psPolynomial.c	(revision 5066)
@@ -7,6 +7,6 @@
 *  polynomials.  It also contains a Gaussian functions.
 *
-*  @version $Revision: 1.120 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-09-11 22:18:40 $
+*  @version $Revision: 1.121 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-09-19 19:53:13 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -162,5 +162,5 @@
 outer coefficients of the Chebyshev polynomials.
  *****************************************************************************/
-static psPolynomial1D **createChebyshevPolys(psS32 maxChebyPoly)
+static psPolynomial1D **createChebyshevPolys(unsigned int maxChebyPoly)
 {
     PS_ASSERT_INT_NONNEGATIVE(maxChebyPoly, NULL);
@@ -169,5 +169,5 @@
 
     chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *));
-    for (psS32 i = 0; i < maxChebyPoly; i++) {
+    for (unsigned int i = 0; i < maxChebyPoly; i++) {
         chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD);
     }
@@ -191,5 +191,5 @@
     } else {
         // XXX: Code this.
-        printf("WARNING: %d-order chebyshev polynomials not correctly implemented.\n", maxChebyPoly);
+        printf("WARNING: %u-order chebyshev polynomials not correctly implemented.\n", maxChebyPoly);
     }
 
@@ -200,7 +200,8 @@
     Polynomial coefficients will be accessed in [w][x][y][z] fashion.
  *****************************************************************************/
-static psF64 ordPolynomial1DEval(psF64 x, const psPolynomial1D* poly)
-{
-    psS32 loop_x = 0;
+static psF64 ordPolynomial1DEval(psF64 x,
+                                 const psPolynomial1D* poly)
+{
+    unsigned int loop_x = 0;
     psF64 polySum = 0.0;
     psF64 xSum = 1.0;
@@ -209,8 +210,8 @@
             "---- Calling ordPolynomial1DEval(%lf)\n", x);
     psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 4,
-            "Polynomial order is %d\n", poly->n);
+            "Polynomial order is %u\n", poly->n);
     for (loop_x = 0; loop_x < poly->n; loop_x++) {
         psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 4,
-                "Polynomial coeff[%d] is %lf\n", loop_x, poly->coeff[loop_x]);
+                "Polynomial coeff[%u] is %lf\n", loop_x, poly->coeff[loop_x]);
     }
 
@@ -230,15 +231,16 @@
 // XXX: How does the mask vector effect Crenshaw's formula?
 // XXX: We assume that x is scaled between -1.0 and 1.0;
-static psF64 chebPolynomial1DEval(psF64 x, const psPolynomial1D* poly)
+static psF64 chebPolynomial1DEval(psF64 x,
+                                  const psPolynomial1D* poly)
 {
     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 %d.", poly->n);
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error: Chebyshev polynomial as order %u.", poly->n);
         return(NAN);
     }
     psVector *d;
-    psS32 n = poly->n;
-    psS32 i;
+    unsigned int n = poly->n;
+    unsigned int i;
     psF64 tmp = 0.0;
 
@@ -318,6 +320,6 @@
     PS_ASSERT_POLY_NON_NULL(poly, NAN);
 
-    psS32 loop_x = 0;
-    psS32 loop_y = 0;
+    unsigned int loop_x = 0;
+    unsigned int loop_y = 0;
     psF64 polySum = 0.0;
     psF64 xSum = 1.0;
@@ -338,5 +340,7 @@
 }
 
-static psF64 chebPolynomial2DEval(psF64 x, psF64 y, const psPolynomial2D* poly)
+static psF64 chebPolynomial2DEval(psF64 x,
+                                  psF64 y,
+                                  const psPolynomial2D* poly)
 {
     PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
@@ -344,10 +348,10 @@
     PS_ASSERT_POLY_NON_NULL(poly, NAN);
 
-    psS32 loop_x = 0;
-    psS32 loop_y = 0;
-    psS32 i = 0;
+    unsigned int loop_x = 0;
+    unsigned int loop_y = 0;
+    unsigned int i = 0;
     psF64 polySum = 0.0;
     psPolynomial1D* *chebPolys = NULL;
-    psS32 maxChebyPoly = 0;
+    unsigned int maxChebyPoly = 0;
 
     // Determine how many Chebyshev polynomials
@@ -375,9 +379,12 @@
 }
 
-static psF64 ordPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psPolynomial3D* poly)
-{
-    psS32 loop_x = 0;
-    psS32 loop_y = 0;
-    psS32 loop_z = 0;
+static psF64 ordPolynomial3DEval(psF64 x,
+                                 psF64 y,
+                                 psF64 z,
+                                 const psPolynomial3D* poly)
+{
+    unsigned int loop_x = 0;
+    unsigned int loop_y = 0;
+    unsigned int loop_z = 0;
     psF64 polySum = 0.0;
     psF64 xSum = 1.0;
@@ -403,16 +410,19 @@
 }
 
-static psF64 chebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psPolynomial3D* poly)
+static psF64 chebPolynomial3DEval(psF64 x,
+                                  psF64 y,
+                                  psF64 z,
+                                  const psPolynomial3D* poly)
 {
     PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
     PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0);
     PS_ASSERT_DOUBLE_WITHIN_RANGE(z, -1.0, 1.0, 0.0);
-    psS32 loop_x = 0;
-    psS32 loop_y = 0;
-    psS32 loop_z = 0;
-    psS32 i = 0;
+    unsigned int loop_x = 0;
+    unsigned int loop_y = 0;
+    unsigned int loop_z = 0;
+    unsigned int i = 0;
     psF64 polySum = 0.0;
     psPolynomial1D* *chebPolys = NULL;
-    psS32 maxChebyPoly = 0;
+    unsigned int maxChebyPoly = 0;
 
     // Determine how many Chebyshev polynomials
@@ -447,10 +457,14 @@
 }
 
-static psF64 ordPolynomial4DEval(psF64 x, psF64 y, psF64 z, psF64 t, const psPolynomial4D* poly)
-{
-    psS32 loop_x = 0;
-    psS32 loop_y = 0;
-    psS32 loop_z = 0;
-    psS32 loop_t = 0;
+static psF64 ordPolynomial4DEval(psF64 x,
+                                 psF64 y,
+                                 psF64 z,
+                                 psF64 t,
+                                 const psPolynomial4D* poly)
+{
+    unsigned int loop_x = 0;
+    unsigned int loop_y = 0;
+    unsigned int loop_z = 0;
+    unsigned int loop_t = 0;
     psF64 polySum = 0.0;
     psF64 xSum = 1.0;
@@ -481,5 +495,9 @@
 }
 
-static psF64 chebPolynomial4DEval(psF64 x, psF64 y, psF64 z, psF64 t, const psPolynomial4D* poly)
+static psF64 chebPolynomial4DEval(psF64 x,
+                                  psF64 y,
+                                  psF64 z,
+                                  psF64 t,
+                                  const psPolynomial4D* poly)
 {
     PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
@@ -487,12 +505,12 @@
     PS_ASSERT_DOUBLE_WITHIN_RANGE(z, -1.0, 1.0, 0.0);
     PS_ASSERT_DOUBLE_WITHIN_RANGE(t, -1.0, 1.0, 0.0);
-    psS32 loop_x = 0;
-    psS32 loop_y = 0;
-    psS32 loop_z = 0;
-    psS32 loop_t = 0;
-    psS32 i = 0;
+    unsigned int loop_x = 0;
+    unsigned int loop_y = 0;
+    unsigned int loop_z = 0;
+    unsigned int loop_t = 0;
+    unsigned int i = 0;
     psF64 polySum = 0.0;
     psPolynomial1D* *chebPolys = NULL;
-    psS32 maxChebyPoly = 0;
+    unsigned int maxChebyPoly = 0;
 
     // Determine how many Chebyshev polynomials
@@ -542,5 +560,8 @@
     evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f]
  *****************************************************************************/
-float psGaussian(float x, float mean, float sigma, bool normal)
+float psGaussian(float x,
+                 float mean,
+                 float sigma,
+                 bool normal)
 {
     psF32 tmp = 1.0;
@@ -568,5 +589,7 @@
  *****************************************************************************/
 #define PS_XXX_GAUSSIAN_SEED 1995
-psVector* p_psGaussianDev(psF32 mean, psF32 sigma, psS32 Npts)
+psVector* p_psGaussianDev(psF32 mean,
+                          psF32 sigma,
+                          unsigned int Npts)
 {
     PS_ASSERT_INT_NONNEGATIVE(Npts, NULL);
@@ -575,5 +598,5 @@
     psRandom *r = psRandomAlloc(PS_RANDOM_TAUS, PS_XXX_GAUSSIAN_SEED);
     psVector* gauss = psVectorAlloc(Npts, PS_TYPE_F32);
-    for (psS32 i = 0; i < Npts; i++) {
+    for (unsigned int i = 0; i < Npts; i++) {
         gauss->data.F32[i] = mean + p_psRandomGaussian(r, sigma);
     }
@@ -586,10 +609,10 @@
     This routine must allocate memory for the polynomial structures.
  *****************************************************************************/
-psPolynomial1D* psPolynomial1DAlloc(int n,
+psPolynomial1D* psPolynomial1DAlloc(unsigned int n,
                                     psPolynomialType type)
 {
     PS_ASSERT_INT_POSITIVE(n, NULL);
 
-    int i = 0;
+    unsigned int i = 0;
     psPolynomial1D* newPoly = NULL;
 
@@ -601,5 +624,5 @@
     newPoly->coeff = psAlloc(n * sizeof(psF64));
     newPoly->coeffErr = psAlloc(n * sizeof(psF64));
-    newPoly->mask = (char *)psAlloc(n * sizeof(char));
+    newPoly->mask = (psMaskType *)psAlloc(n * sizeof(psMaskType));
     for (i = 0; i < n; i++) {
         newPoly->coeff[i] = 0.0;
@@ -611,5 +634,6 @@
 }
 
-psPolynomial2D* psPolynomial2DAlloc( int nX,  int nY,
+psPolynomial2D* psPolynomial2DAlloc( unsigned int nX,
+                                     unsigned int nY,
                                      psPolynomialType type)
 {
@@ -617,6 +641,6 @@
     PS_ASSERT_INT_POSITIVE(nY, NULL);
 
-    int x = 0;
-    int y = 0;
+    unsigned int x = 0;
+    unsigned int y = 0;
     psPolynomial2D* newPoly = NULL;
 
@@ -630,9 +654,9 @@
     newPoly->coeff = psAlloc(nX * sizeof(psF64 *));
     newPoly->coeffErr = psAlloc(nX * sizeof(psF64 *));
-    newPoly->mask = (char **)psAlloc(nX * sizeof(char *));
+    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] = (char *)psAlloc(nY * sizeof(char));
+        newPoly->mask[x] = (psMaskType *)psAlloc(nY * sizeof(psMaskType));
     }
     for (x = 0; x < nX; x++) {
@@ -647,5 +671,7 @@
 }
 
-psPolynomial3D* psPolynomial3DAlloc( int nX,  int nY,  int nZ,
+psPolynomial3D* psPolynomial3DAlloc( unsigned int nX,
+                                     unsigned int nY,
+                                     unsigned int nZ,
                                      psPolynomialType type)
 {
@@ -654,7 +680,7 @@
     PS_ASSERT_INT_POSITIVE(nZ, NULL);
 
-    psS32 x = 0;
-    psS32 y = 0;
-    psS32 z = 0;
+    unsigned int x = 0;
+    unsigned int y = 0;
+    unsigned int z = 0;
     psPolynomial3D* newPoly = NULL;
 
@@ -669,13 +695,13 @@
     newPoly->coeff = psAlloc(nX * sizeof(psF64 **));
     newPoly->coeffErr = psAlloc(nX * sizeof(psF64 **));
-    newPoly->mask = (char ***)psAlloc(nX * sizeof(char **));
+    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] = (char **)psAlloc(nY * sizeof(char *));
+        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] = (char *)psAlloc(nZ * sizeof(char));
+            newPoly->mask[x][y] = (psMaskType *)psAlloc(nZ * sizeof(psMaskType));
         }
     }
@@ -693,5 +719,8 @@
 }
 
-psPolynomial4D* psPolynomial4DAlloc( int nX,  int nY,  int nZ,  int nT,
+psPolynomial4D* psPolynomial4DAlloc( unsigned int nX,
+                                     unsigned int nY,
+                                     unsigned int nZ,
+                                     unsigned int nT,
                                      psPolynomialType type)
 {
@@ -701,8 +730,8 @@
     PS_ASSERT_INT_POSITIVE(nT, NULL);
 
-    psS32 x = 0;
-    psS32 y = 0;
-    psS32 z = 0;
-    psS32 t = 0;
+    unsigned int x = 0;
+    unsigned int y = 0;
+    unsigned int z = 0;
+    unsigned int t = 0;
     psPolynomial4D* newPoly = NULL;
 
@@ -718,17 +747,17 @@
     newPoly->coeff = psAlloc(nX * sizeof(psF64 ***));
     newPoly->coeffErr = psAlloc(nX * sizeof(psF64 ***));
-    newPoly->mask = (char ****)psAlloc(nX * sizeof(char ***));
+    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] = (char ***)psAlloc(nY * sizeof(char **));
+        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] = (char **)psAlloc(nZ * sizeof(char *));
+            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] = (char *)psAlloc(nT * sizeof(char));
+                newPoly->mask[x][y][z] = (psMaskType *)psAlloc(nT * sizeof(psMaskType));
             }
         }
@@ -749,5 +778,6 @@
 }
 
-psF64 psPolynomial1DEval(const psPolynomial1D* poly, psF64 x)
+psF64 psPolynomial1DEval(const psPolynomial1D* poly,
+                         psF64 x)
 {
     PS_ASSERT_POLY_NON_NULL(poly, NAN);
@@ -775,5 +805,5 @@
 
     tmp = psVectorAlloc(x->n, PS_TYPE_F64);
-    for (psS32 i=0;i<x->n;i++) {
+    for (unsigned int i=0;i<x->n;i++) {
         tmp->data.F64[i] = psPolynomial1DEval(poly, x->data.F64[i]);
     }
@@ -782,5 +812,7 @@
 }
 
-psF64 psPolynomial2DEval(const psPolynomial2D* poly, psF64 x, psF64 y)
+psF64 psPolynomial2DEval(const psPolynomial2D* poly,
+                         psF64 x,
+                         psF64 y)
 {
     PS_ASSERT_POLY_NON_NULL(poly, NAN);
@@ -810,5 +842,5 @@
 
     psVector *tmp;
-    psS32 vecLen=x->n;
+    unsigned int vecLen=x->n;
 
     // Determine the length of the output vector to by the minimum of the x,y vectors
@@ -821,5 +853,5 @@
 
     // Evaluate the polynomial at the specified points
-    for (psS32 i=0; i<vecLen; i++) {
+    for (unsigned int i=0; i<vecLen; i++) {
         tmp->data.F64[i] = psPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]);
     }
@@ -829,5 +861,8 @@
 }
 
-psF64 psPolynomial3DEval(const psPolynomial3D* poly, psF64 x, psF64 y, psF64 z)
+psF64 psPolynomial3DEval(const psPolynomial3D* poly,
+                         psF64 x,
+                         psF64 y,
+                         psF64 z)
 {
     PS_ASSERT_POLY_NON_NULL(poly, NAN);
@@ -860,5 +895,5 @@
 
     psVector *tmp;
-    psS32 vecLen=x->n;
+    unsigned int vecLen=x->n;
 
     // Determine the length of output vector from min of the input vectors
@@ -874,5 +909,5 @@
 
     // Evaluate polynomial
-    for (psS32 i = 0; i < vecLen; i++) {
+    for (unsigned int i = 0; i < vecLen; i++) {
         tmp->data.F64[i] = psPolynomial3DEval(poly,
                                               x->data.F64[i],
@@ -885,5 +920,9 @@
 }
 
-psF64 psPolynomial4DEval(const psPolynomial4D* poly, psF64 x, psF64 y, psF64 z, psF64 t)
+psF64 psPolynomial4DEval(const psPolynomial4D* poly,
+                         psF64 x,
+                         psF64 y,
+                         psF64 z,
+                         psF64 t)
 {
     PS_ASSERT_POLY_NON_NULL(poly, NAN);
@@ -918,5 +957,5 @@
 
     psVector *tmp;
-    psS32 vecLen=x->n;
+    unsigned int vecLen=x->n;
 
     // Determine output vector size from min of input vectors
@@ -935,5 +974,5 @@
 
     // Evaluate polynomial
-    for (psS32 i = 0; i < vecLen; i++) {
+    for (unsigned int i = 0; i < vecLen; i++) {
         tmp->data.F64[i] = psPolynomial4DEval(poly,
                                               x->data.F64[i],
Index: /trunk/psLib/src/math/psPolynomial.h
===================================================================
--- /trunk/psLib/src/math/psPolynomial.h	(revision 5065)
+++ /trunk/psLib/src/math/psPolynomial.h	(revision 5066)
@@ -11,6 +11,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-08 00:02:48 $
+ *  @version $Revision: 1.56 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-19 19:53:13 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -56,5 +56,5 @@
     psF32 mean,                        ///< The mean of the Gaussian
     psF32 sigma,                       ///< The sigma of the Gaussian
-    psS32 Npts                         ///< The size of the vector
+    unsigned int Npts                  ///< The size of the vector
 );
 
@@ -73,8 +73,8 @@
 {
     psPolynomialType type;             ///< Polynomial type
-    int n;                             ///< Number of terms
+    unsigned int n;                    ///< Number of terms
     psF64 *coeff;                      ///< Coefficients
     psF64 *coeffErr;                   ///< Error in coefficients
-    char *mask;                        ///< Coefficient mask
+    psMaskType *mask;                  ///< Coefficient mask
 }
 psPolynomial1D;
@@ -84,9 +84,9 @@
 {
     psPolynomialType type;             ///< Polynomial type
-    int nX;                            ///< Number of terms in x
-    int nY;                            ///< Number of terms in y
+    unsigned int nX;                   ///< Number of terms in x
+    unsigned int nY;                   ///< Number of terms in y
     psF64 **coeff;                     ///< Coefficients
     psF64 **coeffErr;                  ///< Error in coefficients
-    char **mask;                       ///< Coefficients mask
+    psMaskType **mask;                 ///< Coefficients mask
 }
 psPolynomial2D;
@@ -96,10 +96,10 @@
 {
     psPolynomialType type;             ///< Polynomial type
-    int nX;                           ///< Number of terms in x
-    int nY;                            ///< Number of terms in y
-    int nZ;                           ///< Number of terms in z
+    unsigned int nX;                   ///< Number of terms in x
+    unsigned int nY;                   ///< Number of terms in y
+    unsigned int nZ;                   ///< Number of terms in z
     psF64 ***coeff;                    ///< Coefficients
     psF64 ***coeffErr;                 ///< Error in coefficients
-    char ***mask;                      ///< Coefficients mask
+    psMaskType ***mask;                ///< Coefficients mask
 }
 psPolynomial3D;
@@ -109,11 +109,11 @@
 {
     psPolynomialType type;             ///< Polynomial type
-    int nX;                            ///< Number of terms in x
-    int nY;                            ///< Number of terms in y
-    int nZ;                            ///< Number of terms in z
-    int nT;                            ///< Number of terms in t
+    unsigned int nX;                   ///< Number of terms in x
+    unsigned int nY;                   ///< Number of terms in y
+    unsigned int nZ;                   ///< Number of terms in z
+    unsigned int nT;                   ///< Number of terms in t
     psF64 ****coeff;                   ///< Coefficients
     psF64 ****coeffErr;                ///< Error in coefficients
-    char ****mask;                     ///< Coefficients mask
+    psMaskType ****mask;               ///< Coefficients mask
 }
 psPolynomial4D;
@@ -125,5 +125,5 @@
  */
 psPolynomial1D* psPolynomial1DAlloc(
-    int n,                             ///< Number of terms
+    unsigned int n,                    ///< Number of terms
     psPolynomialType type              ///< Polynomial Type
 );
@@ -134,6 +134,6 @@
  */
 psPolynomial2D* psPolynomial2DAlloc(
-    int nX,                   ///< Number of terms in x
-    int nY,                   ///< Number of terms in y
+    unsigned int nX,                   ///< Number of terms in x
+    unsigned int nY,                   ///< Number of terms in y
     psPolynomialType type              ///< Polynomial Type
 );
@@ -144,7 +144,7 @@
  */
 psPolynomial3D* psPolynomial3DAlloc(
-    int nX,                            ///< Number of terms in x
-    int nY,                            ///< Number of terms in y
-    int nZ,                            ///< Number of terms in z
+    unsigned int nX,                   ///< Number of terms in x
+    unsigned int nY,                   ///< Number of terms in y
+    unsigned int nZ,                   ///< Number of terms in z
     psPolynomialType type              ///< Polynomial Type
 );
@@ -155,8 +155,8 @@
  */
 psPolynomial4D* psPolynomial4DAlloc(
-    int nX,                            ///< Number of terms in x
-    int nY,                            ///< Number of terms in y
-    int nZ,                            ///< Number of terms in z
-    int nT,                            ///< Number of terms in t
+    unsigned int nX,                   ///< Number of terms in x
+    unsigned int nY,                   ///< Number of terms in y
+    unsigned int nZ,                   ///< Number of terms in z
+    unsigned int nT,                   ///< Number of terms in t
     psPolynomialType type              ///< Polynomial Type
 );
Index: /trunk/psLib/src/math/psSpline.c
===================================================================
--- /trunk/psLib/src/math/psSpline.c	(revision 5065)
+++ /trunk/psLib/src/math/psSpline.c	(revision 5066)
@@ -7,6 +7,6 @@
 *  splines.
 *
-*  @version $Revision: 1.123 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-09-11 22:18:40 $
+*  @version $Revision: 1.124 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-09-19 19:53:13 $
 *
 *
@@ -43,6 +43,6 @@
 /*****************************************************************************/
 static void spline1DFree(psSpline1D *tmpSpline);
-static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x);
-static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x);
+static unsigned int vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x);
+static unsigned int vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x);
 
 /*****************************************************************************/
@@ -69,5 +69,5 @@
 static void spline1DFree(psSpline1D *tmpSpline)
 {
-    psS32 i;
+    unsigned int i;
 
     if (tmpSpline == NULL) {
@@ -101,10 +101,10 @@
 static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \
                                      ps##TYPE *range, \
-                                     psS32 n, \
+                                     unsigned int n, \
                                      ps##TYPE x) \
 { \
     \
-    psS32 i; \
-    psS32 m; \
+    unsigned int i; \
+    unsigned int m; \
     static psVector *p = NULL; \
     p = psVectorRecycle(p, n, PS_TYPE_##TYPE); \
@@ -169,6 +169,6 @@
 static psF32 interpolate1DF32(psF32 *domain,
                               psF32 *range,
-                              psS32 n,
-                              psS32 order,
+                              unsigned int n,
+                              unsigned int order,
                               psF32 x)
 {
@@ -178,5 +178,5 @@
 
     psS32 binNum;
-    psS32 numIntPoints = order+1;
+    unsigned int numIntPoints = order+1;
     psS32 origin;
 
@@ -229,9 +229,9 @@
             "---- calculateSecondDerivs() begin ----\n");
 
-    psS32 i;
-    psS32 k;
+    unsigned int i;
+    unsigned int k;
     psF32 sig;
     psF32 p;
-    psS32 n = y->n;
+    unsigned int n = y->n;
     psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32));
     psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32));
@@ -253,9 +253,9 @@
 
         psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
-                "X[%d] is %f\n", i, X[i]);
+                "X[%u] is %f\n", i, X[i]);
         psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
-                "Y[%d] is %f\n", i, Y[i]);
+                "Y[%u] is %f\n", i, Y[i]);
         psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
-                "u[%d] is %f\n", i, u[i]);
+                "u[%u] is %f\n", i, u[i]);
     }
 
@@ -268,5 +268,5 @@
 
         psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
-                "derivs2[%d] is %f\n", k, derivs2[k]);
+                "derivs2[%u] is %f\n", k, derivs2[k]);
     }
 
@@ -314,5 +314,5 @@
 XXX: Assumes mySpline->knots is psF32.  Must add psU32 and psF64.
  *****************************************************************************/
-psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline,     ///< The spline which will be generated.
+psSpline1D *psVectorFitSpline1D(psSpline1D *spline,     ///< The spline which will be generated.
                                 const psVector* x,        ///< Ordinates (or NULL to just use the indices)
                                 const psVector* y,        ///< Coordinates
@@ -321,13 +321,13 @@
     PS_ASSERT_VECTOR_NON_NULL(y, NULL);
     PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
-    if (mySpline != NULL) {
-        PS_ASSERT_VECTOR_TYPE(mySpline->knots, PS_TYPE_F32, NULL);
+    if (spline != NULL) {
+        PS_ASSERT_VECTOR_TYPE(spline->knots, PS_TYPE_F32, NULL);
     }
     psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
             "---- psVectorFitSpline1D() begin ----\n");
-    psS32 numSplines = (y->n)-1;
+    unsigned int numSplines = (y->n)-1;
     psF32 tmp;
     psF32 H;
-    psS32 i;
+    unsigned int i;
     psF32 slope;
     psVector *x32 = NULL;
@@ -364,17 +364,17 @@
         This can not be implemented until SDR states what order spline should be
         created.
-        Should we error if mySpline is not NULL?
+        Should we error if spline is not NULL?
         Should we error if mySPline is not NULL?
     */
-    if (mySpline == NULL) {
-        mySpline = psSpline1DAllocGeneric(x32, 3);
-    }
-    PS_ASSERT_PTR_NON_NULL(mySpline, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(mySpline->n, NULL);
-
-    if (y32->n != (1 + mySpline->n)) {
+    if (spline == NULL) {
+        spline = psSpline1DAllocGeneric(x32, 3);
+    }
+    PS_ASSERT_PTR_NON_NULL(spline, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(spline->n, NULL);
+
+    if (y32->n != (1 + spline->n)) {
         psError(PS_ERR_BAD_PARAMETER_SIZE, true,
-                "data size / spline size mismatch (%d %d)\n",
-                y32->n, mySpline->n);
+                "data size / spline size mismatch (%u %u)\n",
+                y32->n, spline->n);
         return(NULL);
     }
@@ -382,5 +382,205 @@
     // If these are linear splines, which means their polynomials will have
     // two coefficients, then we do the simple calculation.
-    if (2 == (mySpline->spline[0])->n) {
+    if (2 == (spline->spline[0])->n) {
+        for (i=0;i<spline->n;i++) {
+            slope = (y32->data.F32[i+1] - y32->data.F32[i]) /
+                    (spline->knots->data.F32[i+1] - spline->knots->data.F32[i]);
+            (spline->spline[i])->coeff[0] = y32->data.F32[i] -
+                                            (slope * spline->knots->data.F32[i]);
+
+            (spline->spline[i])->coeff[1] = slope;
+            psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
+                    "---- spline %u coeffs are (%f, %f)\n", i,
+                    (spline->spline[i])->coeff[0],
+                    (spline->spline[i])->coeff[1]);
+        }
+        psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
+                "---- Exiting psVectorFitSpline1D()()\n");
+        return((psSpline1D *) spline);
+    }
+
+    // Check if these are cubic splines (n==4).  If not, psError.
+    if (4 != (spline->spline[0])->n) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
+                "Don't know how to generate %u-order splines.",
+                (spline->spline[0])->n-1);
+        return(NULL);
+    }
+
+    // If we get here, then we know these are cubic splines.  We first
+    // generate the second derivatives at each data point.
+    spline->p_psDeriv2 = calculateSecondDerivs(x32, y32);
+    for (i=0;i<y32->n;i++)
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "Second deriv[%u] is %f\n", i, spline->p_psDeriv2[i]);
+
+    // We generate the coefficients of the spline polynomials.  I can't
+    // concisely explain how this code works.  See above function comments
+    // and Numerical Recipes in C.
+    for (i=0;i<numSplines;i++) {
+        H = x32->data.F32[i+1] - x32->data.F32[i];
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
+                "x data (%f - %f) (%f)\n",
+                x32->data.F32[i],
+                x32->data.F32[i+1], H);
+        //
+        // ******** Calculate 0-order term ********
+        //
+        // From (1)
+        (spline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H);
+        // From (2)
+        ((spline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H);
+        // From (3)
+        tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
+        tmp-= (x32->data.F32[i+1] / H);
+        tmp*= (spline->p_psDeriv2)[i] * H * H / 6.0;
+        ((spline->spline[i])->coeff[0])+= tmp;
+        // From (4)
+        tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
+        tmp+= (x32->data.F32[i] / H);
+        tmp*= (spline->p_psDeriv2)[i+1] * H * H / 6.0;
+        ((spline->spline[i])->coeff[0])+= tmp;
+
+        //
+        // ******** Calculate 1-order term ********
+        //
+        // From (1)
+        (spline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H;
+        // From (2)
+        ((spline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H);
+        // From (3)
+        tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
+        tmp+= (1.0 / H);
+        tmp*= ((spline->p_psDeriv2)[i]) * H * H / 6.0;
+        ((spline->spline[i])->coeff[1])+= tmp;
+        // From (4)
+        tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
+        tmp-= (1.0 / H);
+        tmp*= ((spline->p_psDeriv2)[i+1]) * H * H / 6.0;
+        ((spline->spline[i])->coeff[1])+= tmp;
+
+        //
+        // ******** Calculate 2-order term ********
+        //
+        // From (3)
+        (spline->spline[i])->coeff[2] = ((spline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H);
+        // From (4)
+        ((spline->spline[i])->coeff[2])-= (((spline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H));
+
+        //
+        // ******** Calculate 3-order term ********
+        //
+        // From (3)
+        (spline->spline[i])->coeff[3] = -((spline->p_psDeriv2)[i]) / (6.0 * H);
+        // From (4)
+        ((spline->spline[i])->coeff[3])+=  ((spline->p_psDeriv2)[i+1]) / (6.0 * H);
+
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "(spline->spline[%u])->coeff[0] is %f\n", i, (spline->spline[i])->coeff[0]);
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "(spline->spline[%u])->coeff[1] is %f\n", i, (spline->spline[i])->coeff[1]);
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "(spline->spline[%u])->coeff[2] is %f\n", i, (spline->spline[i])->coeff[2]);
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "(spline->spline[%u])->coeff[3] is %f\n", i, (spline->spline[i])->coeff[3]);
+
+    }
+
+    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
+            "---- psVectorFitSpline1D() end ----\n");
+    return(spline);
+}
+
+
+
+
+
+
+
+
+
+/*****************************************************************************
+psVectorFitSpline1DNEW(): given a psSpline1D data structure and a set of x/y
+ 
+xF32 and yF32 are internal psVectors which are used to hold the psF32 versions
+of the input data, if necessary.  xPtr and yPtr are pointers to either xF32 or
+the x argument.  All computation is done on xPtr and yPtr.  xF32 and yF32 will
+simply be psFree() at the end.
+ 
+XXX: nKnots makes no sense.  This number is always equal to the size of the x
+an y vectors.
+ 
+XXX: How do we specify the spline order?  For now, order=3.
+ *****************************************************************************/
+#define PS_XXX_SPLINE_ORDER 3
+psSpline1D *psVectorFitSpline1DNEW(const psVector* x,        ///< Ordinates (or NULL to just use the indices)
+                                   const psVector* y,        ///< Coordinates
+                                   unsigned int nKnots)
+{
+    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, "---- psVectorFitSpline1DNEW() begin ----\n");
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
+    //    PS_ASSERT_INT_EQUAL(y->n, nKnots);
+
+    //
+    // The following code ensures that xPtr points to a psF32 version of the
+    // ordinate data.
+    //
+    psVector *xF32 = NULL;
+    psVector *xPtr = NULL;
+    if (x != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
+        if (PS_TYPE_F64 == x->type.type) {
+            xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
+            for (unsigned int i = 0 ; i < x->n ; i++) {
+                xF32->data.F32[i] = (psF32) x->data.F64[i];
+            }
+            xPtr = xF32;
+        } else if (PS_TYPE_F32 == x->type.type) {
+            xPtr = (psVector *) x;
+        } else {
+            psError(PS_ERR_UNKNOWN, true, "psVector x is wrong type.\n");
+            return(NULL);
+        }
+    } else {
+        // If x==NULL, create an x32 vector with x values set to (0:n).
+        xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
+        for (unsigned int i = 0 ; i < x->n ; i++) {
+            xF32->data.F32[i] = (psF32) i;
+        }
+        xPtr = xF32;
+    }
+
+    //
+    // If y is of type psF64, then create a new vector yF32 and convert the
+    // y elements.  Regardless of y's type, we create a yPtr which will be
+    // used in the remainder of this function.
+    //
+    psVector *yF32 = NULL;
+    psVector *yPtr = NULL;
+    if (PS_TYPE_F64 == y->type.type) {
+        yF32 = psVectorAlloc(y->n, PS_TYPE_F32);
+        for (unsigned int i = 0 ; i < y->n ; i++) {
+            yF32->data.F32[i] = (psF32) y->data.F64[i];
+        }
+        yPtr = yF32;
+    } else {
+        yPtr = (psVector *) y;
+    }
+
+    psSpline1D *mySpline = psSpline1DAllocGeneric(xPtr, PS_XXX_SPLINE_ORDER);
+
+    unsigned int numSplines = nKnots - 1;
+    psF32 tmp;
+    psF32 H;
+    unsigned int i;
+    psF32 slope;
+    // XXX: get rid of x32 and y32 (this is from old code)
+    psVector *x32 = xPtr;
+    psVector *y32 = yPtr;
+
+    // If these are linear splines, which means their polynomials will have
+    // two coefficients, then we do the simple calculation.
+    if (1 == PS_XXX_SPLINE_ORDER) {
         for (i=0;i<mySpline->n;i++) {
             slope = (y32->data.F32[i+1] - y32->data.F32[i]) /
@@ -391,5 +591,5 @@
             (mySpline->spline[i])->coeff[1] = slope;
             psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
-                    "---- mySpline %d coeffs are (%f, %f)\n", i,
+                    "---- mySpline %u coeffs are (%f, %f)\n", i,
                     (mySpline->spline[i])->coeff[0],
                     (mySpline->spline[i])->coeff[1]);
@@ -400,22 +600,27 @@
     }
 
+    //
     // Check if these are cubic splines (n==4).  If not, psError.
-    if (4 != (mySpline->spline[0])->n) {
+    //
+    if (3 != PS_XXX_SPLINE_ORDER) {
         psError(PS_ERR_BAD_PARAMETER_SIZE, true,
-                "Don't know how to generate %d-order splines.",
-                (mySpline->spline[0])->n-1);
+                "Don't know how to generate %d-order splines.", PS_XXX_SPLINE_ORDER);
         return(NULL);
     }
 
+    //
     // If we get here, then we know these are cubic splines.  We first
     // generate the second derivatives at each data point.
+    //
     mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32);
     for (i=0;i<y32->n;i++)
         psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
-
+                "Second deriv[%u] is %f\n", i, mySpline->p_psDeriv2[i]);
+
+    //
     // We generate the coefficients of the spline polynomials.  I can't
     // concisely explain how this code works.  See above function comments
     // and Numerical Recipes in C.
+    //
     for (i=0;i<numSplines;i++) {
         H = x32->data.F32[i+1] - x32->data.F32[i];
@@ -477,216 +682,11 @@
 
         psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
+                "(mySpline->spline[%u])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
         psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
+                "(mySpline->spline[%u])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
         psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
+                "(mySpline->spline[%u])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
         psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
-
-    }
-
-    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
-            "---- psVectorFitSpline1D() end ----\n");
-    return(mySpline);
-}
-
-
-
-
-
-
-
-
-
-/*****************************************************************************
-psVectorFitSpline1DNEW(): given a psSpline1D data structure and a set of x/y
- 
-xF32 and yF32 are internal psVectors which are used to hold the psF32 versions
-of the input data, if necessary.  xPtr and yPtr are pointers to either xF32 or
-the x argument.  All computation is done on xPtr and yPtr.  xF32 and yF32 will
-simply be psFree() at the end.
- 
-XXX: nKnots makes no sense.  This number is always equal to the size of the x
-an y vectors.
- 
-XXX: How do we specify the spline order?  For now, order=3.
- *****************************************************************************/
-#define PS_XXX_SPLINE_ORDER 3
-psSpline1D *psVectorFitSpline1DNEW(const psVector* x,        ///< Ordinates (or NULL to just use the indices)
-                                   const psVector* y,        ///< Coordinates
-                                   int nKnots)
-{
-    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, "---- psVectorFitSpline1DNEW() begin ----\n");
-    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
-    //    PS_ASSERT_INT_EQUAL(y->n, nKnots);
-
-    //
-    // The following code ensures that xPtr points to a psF32 version of the
-    // ordinate data.
-    //
-    psVector *xF32 = NULL;
-    psVector *xPtr = NULL;
-    if (x != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
-        if (PS_TYPE_F64 == x->type.type) {
-            xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
-            for (psS32 i = 0 ; i < x->n ; i++) {
-                xF32->data.F32[i] = (psF32) x->data.F64[i];
-            }
-            xPtr = xF32;
-        } else if (PS_TYPE_F32 == x->type.type) {
-            xPtr = (psVector *) x;
-        } else {
-            psError(PS_ERR_UNKNOWN, true, "psVector x is wrong type.\n");
-            return(NULL);
-        }
-    } else {
-        // If x==NULL, create an x32 vector with x values set to (0:n).
-        xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
-        for (psS32 i = 0 ; i < x->n ; i++) {
-            xF32->data.F32[i] = (psF32) i;
-        }
-        xPtr = xF32;
-    }
-
-    //
-    // If y is of type psF64, then create a new vector yF32 and convert the
-    // y elements.  Regardless of y's type, we create a yPtr which will be
-    // used in the remainder of this function.
-    //
-    psVector *yF32 = NULL;
-    psVector *yPtr = NULL;
-    if (PS_TYPE_F64 == y->type.type) {
-        yF32 = psVectorAlloc(y->n, PS_TYPE_F32);
-        for (psS32 i = 0 ; i < y->n ; i++) {
-            yF32->data.F32[i] = (psF32) y->data.F64[i];
-        }
-        yPtr = yF32;
-    } else {
-        yPtr = (psVector *) y;
-    }
-
-    psSpline1D *mySpline = psSpline1DAllocGeneric(xPtr, PS_XXX_SPLINE_ORDER);
-
-    psS32 numSplines = nKnots - 1;
-    psF32 tmp;
-    psF32 H;
-    psS32 i;
-    psF32 slope;
-    // XXX: get rid of x32 and y32 (this is from old code)
-    psVector *x32 = xPtr;
-    psVector *y32 = yPtr;
-
-    // If these are linear splines, which means their polynomials will have
-    // two coefficients, then we do the simple calculation.
-    if (1 == PS_XXX_SPLINE_ORDER) {
-        for (i=0;i<mySpline->n;i++) {
-            slope = (y32->data.F32[i+1] - y32->data.F32[i]) /
-                    (mySpline->knots->data.F32[i+1] - mySpline->knots->data.F32[i]);
-            (mySpline->spline[i])->coeff[0] = y32->data.F32[i] -
-                                              (slope * mySpline->knots->data.F32[i]);
-
-            (mySpline->spline[i])->coeff[1] = slope;
-            psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
-                    "---- mySpline %d coeffs are (%f, %f)\n", i,
-                    (mySpline->spline[i])->coeff[0],
-                    (mySpline->spline[i])->coeff[1]);
-        }
-        psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
-                "---- Exiting psVectorFitSpline1D()()\n");
-        return((psSpline1D *) mySpline);
-    }
-
-    //
-    // Check if these are cubic splines (n==4).  If not, psError.
-    //
-    if (3 != PS_XXX_SPLINE_ORDER) {
-        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
-                "Don't know how to generate %d-order splines.", PS_XXX_SPLINE_ORDER);
-        return(NULL);
-    }
-
-    //
-    // If we get here, then we know these are cubic splines.  We first
-    // generate the second derivatives at each data point.
-    //
-    mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32);
-    for (i=0;i<y32->n;i++)
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
-
-    //
-    // We generate the coefficients of the spline polynomials.  I can't
-    // concisely explain how this code works.  See above function comments
-    // and Numerical Recipes in C.
-    //
-    for (i=0;i<numSplines;i++) {
-        H = x32->data.F32[i+1] - x32->data.F32[i];
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
-                "x data (%f - %f) (%f)\n",
-                x32->data.F32[i],
-                x32->data.F32[i+1], H);
-        //
-        // ******** Calculate 0-order term ********
-        //
-        // From (1)
-        (mySpline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H);
-        // From (2)
-        ((mySpline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H);
-        // From (3)
-        tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
-        tmp-= (x32->data.F32[i+1] / H);
-        tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0;
-        ((mySpline->spline[i])->coeff[0])+= tmp;
-        // From (4)
-        tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
-        tmp+= (x32->data.F32[i] / H);
-        tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0;
-        ((mySpline->spline[i])->coeff[0])+= tmp;
-
-        //
-        // ******** Calculate 1-order term ********
-        //
-        // From (1)
-        (mySpline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H;
-        // From (2)
-        ((mySpline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H);
-        // From (3)
-        tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
-        tmp+= (1.0 / H);
-        tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0;
-        ((mySpline->spline[i])->coeff[1])+= tmp;
-        // From (4)
-        tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
-        tmp-= (1.0 / H);
-        tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
-        ((mySpline->spline[i])->coeff[1])+= tmp;
-
-        //
-        // ******** Calculate 2-order term ********
-        //
-        // From (3)
-        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H);
-        // From (4)
-        ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H));
-
-        //
-        // ******** Calculate 3-order term ********
-        //
-        // From (3)
-        (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
-        // From (4)
-        ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
-
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
-        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
-                "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
+                "(mySpline->spline[%u])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
 
     }
@@ -726,6 +726,6 @@
 XXX: What should be the default type for knots be?  psF32 is assumed.
  *****************************************************************************/
-psSpline1D *psSpline1DAlloc(int numSplines,
-                            int order,
+psSpline1D *psSpline1DAlloc(unsigned int numSplines,
+                            unsigned int order,
                             float min,
                             float max)
@@ -743,5 +743,5 @@
     //
     tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
-    for (psS32 i=0;i<numSplines;i++) {
+    for (unsigned int i=0;i<numSplines;i++) {
         (tmpSpline->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
     }
@@ -756,5 +756,5 @@
     psF32 width = (max - min) / ((psF32) numSplines);
     tmpSpline->knots->data.F32[0] = min;
-    for (psS32 i=1;i<numSplines;i++) {
+    for (unsigned int i=1;i<numSplines;i++) {
         tmpSpline->knots->data.F32[i] = min + (width * (psF32) i);
     }
@@ -773,9 +773,9 @@
 
     if (in->type.type == PS_TYPE_F32) {
-        for (psS32 i = 0 ; i < in->n ; i++) {
+        for (unsigned int i = 0 ; i < in->n ; i++) {
             out->data.F32[i] = in->data.F32[i];
         }
     } else if (in->type.type == PS_TYPE_F64) {
-        for (psS32 i = 0 ; i < in->n ; i++) {
+        for (unsigned int i = 0 ; i < in->n ; i++) {
             out->data.F64[i] = in->data.F64[i];
         }
@@ -791,5 +791,5 @@
  *****************************************************************************/
 psSpline1D *psSpline1DAllocGeneric(const psVector *bounds,
-                                   int order)
+                                   unsigned int order)
 {
     PS_ASSERT_VECTOR_NON_NULL(bounds, NULL);
@@ -799,5 +799,5 @@
 
     psSpline1D *tmpSpline = (psSpline1D *) psAlloc(sizeof(psSpline1D));
-    psS32 numSplines = bounds->n - 1;
+    unsigned int numSplines = bounds->n - 1;
     tmpSpline->n = numSplines;
 
@@ -807,5 +807,5 @@
     //
     tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
-    for (psS32 i=0;i<numSplines;i++) {
+    for (unsigned int i=0;i<numSplines;i++) {
         (tmpSpline->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
     }
@@ -818,5 +818,5 @@
     // XXX:Ensure that the knots are monotonic.
     //
-    for (psS32 i=0;i<bounds->n-1;i++) {
+    for (unsigned int i=0;i<bounds->n-1;i++) {
         if (FLT_EPSILON >= fabs(bounds->data.F32[i+1]-bounds->data.F32[i])) {
             psError(PS_ERR_UNKNOWN, true, "data points must be distinct ([%d] %f %f)\n", i, bounds->data.F32[i], bounds->data.F32[i+1]);
@@ -839,7 +839,7 @@
  *****************************************************************************/
 #define FUNC_MACRO_VECTOR_BIN_DISECT(TYPE) \
-static psS32 vectorBinDisect##TYPE(ps##TYPE *bins, \
-                                   psS32 numBins, \
-                                   ps##TYPE x) \
+static unsigned int vectorBinDisect##TYPE(ps##TYPE *bins, \
+        psS32 numBins, \
+        ps##TYPE x) \
 { \
     psS32 min; \
@@ -906,6 +906,6 @@
 XXX: Assert that the psVector and psScalar have the same type.
  *****************************************************************************/
-psS32 p_psVectorBinDisect(psVector *bins,
-                          psScalar *x)
+unsigned int p_psVectorBinDisect(psVector *bins,
+                                 psScalar *x)
 {
     PS_ASSERT_VECTOR_NON_NULL(bins, -4);
@@ -972,5 +972,5 @@
 psScalar *p_psVectorInterpolate(psVector *domain,
                                 psVector *range,
-                                int order,
+                                unsigned int order,
                                 psScalar *x)
 {
Index: /trunk/psLib/src/math/psSpline.h
===================================================================
--- /trunk/psLib/src/math/psSpline.h	(revision 5065)
+++ /trunk/psLib/src/math/psSpline.h	(revision 5066)
@@ -10,6 +10,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-11 22:18:40 $
+ *  @version $Revision: 1.56 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-19 19:53:13 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -35,5 +35,5 @@
 typedef struct
 {
-    int n;                             ///< The number of spline pieces
+    unsigned int n;                    ///< The number of spline pieces
     psPolynomial1D **spline;           ///< An array of n pointers to the spline polynomials
     psVector *knots;                   ///< The boundaries between each spline piece.  Size is n+1.
@@ -51,6 +51,6 @@
  */
 psSpline1D *psSpline1DAlloc(
-    int n,                             ///< Number of spline polynomials
-    int order,                         ///< Order of spline polynomials
+    unsigned int n,                    ///< Number of spline polynomials
+    unsigned int order,                ///< Order of spline polynomials
     float min,                         ///< Lower boundary value of spline polynomials
     float max                          ///< Upper boundary value of spline polynomials
@@ -65,5 +65,5 @@
 psSpline1D *psSpline1DAllocGeneric(
     const psVector *bounds,            ///< Bounds for spline polynomials
-    int order                          ///< Order of spline polynomials
+    unsigned int order                 ///< Order of spline polynomials
 );
 
@@ -91,5 +91,5 @@
  *  @return psS32    corresponding index number of specified value
  */
-psS32 p_psVectorBinDisect(
+unsigned int p_psVectorBinDisect(
     psVector *bins,                    ///< Array of non-decreasing values
     psScalar *x                        ///< Target value to find
@@ -104,5 +104,5 @@
     psVector *domain,                  ///< Domain (x coords) for interpolation
     psVector *range,                   ///< Range (y coords) for interpolation
-    int order,                         ///< Order of interpolation function
+    unsigned int order,                ///< Order of interpolation function
     psScalar *x                        ///< Location at which to evaluate
 );
@@ -126,5 +126,5 @@
  */
 psSpline1D *psVectorFitSpline1D(
-    psSpline1D *mySpline,              ///< The spline which will be generated.
+    psSpline1D *spline,                ///< The spline which will be generated.
     const psVector* x,                 ///< Ordinates (or NULL to just use the indices)
     const psVector* y,                 ///< Coordinates
@@ -135,5 +135,5 @@
     const psVector* x,                 ///< Ordinates (or NULL to just use the indices)
     const psVector* y,                 ///< Coordinates
-    int nKnots
+    unsigned int nKnots                ///< Number of Knots
 );