Index: trunk/psLib/src/math/psPolynomial.c
===================================================================
--- trunk/psLib/src/math/psPolynomial.c	(revision 2778)
+++ trunk/psLib/src/math/psPolynomial.c	(revision 2788)
@@ -7,6 +7,6 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.74 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-12-21 20:42:07 $
+ *  @version $Revision: 1.75 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-12-22 05:09:32 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -55,5 +55,5 @@
 static void dPolynomial4DFree(psDPolynomial4D* myPoly);
 static void spline1DFree(psSpline1D *tmpSpline);
-static psS32 vectorBinDisectF32(float *bins,psS32 numBins,float x);
+static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x);
 static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x);
 
@@ -277,9 +277,9 @@
  
  *****************************************************************************/
-static float ordPolynomial1DEval(float x, const psPolynomial1D* myPoly)
+static psF32 ordPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly)
 {
     psS32 loop_x = 0;
-    float polySum = 0.0;
-    float xSum = 1.0;
+    psF32 polySum = 0.0;
+    psF32 xSum = 1.0;
 
     psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4,
@@ -307,5 +307,5 @@
 // XXX: How does the mask vector effect Crenshaw's formula?
 // XXX: We assume that x is scaled between -1.0 and 1.0;
-static float chebPolynomial1DEval(float x, const psPolynomial1D* myPoly)
+static psF32 chebPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly)
 {
     //    PS_FLOAT_CHECK_RANGE(x, -1.0, 1.0, 0.0);
@@ -314,5 +314,5 @@
     psS32 n;
     psS32 i;
-    float tmp;
+    psF32 tmp;
 
     n = myPoly->n;
@@ -337,5 +337,5 @@
     psS32 n;
     psS32 i;
-    float tmp;
+    psF32 tmp;
     psPolynomial1D **chebPolys = NULL;
 
@@ -355,5 +355,5 @@
 }
 
-static float ordPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
+static psF32 ordPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly)
 {
     PS_POLY_CHECK_NULL(myPoly, NAN);
@@ -361,7 +361,7 @@
     psS32 loop_x = 0;
     psS32 loop_y = 0;
-    float polySum = 0.0;
-    float xSum = 1.0;
-    float ySum = 1.0;
+    psF32 polySum = 0.0;
+    psF32 xSum = 1.0;
+    psF32 ySum = 1.0;
 
     for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
@@ -379,5 +379,5 @@
 }
 
-static float chebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
+static psF32 chebPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly)
 {
     PS_POLY_CHECK_NULL(myPoly, NAN);
@@ -386,5 +386,5 @@
     psS32 loop_y = 0;
     psS32 i = 0;
-    float polySum = 0.0;
+    psF32 polySum = 0.0;
     psPolynomial1D* *chebPolys = NULL;
     psS32 maxChebyPoly = 0;
@@ -414,13 +414,13 @@
 }
 
-static float ordPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
+static psF32 ordPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly)
 {
     psS32 loop_x = 0;
     psS32 loop_y = 0;
     psS32 loop_z = 0;
-    float polySum = 0.0;
-    float xSum = 1.0;
-    float ySum = 1.0;
-    float zSum = 1.0;
+    psF32 polySum = 0.0;
+    psF32 xSum = 1.0;
+    psF32 ySum = 1.0;
+    psF32 zSum = 1.0;
 
     for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
@@ -442,5 +442,5 @@
 }
 
-static float chebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
+static psF32 chebPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly)
 {
     psS32 loop_x = 0;
@@ -448,5 +448,5 @@
     psS32 loop_z = 0;
     psS32 i = 0;
-    float polySum = 0.0;
+    psF32 polySum = 0.0;
     psPolynomial1D* *chebPolys = NULL;
     psS32 maxChebyPoly = 0;
@@ -483,5 +483,5 @@
 }
 
-static float ordPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
+static psF32 ordPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly)
 {
     psS32 loop_w = 0;
@@ -489,9 +489,9 @@
     psS32 loop_y = 0;
     psS32 loop_z = 0;
-    float polySum = 0.0;
-    float wSum = 1.0;
-    float xSum = 1.0;
-    float ySum = 1.0;
-    float zSum = 1.0;
+    psF32 polySum = 0.0;
+    psF32 wSum = 1.0;
+    psF32 xSum = 1.0;
+    psF32 ySum = 1.0;
+    psF32 zSum = 1.0;
 
     for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
@@ -517,5 +517,5 @@
 }
 
-static float chebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
+static psF32 chebPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly)
 {
     psS32 loop_w = 0;
@@ -524,5 +524,5 @@
     psS32 loop_z = 0;
     psS32 i = 0;
-    float polySum = 0.0;
+    psF32 polySum = 0.0;
     psPolynomial1D* *chebPolys = NULL;
     psS32 maxChebyPoly = 0;
@@ -568,9 +568,9 @@
     Polynomial coefficients will be accessed in [w][x][y][z] fashion.
  *****************************************************************************/
-static double dOrdPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
+static psF64 dOrdPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly)
 {
     psS32 loop_x = 0;
-    double polySum = 0.0;
-    double xSum = 1.0;
+    psF64 polySum = 0.0;
+    psF64 xSum = 1.0;
 
     for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
@@ -586,10 +586,10 @@
 // XXX: You can do this without having to psAlloc() vector d.
 // XXX: How does the mask vector effect Crenshaw's formula?
-static double dChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
+static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly)
 {
     psVector *d;
     psS32 n;
     psS32 i;
-    double tmp;
+    psF64 tmp;
 
     n = myPoly->n;
@@ -611,11 +611,11 @@
 }
 
-static double dOrdPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
+static psF64 dOrdPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly)
 {
     psS32 loop_x = 0;
     psS32 loop_y = 0;
-    double polySum = 0.0;
-    double xSum = 1.0;
-    double ySum = 1.0;
+    psF64 polySum = 0.0;
+    psF64 xSum = 1.0;
+    psF64 ySum = 1.0;
 
     for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
@@ -633,10 +633,10 @@
 }
 
-static double dChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
+static psF64 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly)
 {
     psS32 loop_x = 0;
     psS32 loop_y = 0;
     psS32 i = 0;
-    double polySum = 0.0;
+    psF64 polySum = 0.0;
     psPolynomial1D* *chebPolys = NULL;
     psS32 maxChebyPoly = 0;
@@ -667,13 +667,13 @@
 }
 
-static double dOrdPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
+static psF64 dOrdPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly)
 {
     psS32 loop_x = 0;
     psS32 loop_y = 0;
     psS32 loop_z = 0;
-    double polySum = 0.0;
-    double xSum = 1.0;
-    double ySum = 1.0;
-    double zSum = 1.0;
+    psF64 polySum = 0.0;
+    psF64 xSum = 1.0;
+    psF64 ySum = 1.0;
+    psF64 zSum = 1.0;
 
     for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
@@ -695,5 +695,5 @@
 }
 
-static double dChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
+static psF64 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly)
 {
     psS32 loop_x = 0;
@@ -701,5 +701,5 @@
     psS32 loop_z = 0;
     psS32 i = 0;
-    double polySum = 0.0;
+    psF64 polySum = 0.0;
     psPolynomial1D* *chebPolys = NULL;
     psS32 maxChebyPoly = 0;
@@ -736,5 +736,5 @@
 }
 
-static double dOrdPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
+static psF64 dOrdPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly)
 {
     psS32 loop_w = 0;
@@ -742,9 +742,9 @@
     psS32 loop_y = 0;
     psS32 loop_z = 0;
-    double polySum = 0.0;
-    double wSum = 1.0;
-    double xSum = 1.0;
-    double ySum = 1.0;
-    double zSum = 1.0;
+    psF64 polySum = 0.0;
+    psF64 wSum = 1.0;
+    psF64 xSum = 1.0;
+    psF64 ySum = 1.0;
+    psF64 zSum = 1.0;
 
     for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
@@ -770,5 +770,5 @@
 }
 
-static double dChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
+static psF64 dChebPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly)
 {
     psS32 loop_w = 0;
@@ -777,5 +777,5 @@
     psS32 loop_z = 0;
     psS32 i = 0;
-    double polySum = 0.0;
+    psF64 polySum = 0.0;
     psPolynomial1D* *chebPolys = NULL;
     psS32 maxChebyPoly = 0;
@@ -828,5 +828,5 @@
  *****************************************************************************/
 #define FUNC_MACRO_FULL_INTERPOLATE_1D(TYPE) \
-static float fullInterpolate1D##TYPE(ps##TYPE *domain, \
+static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \
                                      ps##TYPE *range, \
                                      psS32 n, \
@@ -896,9 +896,9 @@
 LaGrange interpolation.
  *****************************************************************************/
-static float interpolate1DF32(float *domain,
-                              float *range,
+static psF32 interpolate1DF32(psF32 *domain,
+                              psF32 *range,
                               psS32 n,
                               psS32 order,
-                              float x)
+                              psF32 x)
 {
     psS32 binNum;
@@ -941,7 +941,7 @@
     evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f]
  *****************************************************************************/
-float psGaussian(float x, float mean, float sigma, psBool normal)
-{
-    float tmp = 1.0;
+psF32 psGaussian(psF32 x, psF32 mean, psF32 sigma, psBool normal)
+{
+    psF32 tmp = 1.0;
 
     psTrace(".psLib.dataManip.psFunctions.psGaussian", 4,
@@ -967,5 +967,5 @@
 XXX: There is no way to seed the random generator.
  *****************************************************************************/
-psVector* p_psGaussianDev(float mean, float sigma, psS32 Npts)
+psVector* p_psGaussianDev(psF32 mean, psF32 sigma, psS32 Npts)
 {
     PS_INT_CHECK_NON_NEGATIVE(Npts, NULL);
@@ -1007,7 +1007,7 @@
     newPoly->type = type;
     newPoly->n = n;
-    newPoly->coeff = (float *)psAlloc(n * sizeof(float));
-    newPoly->coeffErr = (float *)psAlloc(n * sizeof(float));
-    newPoly->mask = (char *)psAlloc(n * sizeof(char));
+    newPoly->coeff = (psF32 *)psAlloc(n * sizeof(psF32));
+    newPoly->coeffErr = (psF32 *)psAlloc(n * sizeof(psF32));
+    newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8));
     for (i = 0; i < n; i++) {
         newPoly->coeff[i] = 0.0;
@@ -1036,11 +1036,11 @@
     newPoly->nY = nY;
 
-    newPoly->coeff = (float **)psAlloc(nX * sizeof(float *));
-    newPoly->coeffErr = (float **)psAlloc(nX * sizeof(float *));
-    newPoly->mask = (char **)psAlloc(nX * sizeof(char *));
+    newPoly->coeff = (psF32 **)psAlloc(nX * sizeof(psF32 *));
+    newPoly->coeffErr = (psF32 **)psAlloc(nX * sizeof(psF32 *));
+    newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *));
     for (x = 0; x < nX; x++) {
-        newPoly->coeff[x] = (float *)psAlloc(nY * sizeof(float));
-        newPoly->coeffErr[x] = (float *)psAlloc(nY * sizeof(float));
-        newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char));
+        newPoly->coeff[x] = (psF32 *)psAlloc(nY * sizeof(psF32));
+        newPoly->coeffErr[x] = (psF32 *)psAlloc(nY * sizeof(psF32));
+        newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8));
     }
     for (x = 0; x < nX; x++) {
@@ -1075,15 +1075,15 @@
     newPoly->nZ = nZ;
 
-    newPoly->coeff = (float ***)psAlloc(nX * sizeof(float **));
-    newPoly->coeffErr = (float ***)psAlloc(nX * sizeof(float **));
-    newPoly->mask = (char ***)psAlloc(nX * sizeof(char **));
+    newPoly->coeff = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
+    newPoly->coeffErr = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
+    newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **));
     for (x = 0; x < nX; x++) {
-        newPoly->coeff[x] = (float **)psAlloc(nY * sizeof(float *));
-        newPoly->coeffErr[x] = (float **)psAlloc(nY * sizeof(float *));
-        newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *));
+        newPoly->coeff[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
+        newPoly->coeffErr[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
+        newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *));
         for (y = 0; y < nY; y++) {
-            newPoly->coeff[x][y] = (float *)psAlloc(nZ * sizeof(float));
-            newPoly->coeffErr[x][y] = (float *)psAlloc(nZ * sizeof(float));
-            newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char));
+            newPoly->coeff[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
+            newPoly->coeffErr[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
+            newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8));
         }
     }
@@ -1124,19 +1124,19 @@
     newPoly->nZ = nZ;
 
-    newPoly->coeff = (float ****)psAlloc(nW * sizeof(float ***));
-    newPoly->coeffErr = (float ****)psAlloc(nW * sizeof(float ***));
-    newPoly->mask = (char ****)psAlloc(nW * sizeof(char ***));
+    newPoly->coeff = (psF32 ****)psAlloc(nW * sizeof(psF32 ***));
+    newPoly->coeffErr = (psF32 ****)psAlloc(nW * sizeof(psF32 ***));
+    newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***));
     for (w = 0; w < nW; w++) {
-        newPoly->coeff[w] = (float ***)psAlloc(nX * sizeof(float **));
-        newPoly->coeffErr[w] = (float ***)psAlloc(nX * sizeof(float **));
-        newPoly->mask[w] = (char ***)psAlloc(nX * sizeof(char **));
+        newPoly->coeff[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
+        newPoly->coeffErr[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
+        newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **));
         for (x = 0; x < nX; x++) {
-            newPoly->coeff[w][x] = (float **)psAlloc(nY * sizeof(float *));
-            newPoly->coeffErr[w][x] = (float **)psAlloc(nY * sizeof(float *));
-            newPoly->mask[w][x] = (char **)psAlloc(nY * sizeof(char *));
+            newPoly->coeff[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
+            newPoly->coeffErr[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
+            newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *));
             for (y = 0; y < nY; y++) {
-                newPoly->coeff[w][x][y] = (float *)psAlloc(nZ * sizeof(float));
-                newPoly->coeffErr[w][x][y] = (float *)psAlloc(nZ * sizeof(float));
-                newPoly->mask[w][x][y] = (char *)psAlloc(nZ * sizeof(char));
+                newPoly->coeff[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
+                newPoly->coeffErr[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
+                newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8));
             }
         }
@@ -1157,5 +1157,5 @@
 }
 
-float psPolynomial1DEval(const psPolynomial1D* myPoly, float x)
+psF32 psPolynomial1DEval(const psPolynomial1D* myPoly, psF32 x)
 {
     PS_POLY_CHECK_NULL(myPoly, NAN);
@@ -1195,5 +1195,5 @@
 }
 
-float psPolynomial2DEval(const psPolynomial2D* myPoly, float x, float y)
+psF32 psPolynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y)
 {
     PS_POLY_CHECK_NULL(myPoly, NAN);
@@ -1250,5 +1250,5 @@
 }
 
-float psPolynomial3DEval(const psPolynomial3D* myPoly, float x, float y, float z)
+psF32 psPolynomial3DEval(const psPolynomial3D* myPoly, psF32 x, psF32 y, psF32 z)
 {
     PS_POLY_CHECK_NULL(myPoly, NAN);
@@ -1316,5 +1316,5 @@
 }
 
-float psPolynomial4DEval(const psPolynomial4D* myPoly, float w, float x, float y, float z)
+psF32 psPolynomial4DEval(const psPolynomial4D* myPoly, psF32 w, psF32 x, psF32 y, psF32 z)
 {
     PS_POLY_CHECK_NULL(myPoly, NAN);
@@ -1406,7 +1406,7 @@
     newPoly->type = type;
     newPoly->n = n;
-    newPoly->coeff = (double *)psAlloc(n * sizeof(double));
-    newPoly->coeffErr = (double *)psAlloc(n * sizeof(double));
-    newPoly->mask = (char *)psAlloc(n * sizeof(char));
+    newPoly->coeff = (psF64 *)psAlloc(n * sizeof(psF64));
+    newPoly->coeffErr = (psF64 *)psAlloc(n * sizeof(psF64));
+    newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8));
     for (i = 0; i < n; i++) {
         newPoly->coeff[i] = 0.0;
@@ -1435,11 +1435,11 @@
     newPoly->nY = nY;
 
-    newPoly->coeff = (double **)psAlloc(nX * sizeof(double *));
-    newPoly->coeffErr = (double **)psAlloc(nX * sizeof(double *));
-    newPoly->mask = (char **)psAlloc(nX * sizeof(char *));
+    newPoly->coeff = (psF64 **)psAlloc(nX * sizeof(psF64 *));
+    newPoly->coeffErr = (psF64 **)psAlloc(nX * sizeof(psF64 *));
+    newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *));
     for (x = 0; x < nX; x++) {
-        newPoly->coeff[x] = (double *)psAlloc(nY * sizeof(double));
-        newPoly->coeffErr[x] = (double *)psAlloc(nY * sizeof(double));
-        newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char));
+        newPoly->coeff[x] = (psF64 *)psAlloc(nY * sizeof(psF64));
+        newPoly->coeffErr[x] = (psF64 *)psAlloc(nY * sizeof(psF64));
+        newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8));
     }
     for (x = 0; x < nX; x++) {
@@ -1474,15 +1474,15 @@
     newPoly->nZ = nZ;
 
-    newPoly->coeff = (double ***)psAlloc(nX * sizeof(double **));
-    newPoly->coeffErr = (double ***)psAlloc(nX * sizeof(double **));
-    newPoly->mask = (char ***)psAlloc(nX * sizeof(char **));
+    newPoly->coeff = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
+    newPoly->coeffErr = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
+    newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **));
     for (x = 0; x < nX; x++) {
-        newPoly->coeff[x] = (double **)psAlloc(nY * sizeof(double *));
-        newPoly->coeffErr[x] = (double **)psAlloc(nY * sizeof(double *));
-        newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *));
+        newPoly->coeff[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
+        newPoly->coeffErr[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
+        newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *));
         for (y = 0; y < nY; y++) {
-            newPoly->coeff[x][y] = (double *)psAlloc(nZ * sizeof(double));
-            newPoly->coeffErr[x][y] = (double *)psAlloc(nZ * sizeof(double));
-            newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char));
+            newPoly->coeff[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
+            newPoly->coeffErr[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
+            newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8));
         }
     }
@@ -1523,19 +1523,19 @@
     newPoly->nZ = nZ;
 
-    newPoly->coeff = (double ****)psAlloc(nW * sizeof(double ***));
-    newPoly->coeffErr = (double ****)psAlloc(nW * sizeof(double ***));
-    newPoly->mask = (char ****)psAlloc(nW * sizeof(char ***));
+    newPoly->coeff = (psF64 ****)psAlloc(nW * sizeof(psF64 ***));
+    newPoly->coeffErr = (psF64 ****)psAlloc(nW * sizeof(psF64 ***));
+    newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***));
     for (w = 0; w < nW; w++) {
-        newPoly->coeff[w] = (double ***)psAlloc(nX * sizeof(double **));
-        newPoly->coeffErr[w] = (double ***)psAlloc(nX * sizeof(double **));
-        newPoly->mask[w] = (char ***)psAlloc(nX * sizeof(char **));
+        newPoly->coeff[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
+        newPoly->coeffErr[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
+        newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **));
         for (x = 0; x < nX; x++) {
-            newPoly->coeff[w][x] = (double **)psAlloc(nY * sizeof(double *));
-            newPoly->coeffErr[w][x] = (double **)psAlloc(nY * sizeof(double *));
-            newPoly->mask[w][x] = (char **)psAlloc(nY * sizeof(char *));
+            newPoly->coeff[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
+            newPoly->coeffErr[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
+            newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *));
             for (y = 0; y < nY; y++) {
-                newPoly->coeff[w][x][y] = (double *)psAlloc(nZ * sizeof(double));
-                newPoly->coeffErr[w][x][y] = (double *)psAlloc(nZ * sizeof(double));
-                newPoly->mask[w][x][y] = (char *)psAlloc(nZ * sizeof(char));
+                newPoly->coeff[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
+                newPoly->coeffErr[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
+                newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8));
             }
         }
@@ -1557,5 +1557,5 @@
 
 
-double psDPolynomial1DEval(const psDPolynomial1D* myPoly, double x)
+psF64 psDPolynomial1DEval(const psDPolynomial1D* myPoly, psF64 x)
 {
     PS_POLY_CHECK_NULL(myPoly, NAN);
@@ -1598,7 +1598,7 @@
 
 
-double psDPolynomial2DEval(const psDPolynomial2D* myPoly,
-                           double x,
-                           double y)
+psF64 psDPolynomial2DEval(const psDPolynomial2D* myPoly,
+                          psF64 x,
+                          psF64 y)
 {
     PS_POLY_CHECK_NULL(myPoly, NAN);
@@ -1655,8 +1655,8 @@
 
 
-double psDPolynomial3DEval(const psDPolynomial3D* myPoly,
-                           double x,
-                           double y,
-                           double z)
+psF64 psDPolynomial3DEval(const psDPolynomial3D* myPoly,
+                          psF64 x,
+                          psF64 y,
+                          psF64 z)
 {
     PS_POLY_CHECK_NULL(myPoly, NAN);
@@ -1724,9 +1724,9 @@
 }
 
-double psDPolynomial4DEval(const psDPolynomial4D* myPoly,
-                           double w,
-                           double x,
-                           double y,
-                           double z)
+psF64 psDPolynomial4DEval(const psDPolynomial4D* myPoly,
+                          psF64 w,
+                          psF64 x,
+                          psF64 y,
+                          psF64 z)
 {
     PS_POLY_CHECK_NULL(myPoly, NAN);
@@ -1810,6 +1810,6 @@
 //    psS32 n;
 //    psPolynomial1D **spline;
-//    float *p_psDeriv2;
-//    float *domains;
+//    psF32 *p_psDeriv2;
+//    psF32 *domains;
 //} psSpline1D;
 
@@ -1822,6 +1822,6 @@
 psSpline1D *psSpline1DAlloc(psS32 numSplines,
                             psS32 order,
-                            float min,
-                            float max)
+                            psF32 min,
+                            psF32 max)
 {
     PS_INT_CHECK_NON_NEGATIVE(numSplines, NULL);
@@ -1831,6 +1831,6 @@
     psSpline1D *tmp = NULL;
     psS32 i;
-    float tmpDomain;
-    float width;
+    psF32 tmpDomain;
+    psF32 width;
 
     tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D));
@@ -1845,6 +1845,6 @@
     tmp->p_psDeriv2 = NULL;
 
-    tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float));
-    width = (max - min) / ((float) numSplines);
+    tmp->domains = (psF32 *) psAlloc((numSplines+1) * sizeof(psF32));
+    width = (max - min) / ((psF32) numSplines);
 
     (tmp->domains)[0] = min;
@@ -1884,5 +1884,5 @@
     }
 
-    tmp->domains = (float *) psAlloc((bounds->n) * sizeof(float));
+    tmp->domains = (psF32 *) psAlloc((bounds->n) * sizeof(psF32));
 
     for (i=0;i<bounds->n;i++) {
@@ -1977,5 +1977,5 @@
 {
     PS_PTR_CHECK_TYPE_EQUAL(x, bins, -3);
-    char* strType;
+    psS8* strType;
 
     switch (x->type.type) {
@@ -2074,10 +2074,10 @@
         domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32);
 
-        psScalar *tmpScalar = psScalarAlloc((double)
+        psScalar *tmpScalar = psScalarAlloc((psF64)
                                             interpolate1DF32(domain32->data.F32,
                                                              range32->data.F32,
                                                              domain32->n,
                                                              order,
-                                                             (float) x->data.F64), PS_TYPE_F64);
+                                                             (psF32) x->data.F64), PS_TYPE_F64);
         psFree(range32);
         psFree(domain32);
@@ -2089,5 +2089,5 @@
 
     } else {
-        char* strType;
+        psS8* strType;
         PS_TYPE_NAME(strType,x->type.type);
         psError(PS_ERR_BAD_PARAMETER_TYPE,
@@ -2115,6 +2115,6 @@
      the spline fit functions require F32 and F64.
  *****************************************************************************/
-float psSpline1DEval(
-    float x,
+psF32 psSpline1DEval(
+    psF32 x,
     const psSpline1D *spline
 )
@@ -2172,10 +2172,10 @@
         for (i=0;i<x->n;i++) {
             tmpVector->data.F32[i] = psSpline1DEval(
-                                         (float) x->data.F64[i],
+                                         (psF32) x->data.F64[i],
                                          spline
                                      );
         }
     } else {
-        char* strType;
+        psS8* strType;
         PS_TYPE_NAME(strType,x->type.type);
         psError(PS_ERR_BAD_PARAMETER_TYPE,