Changeset 2273 for trunk/psLib/src/math/psSpline.c

Timestamp:

Nov 3, 2004, 3:05:00 PM (22 years ago)

Author:

desonia

Message:

changed the psError signature to match SDRS. Also made misc. cleanups as
I was combing the files.

File:

• trunk/psLib/src/math/psSpline.c (modified) (40 diffs)

trunk/psLib/src/math/psSpline.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.58 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-28 00:22:53 $
+ *  @version $Revision: 1.59 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-11-04 01:04:59 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -35 +35 @@
 #include "psFunctions.h"
 #include "psConstants.h"
+
+#include "psDataManipErrors.h"
+
 /*****************************************************************************/
 /* DEFINE STATEMENTS                                                         */
@@ -50 +53 @@
 static void dPolynomial3DFree(psDPolynomial3D* myPoly);
 static void dPolynomial4DFree(psDPolynomial4D* myPoly);
+static void spline1DFree(psSpline1D *tmpSpline);
+static psS32 vectorBinDisectF32(float *bins,psS32 numBins,float x);
+static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x);
 
 /*****************************************************************************/
@@ -67 +73 @@
 /*****************************************************************************/
 
+static void spline1DFree(psSpline1D *tmpSpline)
+{
+    psS32 i;
+
+    if (tmpSpline == NULL) {
+        return;
+    }
+
+    if (tmpSpline->spline != NULL) {
+        for (i=0;i<tmpSpline->n;i++) {
+            psFree((tmpSpline->spline)[i]);
+        }
+        psFree(tmpSpline->spline);
+    }
+
+    if (tmpSpline->p_psDeriv2 != NULL) {
+        psFree(tmpSpline->p_psDeriv2);
+    }
+    psFree(tmpSpline->domains);
+
+    return;
+}
+
+static void polynomial1DFree(psPolynomial1D* myPoly)
+{
+    psFree(myPoly->coeff);
+    psFree(myPoly->coeffErr);
+    psFree(myPoly->mask);
+}
+
+static void polynomial2DFree(psPolynomial2D* myPoly)
+{
+    psS32 x = 0;
+
+    for (x = 0; x < myPoly->nX; x++) {
+        psFree(myPoly->coeff[x]);
+        psFree(myPoly->coeffErr[x]);
+        psFree(myPoly->mask[x]);
+    }
+    psFree(myPoly->coeff);
+    psFree(myPoly->coeffErr);
+    psFree(myPoly->mask);
+}
+
+static void polynomial3DFree(psPolynomial3D* myPoly)
+{
+    psS32 x = 0;
+    psS32 y = 0;
+
+    for (x = 0; x < myPoly->nX; x++) {
+        for (y = 0; y < myPoly->nY; y++) {
+            psFree(myPoly->coeff[x][y]);
+            psFree(myPoly->coeffErr[x][y]);
+            psFree(myPoly->mask[x][y]);
+        }
+        psFree(myPoly->coeff[x]);
+        psFree(myPoly->coeffErr[x]);
+        psFree(myPoly->mask[x]);
+    }
+
+    psFree(myPoly->coeff);
+    psFree(myPoly->coeffErr);
+    psFree(myPoly->mask);
+}
+
+static void polynomial4DFree(psPolynomial4D* myPoly)
+{
+    psS32 w = 0;
+    psS32 x = 0;
+    psS32 y = 0;
+
+    for (w = 0; w < myPoly->nW; w++) {
+        for (x = 0; x < myPoly->nX; x++) {
+            for (y = 0; y < myPoly->nY; y++) {
+                psFree(myPoly->coeff[w][x][y]);
+                psFree(myPoly->coeffErr[w][x][y]);
+                psFree(myPoly->mask[w][x][y]);
+            }
+            psFree(myPoly->coeff[w][x]);
+            psFree(myPoly->coeffErr[w][x]);
+            psFree(myPoly->mask[w][x]);
+        }
+        psFree(myPoly->coeff[w]);
+        psFree(myPoly->coeffErr[w]);
+        psFree(myPoly->mask[w]);
+    }
+
+    psFree(myPoly->coeff);
+    psFree(myPoly->coeffErr);
+    psFree(myPoly->mask);
+}
+
+static void dPolynomial1DFree(psDPolynomial1D* myPoly)
+{
+    psFree(myPoly->coeff);
+    psFree(myPoly->coeffErr);
+    psFree(myPoly->mask);
+}
+
+static void dPolynomial2DFree(psDPolynomial2D* myPoly)
+{
+    psS32 x = 0;
+
+    for (x = 0; x < myPoly->nX; x++) {
+        psFree(myPoly->coeff[x]);
+        psFree(myPoly->coeffErr[x]);
+        psFree(myPoly->mask[x]);
+    }
+    psFree(myPoly->coeff);
+    psFree(myPoly->coeffErr);
+    psFree(myPoly->mask);
+}
+
+static void dPolynomial3DFree(psDPolynomial3D* myPoly)
+{
+    psS32 x = 0;
+    psS32 y = 0;
+
+    for (x = 0; x < myPoly->nX; x++) {
+        for (y = 0; y < myPoly->nY; y++) {
+            psFree(myPoly->coeff[x][y]);
+            psFree(myPoly->coeffErr[x][y]);
+            psFree(myPoly->mask[x][y]);
+        }
+        psFree(myPoly->coeff[x]);
+        psFree(myPoly->coeffErr[x]);
+        psFree(myPoly->mask[x]);
+    }
+
+    psFree(myPoly->coeff);
+    psFree(myPoly->coeffErr);
+    psFree(myPoly->mask);
+}
+
+static void dPolynomial4DFree(psDPolynomial4D* myPoly)
+{
+    psS32 w = 0;
+    psS32 x = 0;
+    psS32 y = 0;
+
+    for (w = 0; w < myPoly->nW; w++) {
+        for (x = 0; x < myPoly->nX; x++) {
+            for (y = 0; y < myPoly->nY; y++) {
+                psFree(myPoly->coeff[w][x][y]);
+                psFree(myPoly->coeffErr[w][x][y]);
+                psFree(myPoly->mask[w][x][y]);
+            }
+            psFree(myPoly->coeff[w][x]);
+            psFree(myPoly->coeffErr[w][x]);
+            psFree(myPoly->mask[w][x]);
+        }
+        psFree(myPoly->coeff[w]);
+        psFree(myPoly->coeffErr[w]);
+        psFree(myPoly->mask[w]);
+    }
+
+    psFree(myPoly->coeff);
+    psFree(myPoly->coeffErr);
+    psFree(myPoly->mask);
+}
+
 /*****************************************************************************
-CreateChebyshevPolys(n): this routine takes as input the required order n,
+createChebyshevPolys(n): this routine takes as input the required order n,
 and returns as output as a pointer to an array of n psPolynomial1D
 structures, corresponding to the first n Chebyshev polynomials.
@@ -76 +243 @@
 outer coefficients of the Chebyshev polynomials.
  *****************************************************************************/
-static psPolynomial1D **CreateChebyshevPolys(psS32 maxChebyPoly)
+static psPolynomial1D **createChebyshevPolys(psS32 maxChebyPoly)
 {
     PS_INT_CHECK_NON_NEGATIVE(maxChebyPoly, NULL);
@@ -103 +270 @@
 
     return (chebPolys);
+}
+
+/*****************************************************************************
+    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
+ 
+    XXX: Should the "coeffErr[]" should be used as well?
+ *****************************************************************************/
+static float ordPolynomial1DEval(float x, const psPolynomial1D* myPoly)
+{
+    psS32 loop_x = 0;
+    float polySum = 0.0;
+    float xSum = 1.0;
+
+    psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4,
+            "---- Calling ordPolynomial1DEval(%f)\n", x);
+    psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4,
+            "Polynomial order is %d\n", myPoly->n);
+    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
+        psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4,
+                "Polynomial coeff[%d] is %f\n", loop_x, myPoly->coeff[loop_x]);
+    }
+
+    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
+        if (myPoly->mask[loop_x] == 0) {
+            psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 10,
+                    "polysum+= sum*coeff [%f+= (%f * %f)\n", polySum, xSum, myPoly->coeff[loop_x]);
+            polySum += xSum * myPoly->coeff[loop_x];
+            xSum *= x;
+        }
+    }
+
+    return(polySum);
+}
+
+// XXX: You can do this without having to psAlloc() vector d.
+// XXX: How does the mask vector effect Crenshaw's formula?
+static float chebPolynomial1DEval(float x, const psPolynomial1D* myPoly)
+{
+    psVector *d;
+    psS32 n;
+    psS32 i;
+    float tmp;
+
+    n = myPoly->n;
+    d = psVectorAlloc(n, PS_TYPE_F32);
+    d->data.F32[n-1] = myPoly->coeff[n-1];
+    d->data.F32[n-2] = (2.0 * x * d->data.F32[n-1]) + myPoly->coeff[n-2];
+    for (i=n-3;i>=1;i--) {
+        d->data.F32[i] = (2.0 * x * d->data.F32[i+1]) -
+                         (d->data.F32[i+2]) +
+                         (myPoly->coeff[i]);
+    }
+
+    tmp = (x * d->data.F32[1]) -
+          (d->data.F32[2]) +
+          (0.5 * myPoly->coeff[0]);
+
+    psFree(d);
+    return(tmp);
+    /*
+
+    psS32 n;
+    psS32 i;
+    float tmp;
+    psPolynomial1D **chebPolys = NULL;
+
+    n = myPoly->n;
+    chebPolys = createChebyshevPolys(n);
+
+    tmp = 0.0;
+    for (i=0;i<myPoly->n;i++) {
+        tmp+= (myPoly->coeff[i] * psPolynomial1DEval(x, chebPolys[i]));
+        //            printf("HMMM: psPolynomial1DEval(%f, chebPolys[%d]) is %f\n", x, i, psPolynomial1DEval(x, chebPolys[i]));
+    }
+    tmp-= (myPoly->coeff[0]/2.0);
+
+
+    return(tmp);
+    */
+}
+
+static float ordPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
+{
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
+    psS32 loop_x = 0;
+    psS32 loop_y = 0;
+    float polySum = 0.0;
+    float xSum = 1.0;
+    float ySum = 1.0;
+
+    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+        ySum = xSum;
+        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+            if (myPoly->mask[loop_x][loop_y] == 0) {
+                polySum += ySum * myPoly->coeff[loop_x][loop_y];
+                ySum *= y;
+            }
+        }
+        xSum *= x;
+    }
+
+    return(polySum);
+}
+
+static float chebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
+{
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
+    psS32 loop_x = 0;
+    psS32 loop_y = 0;
+    psS32 i = 0;
+    float polySum = 0.0;
+    psPolynomial1D* *chebPolys = NULL;
+    psS32 maxChebyPoly = 0;
+
+    // Determine how many Chebyshev polynomials
+    // are needed, then create them.
+    maxChebyPoly = myPoly->nX;
+    if (myPoly->nY > maxChebyPoly) {
+        maxChebyPoly = myPoly->nY;
+    }
+    chebPolys = createChebyshevPolys(maxChebyPoly);
+
+    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+            if (myPoly->mask[loop_x][loop_y] == 0) {
+                polySum += myPoly->coeff[loop_x][loop_y] *
+                           psPolynomial1DEval(x, chebPolys[loop_x]) *
+                           psPolynomial1DEval(y, chebPolys[loop_y]);
+            }
+        }
+    }
+    for (i=0;i<maxChebyPoly;i++) {
+        psFree(chebPolys[i]);
+    }
+    psFree(chebPolys);
+    return(polySum);
+}
+
+static float ordPolynomial3DEval(float x, float y, float 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;
+
+    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+        ySum = xSum;
+        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+            zSum = ySum;
+            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
+                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
+                    polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z];
+                    zSum *= z;
+                }
+            }
+            ySum *= y;
+        }
+        xSum *= x;
+    }
+
+    return(polySum);
+}
+
+static float chebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
+{
+    psS32 loop_x = 0;
+    psS32 loop_y = 0;
+    psS32 loop_z = 0;
+    psS32 i = 0;
+    float polySum = 0.0;
+    psPolynomial1D* *chebPolys = NULL;
+    psS32 maxChebyPoly = 0;
+
+    // Determine how many Chebyshev polynomials
+    // are needed, then create them.
+    maxChebyPoly = myPoly->nX;
+    if (myPoly->nY > maxChebyPoly) {
+        maxChebyPoly = myPoly->nY;
+    }
+    if (myPoly->nZ > maxChebyPoly) {
+        maxChebyPoly = myPoly->nZ;
+    }
+    chebPolys = createChebyshevPolys(maxChebyPoly);
+
+    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
+                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
+                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
+                               psPolynomial1DEval(x, chebPolys[loop_x]) *
+                               psPolynomial1DEval(y, chebPolys[loop_y]) *
+                               psPolynomial1DEval(z, chebPolys[loop_z]);
+                }
+            }
+        }
+    }
+
+    for (i=0;i<maxChebyPoly;i++) {
+        psFree(chebPolys[i]);
+    }
+    psFree(chebPolys);
+    return(polySum);
+}
+
+static float ordPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
+{
+    psS32 loop_w = 0;
+    psS32 loop_x = 0;
+    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;
+
+    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
+        xSum = wSum;
+        for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+            ySum = xSum;
+            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+                zSum = ySum;
+                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
+                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
+                        polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                        zSum *= z;
+                    }
+                }
+                ySum *= y;
+            }
+            xSum *= x;
+        }
+        wSum *= w;
+    }
+
+    return(polySum);
+}
+
+static float chebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
+{
+    psS32 loop_w = 0;
+    psS32 loop_x = 0;
+    psS32 loop_y = 0;
+    psS32 loop_z = 0;
+    psS32 i = 0;
+    float polySum = 0.0;
+    psPolynomial1D* *chebPolys = NULL;
+    psS32 maxChebyPoly = 0;
+
+    // Determine how many Chebyshev polynomials
+    // are needed, then create them.
+    maxChebyPoly = myPoly->nW;
+    if (myPoly->nX > maxChebyPoly) {
+        maxChebyPoly = myPoly->nX;
+    }
+    if (myPoly->nY > maxChebyPoly) {
+        maxChebyPoly = myPoly->nY;
+    }
+    if (myPoly->nZ > maxChebyPoly) {
+        maxChebyPoly = myPoly->nZ;
+    }
+    chebPolys = createChebyshevPolys(maxChebyPoly);
+
+    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
+        for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
+                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
+                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
+                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
+                                   psPolynomial1DEval(x, chebPolys[loop_x]) *
+                                   psPolynomial1DEval(y, chebPolys[loop_y]) *
+                                   psPolynomial1DEval(z, chebPolys[loop_z]);
+                    }
+                }
+            }
+        }
+    }
+
+    for (i=0;i<maxChebyPoly;i++) {
+        psFree(chebPolys[i]);
+    }
+    psFree(chebPolys);
+    return(polySum);
+}
+
+/*****************************************************************************
+    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
+ *****************************************************************************/
+static double dOrdPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
+{
+    psS32 loop_x = 0;
+    double polySum = 0.0;
+    double xSum = 1.0;
+
+    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
+        if (myPoly->mask[loop_x] == 0) {
+            polySum += xSum * myPoly->coeff[loop_x];
+            xSum *= x;
+        }
+    }
+
+    return(polySum);
+}
+
+// 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)
+{
+    psVector *d;
+    psS32 n;
+    psS32 i;
+    double tmp;
+
+    n = myPoly->n;
+    d = psVectorAlloc(n, PS_TYPE_F64);
+    d->data.F64[n-1] = myPoly->coeff[n-1];
+    d->data.F64[n-2] = (2.0 * x * d->data.F64[n-1]) + myPoly->coeff[n-2];
+    for (i=n-3;i>=1;i--) {
+        d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) -
+                         (d->data.F64[i+2]) +
+                         (myPoly->coeff[i]);
+    }
+
+    tmp = (x * d->data.F64[1]) -
+          (d->data.F64[2]) +
+          (0.5 * myPoly->coeff[0]);
+
+    psFree(d);
+    return(tmp);
+}
+
+static double dOrdPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
+{
+    psS32 loop_x = 0;
+    psS32 loop_y = 0;
+    double polySum = 0.0;
+    double xSum = 1.0;
+    double ySum = 1.0;
+
+    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+        ySum = xSum;
+        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+            if (myPoly->mask[loop_x][loop_y] == 0) {
+                polySum += ySum * myPoly->coeff[loop_x][loop_y];
+                ySum *= y;
+            }
+        }
+        xSum *= x;
+    }
+
+    return(polySum);
+}
+
+static double dChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
+{
+    psS32 loop_x = 0;
+    psS32 loop_y = 0;
+    psS32 i = 0;
+    double polySum = 0.0;
+    psPolynomial1D* *chebPolys = NULL;
+    psS32 maxChebyPoly = 0;
+
+    // Determine how many Chebyshev polynomials
+    // are needed, then create them.
+    maxChebyPoly = myPoly->nX;
+    if (myPoly->nY > maxChebyPoly) {
+        maxChebyPoly = myPoly->nY;
+    }
+    chebPolys = createChebyshevPolys(maxChebyPoly);
+
+    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+            if (myPoly->mask[loop_x][loop_y] == 0) {
+                polySum += myPoly->coeff[loop_x][loop_y] *
+                           psPolynomial1DEval(x, chebPolys[loop_x]) *
+                           psPolynomial1DEval(y, chebPolys[loop_y]);
+            }
+        }
+    }
+
+    for (i=0;i<maxChebyPoly;i++) {
+        psFree(chebPolys[i]);
+    }
+    psFree(chebPolys);
+    return(polySum);
+}
+
+static double dOrdPolynomial3DEval(double x, double y, double 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;
+
+    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+        ySum = xSum;
+        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+            zSum = ySum;
+            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
+                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
+                    polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z];
+                    zSum *= z;
+                }
+            }
+            ySum *= y;
+        }
+        xSum *= x;
+    }
+
+    return(polySum);
+}
+
+static double dChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
+{
+    psS32 loop_x = 0;
+    psS32 loop_y = 0;
+    psS32 loop_z = 0;
+    psS32 i = 0;
+    double polySum = 0.0;
+    psPolynomial1D* *chebPolys = NULL;
+    psS32 maxChebyPoly = 0;
+
+    // Determine how many Chebyshev polynomials
+    // are needed, then create them.
+    maxChebyPoly = myPoly->nX;
+    if (myPoly->nY > maxChebyPoly) {
+        maxChebyPoly = myPoly->nY;
+    }
+    if (myPoly->nZ > maxChebyPoly) {
+        maxChebyPoly = myPoly->nZ;
+    }
+    chebPolys = createChebyshevPolys(maxChebyPoly);
+
+    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
+                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
+                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
+                               psPolynomial1DEval(x, chebPolys[loop_x]) *
+                               psPolynomial1DEval(y, chebPolys[loop_y]) *
+                               psPolynomial1DEval(z, chebPolys[loop_z]);
+                }
+            }
+        }
+    }
+
+    for (i=0;i<maxChebyPoly;i++) {
+        psFree(chebPolys[i]);
+    }
+    psFree(chebPolys);
+    return(polySum);
+}
+
+static double dOrdPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
+{
+    psS32 loop_w = 0;
+    psS32 loop_x = 0;
+    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;
+
+    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
+        xSum = wSum;
+        for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+            ySum = xSum;
+            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+                zSum = ySum;
+                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
+                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
+                        polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                        zSum *= z;
+                    }
+                }
+                ySum *= y;
+            }
+            xSum *= x;
+        }
+        wSum *= w;
+    }
+
+    return(polySum);
+}
+
+static double dChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
+{
+    psS32 loop_w = 0;
+    psS32 loop_x = 0;
+    psS32 loop_y = 0;
+    psS32 loop_z = 0;
+    psS32 i = 0;
+    double polySum = 0.0;
+    psPolynomial1D* *chebPolys = NULL;
+    psS32 maxChebyPoly = 0;
+
+    // Determine how many Chebyshev polynomials
+    // are needed, then create them.
+    maxChebyPoly = myPoly->nW;
+    if (myPoly->nX > maxChebyPoly) {
+        maxChebyPoly = myPoly->nX;
+    }
+    if (myPoly->nY > maxChebyPoly) {
+        maxChebyPoly = myPoly->nY;
+    }
+    if (myPoly->nZ > maxChebyPoly) {
+        maxChebyPoly = myPoly->nZ;
+    }
+    chebPolys = createChebyshevPolys(maxChebyPoly);
+
+    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
+        for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
+            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
+                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
+                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
+                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
+                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
+                                   psPolynomial1DEval(x, chebPolys[loop_x]) *
+                                   psPolynomial1DEval(y, chebPolys[loop_y]) *
+                                   psPolynomial1DEval(z, chebPolys[loop_z]);
+                    }
+                }
+            }
+        }
+    }
+
+    for (i=0;i<maxChebyPoly;i++) {
+        psFree(chebPolys[i]);
+    }
+    psFree(chebPolys);
+    return(polySum);
+}
+
+
+/*****************************************************************************
+p_psInterpolate1D(): This routine will take as input n-element floating
+point arrays domain and range, and the x value, assumed to lie with the
+domain vector.  It produces as output the (n-1)-order LaGrange interpolated
+value of x.
+ 
+XXX: do we error check for non-distinct domain values?
+ *****************************************************************************/
+static float fullInterpolate1DF32(float *domain,
+                                  float *range,
+                                  psS32 n,
+                                  float x)
+{
+    PS_INT_CHECK_NON_NEGATIVE(n, NAN);
+    PS_PTR_CHECK_NULL(domain, NAN);
+    PS_PTR_CHECK_NULL(range, NAN);
+
+    psS32 i;
+    psS32 m;
+    static psVector *p = NULL;
+    p = psVectorRecycle(p, n, PS_TYPE_F32);
+    p_psMemSetPersistent(p, true);
+    p_psMemSetPersistent(p->data.F32, true);
+    /*
+        psVector *p = psVectorAlloc(n, PS_TYPE_F32);
+        float tmp;
+    */
+
+    psTrace(".psLib.dataManip.psFunctions.fullInterpolate1DF32", 4,
+            "---- fullInterpolate1DF32() begin (%d-order at x=%f) (%d data points)----\n", n-1, x, n);
+
+    for (i=0;i<n;i++) {
+        psTrace(".psLib.dataManip.psFunctions.fullInterpolate1DF32", 6,
+                "domain/range is (%f %f)\n", domain[i], range[i]);
+    }
+
+    for (i=0;i<n;i++) {
+        p->data.F32[i] = range[i];
+        psTrace(".psLib.dataManip.psFunctions.fullInterpolate1DF32", 6,
+                "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
+
+    }
+
+    // From NR, during each iteration of the m loop, we are computing the
+    // p_{i ... i+m} terms.
+    for (m=1;m<n;m++) {
+        for (i=0;i<n-m;i++) {
+            // From NR: we are computing P_{i ... i+m}
+            p->data.F32[i] = (((x-domain[i+m]) * p->data.F32[i]) +
+                              ((domain[i]-x) * p->data.F32[i+1])) /
+                             (domain[i] - domain[i+m]);
+            //printf("((%f-%f * %f) + (%f-%f * %f)) / (%f - %f)\n", x, domain[i+m], p->data.F32[i], domain[i], x, p->data.F32[i+1], domain[i], domain[i+m]);
+            psTrace(".psLib.dataManip.psFunctions.fullInterpolate1DF32", 6,
+                    "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
+        }
+    }
+    psTrace(".psLib.dataManip.psFunctions.fullInterpolate1DF32", 4,
+            "---- fullInterpolate1DF32() end ----\n");
+
+    /*
+        tmp = p->data.F32[0];
+        psFree(p);
+        return(tmp);
+    */
+    return(p->data.F32[0]);
+}
+
+
+/*****************************************************************************
+interpolate1DF32(): this is the base 1-D flat memory routine to perform
+LaGrange interpolation.
+ *****************************************************************************/
+static float interpolate1DF32(float *domain,
+                              float *range,
+                              psS32 n,
+                              psS32 order,
+                              float x)
+{
+    psS32 binNum;
+    psS32 numIntPoints = order+1;
+    psS32 origin;
+
+    psTrace(".psLib.dataManip.psFunctions.interpolate1DF32", 4,
+            "---- interpolate1DF32() begin ----\n");
+
+    binNum = vectorBinDisectF32(domain, n, x);
+
+    if (0 == numIntPoints%2) {
+        origin = binNum - ((numIntPoints/2) - 1);
+    } else {
+        origin = binNum - (numIntPoints/2);
+        if ((x-domain[binNum]) > (domain[binNum+1]-x)) {
+            // x is closer to binNum+1.
+            origin = 1 + (binNum - (numIntPoints/2));
+        }
+    }
+    if (origin < 0) {
+        origin = 0;
+    }
+    if ((origin + numIntPoints) > n) {
+        origin = n - numIntPoints;
+    }
+
+    psTrace(".psLib.dataManip.psFunctions.interpolate1DF32", 4,
+            "---- interpolate1DF32() end ----\n");
+    return(fullInterpolate1DF32(&domain[origin], &range[origin], order+1, x));
 }
 
@@ -336 +1154 @@
 }
 
-static void polynomial1DFree(psPolynomial1D* myPoly)
-{
-    psFree(myPoly->coeff);
-    psFree(myPoly->coeffErr);
-    psFree(myPoly->mask);
-}
-
-static void polynomial2DFree(psPolynomial2D* myPoly)
-{
-    psS32 x = 0;
-
-    for (x = 0; x < myPoly->nX; x++) {
-        psFree(myPoly->coeff[x]);
-        psFree(myPoly->coeffErr[x]);
-        psFree(myPoly->mask[x]);
-    }
-    psFree(myPoly->coeff);
-    psFree(myPoly->coeffErr);
-    psFree(myPoly->mask);
-}
-
-static void polynomial3DFree(psPolynomial3D* myPoly)
-{
-    psS32 x = 0;
-    psS32 y = 0;
-
-    for (x = 0; x < myPoly->nX; x++) {
-        for (y = 0; y < myPoly->nY; y++) {
-            psFree(myPoly->coeff[x][y]);
-            psFree(myPoly->coeffErr[x][y]);
-            psFree(myPoly->mask[x][y]);
-        }
-        psFree(myPoly->coeff[x]);
-        psFree(myPoly->coeffErr[x]);
-        psFree(myPoly->mask[x]);
-    }
-
-    psFree(myPoly->coeff);
-    psFree(myPoly->coeffErr);
-    psFree(myPoly->mask);
-}
-
-static void polynomial4DFree(psPolynomial4D* myPoly)
-{
-    psS32 w = 0;
-    psS32 x = 0;
-    psS32 y = 0;
-
-    for (w = 0; w < myPoly->nW; w++) {
-        for (x = 0; x < myPoly->nX; x++) {
-            for (y = 0; y < myPoly->nY; y++) {
-                psFree(myPoly->coeff[w][x][y]);
-                psFree(myPoly->coeffErr[w][x][y]);
-                psFree(myPoly->mask[w][x][y]);
-            }
-            psFree(myPoly->coeff[w][x]);
-            psFree(myPoly->coeffErr[w][x]);
-            psFree(myPoly->mask[w][x]);
-        }
-        psFree(myPoly->coeff[w]);
-        psFree(myPoly->coeffErr[w]);
-        psFree(myPoly->mask[w]);
-    }
-
-    psFree(myPoly->coeff);
-    psFree(myPoly->coeffErr);
-    psFree(myPoly->mask);
-}
-
-/*****************************************************************************
-    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
- 
-    XXX: Should the "coeffErr[]" should be used as well?
- *****************************************************************************/
-float p_psOrdPolynomial1DEval(float x, const psPolynomial1D* myPoly)
-{
-    psS32 loop_x = 0;
-    float polySum = 0.0;
-    float xSum = 1.0;
-
-    psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4,
-            "---- Calling p_psOrdPolynomial1DEval(%f)\n", x);
-    psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4,
-            "Polynomial order is %d\n", myPoly->n);
-    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
-        psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 4,
-                "Polynomial coeff[%d] is %f\n", loop_x, myPoly->coeff[loop_x]);
-    }
-
-    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
-        if (myPoly->mask[loop_x] == 0) {
-            psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 10,
-                    "polysum+= sum*coeff [%f+= (%f * %f)\n", polySum, xSum, myPoly->coeff[loop_x]);
-            polySum += xSum * myPoly->coeff[loop_x];
-            xSum *= x;
-        }
-    }
-
-    return(polySum);
-}
-
-// XXX: You can do this without having to psAlloc() vector d.
-// XXX: How does the mask vector effect Crenshaw's formula?
-float p_psChebPolynomial1DEval(float x, const psPolynomial1D* myPoly)
-{
-    psVector *d;
-    psS32 n;
-    psS32 i;
-    float tmp;
-
-    n = myPoly->n;
-    d = psVectorAlloc(n, PS_TYPE_F32);
-    d->data.F32[n-1] = myPoly->coeff[n-1];
-    d->data.F32[n-2] = (2.0 * x * d->data.F32[n-1]) + myPoly->coeff[n-2];
-    for (i=n-3;i>=1;i--) {
-        d->data.F32[i] = (2.0 * x * d->data.F32[i+1]) -
-                         (d->data.F32[i+2]) +
-                         (myPoly->coeff[i]);
-    }
-
-    tmp = (x * d->data.F32[1]) -
-          (d->data.F32[2]) +
-          (0.5 * myPoly->coeff[0]);
-
-    psFree(d);
-    return(tmp);
-    /*
-
-    psS32 n;
-    psS32 i;
-    float tmp;
-    psPolynomial1D **chebPolys = NULL;
-
-    n = myPoly->n;
-    chebPolys = CreateChebyshevPolys(n);
-
-    tmp = 0.0;
-    for (i=0;i<myPoly->n;i++) {
-        tmp+= (myPoly->coeff[i] * psPolynomial1DEval(x, chebPolys[i]));
-        //            printf("HMMM: psPolynomial1DEval(%f, chebPolys[%d]) is %f\n", x, i, psPolynomial1DEval(x, chebPolys[i]));
-    }
-    tmp-= (myPoly->coeff[0]/2.0);
-
-
-    return(tmp);
-    */
-}
-
 float psPolynomial1DEval(float x, const psPolynomial1D* myPoly)
 {
@@ -489 +1159 @@
 
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
-        return(p_psOrdPolynomial1DEval(x, myPoly));
+        return(ordPolynomial1DEval(x, myPoly));
     } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
-        return(p_psChebPolynomial1DEval(x, myPoly));
+        return(chebPolynomial1DEval(x, myPoly));
     } else {
-        psError(__func__, "Unknown polynomial type 0x%x\n", myPoly->type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+                PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE,
+                myPoly->type);
     }
     return(0.0);
@@ -520 +1192 @@
 }
 
-
-float p_psOrdPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
+float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
 {
     PS_POLY_CHECK_NULL(myPoly, NAN);
 
-    psS32 loop_x = 0;
-    psS32 loop_y = 0;
-    float polySum = 0.0;
-    float xSum = 1.0;
-    float ySum = 1.0;
-
-    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
-        ySum = xSum;
-        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
-            if (myPoly->mask[loop_x][loop_y] == 0) {
-                polySum += ySum * myPoly->coeff[loop_x][loop_y];
-                ySum *= y;
-            }
-        }
-        xSum *= x;
-    }
-
-    return(polySum);
-}
-
-float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
-{
-    PS_POLY_CHECK_NULL(myPoly, NAN);
-
-    psS32 loop_x = 0;
-    psS32 loop_y = 0;
-    psS32 i = 0;
-    float polySum = 0.0;
-    psPolynomial1D* *chebPolys = NULL;
-    psS32 maxChebyPoly = 0;
-
-    // Determine how many Chebyshev polynomials
-    // are needed, then create them.
-    maxChebyPoly = myPoly->nX;
-    if (myPoly->nY > maxChebyPoly) {
-        maxChebyPoly = myPoly->nY;
-    }
-    chebPolys = CreateChebyshevPolys(maxChebyPoly);
-
-    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
-        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
-            if (myPoly->mask[loop_x][loop_y] == 0) {
-                polySum += myPoly->coeff[loop_x][loop_y] *
-                           psPolynomial1DEval(x, chebPolys[loop_x]) *
-                           psPolynomial1DEval(y, chebPolys[loop_y]);
-            }
-        }
-    }
-    for (i=0;i<maxChebyPoly;i++) {
-        psFree(chebPolys[i]);
-    }
-    psFree(chebPolys);
-    return(polySum);
-}
-
-float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
-{
-    PS_POLY_CHECK_NULL(myPoly, NAN);
-
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
-        return(p_psOrdPolynomial2DEval(x, y, myPoly));
+        return(ordPolynomial2DEval(x, y, myPoly));
     } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
-        return(p_psChebPolynomial2DEval(x, y, myPoly));
+        return(chebPolynomial2DEval(x, y, myPoly));
     } else {
-        psError(__func__, "Unknown polynomial type 0x%x\n", myPoly->type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+                PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE,
+                myPoly->type);
     }
     return(0.0);
 }
-
 
 psVector *psPolynomial2DEvalVector(const psVector *x,
@@ -632 +1245 @@
 }
 
-
-
-float p_psOrdPolynomial3DEval(float x, float y, float 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;
-
-    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
-        ySum = xSum;
-        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
-            zSum = ySum;
-            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
-                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
-                    polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z];
-                    zSum *= z;
-                }
-            }
-            ySum *= y;
-        }
-        xSum *= x;
-    }
-
-    return(polySum);
-}
-
-float p_psChebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
-{
-    psS32 loop_x = 0;
-    psS32 loop_y = 0;
-    psS32 loop_z = 0;
-    psS32 i = 0;
-    float polySum = 0.0;
-    psPolynomial1D* *chebPolys = NULL;
-    psS32 maxChebyPoly = 0;
-
-    // Determine how many Chebyshev polynomials
-    // are needed, then create them.
-    maxChebyPoly = myPoly->nX;
-    if (myPoly->nY > maxChebyPoly) {
-        maxChebyPoly = myPoly->nY;
-    }
-    if (myPoly->nZ > maxChebyPoly) {
-        maxChebyPoly = myPoly->nZ;
-    }
-    chebPolys = CreateChebyshevPolys(maxChebyPoly);
-
-    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
-        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
-            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
-                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
-                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
-                               psPolynomial1DEval(x, chebPolys[loop_x]) *
-                               psPolynomial1DEval(y, chebPolys[loop_y]) *
-                               psPolynomial1DEval(z, chebPolys[loop_z]);
-                }
-            }
-        }
-    }
-
-    for (i=0;i<maxChebyPoly;i++) {
-        psFree(chebPolys[i]);
-    }
-    psFree(chebPolys);
-    return(polySum);
-}
-
 float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
 {
@@ -708 +1250 @@
 
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
-        return(p_psOrdPolynomial3DEval(x, y, z, myPoly));
+        return(ordPolynomial3DEval(x, y, z, myPoly));
     } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
-        return(p_psChebPolynomial3DEval(x, y, z, myPoly));
+        return(chebPolynomial3DEval(x, y, z, myPoly));
     } else {
-        psError(__func__, "Unknown polynomial type 0x%x\n", myPoly->type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+                PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE,
+                myPoly->type);
     }
     return(0.0);
@@ -765 +1309 @@
 }
 
-
-
-
-
-
-float p_psOrdPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
-{
-    psS32 loop_w = 0;
-    psS32 loop_x = 0;
-    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;
-
-    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
-        xSum = wSum;
-        for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
-            ySum = xSum;
-            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
-                zSum = ySum;
-                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
-                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
-                        polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
-                        zSum *= z;
-                    }
-                }
-                ySum *= y;
-            }
-            xSum *= x;
-        }
-        wSum *= w;
-    }
-
-    return(polySum);
-}
-
-float p_psChebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
-{
-    psS32 loop_w = 0;
-    psS32 loop_x = 0;
-    psS32 loop_y = 0;
-    psS32 loop_z = 0;
-    psS32 i = 0;
-    float polySum = 0.0;
-    psPolynomial1D* *chebPolys = NULL;
-    psS32 maxChebyPoly = 0;
-
-    // Determine how many Chebyshev polynomials
-    // are needed, then create them.
-    maxChebyPoly = myPoly->nW;
-    if (myPoly->nX > maxChebyPoly) {
-        maxChebyPoly = myPoly->nX;
-    }
-    if (myPoly->nY > maxChebyPoly) {
-        maxChebyPoly = myPoly->nY;
-    }
-    if (myPoly->nZ > maxChebyPoly) {
-        maxChebyPoly = myPoly->nZ;
-    }
-    chebPolys = CreateChebyshevPolys(maxChebyPoly);
-
-    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
-        for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
-            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
-                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
-                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
-                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
-                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
-                                   psPolynomial1DEval(x, chebPolys[loop_x]) *
-                                   psPolynomial1DEval(y, chebPolys[loop_y]) *
-                                   psPolynomial1DEval(z, chebPolys[loop_z]);
-                    }
-                }
-            }
-        }
-    }
-
-    for (i=0;i<maxChebyPoly;i++) {
-        psFree(chebPolys[i]);
-    }
-    psFree(chebPolys);
-    return(polySum);
-}
-
 float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
 {
@@ -857 +1314 @@
 
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
-        return(p_psOrdPolynomial4DEval(w,x,y,z, myPoly));
+        return(ordPolynomial4DEval(w,x,y,z, myPoly));
     } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
-        return(p_psChebPolynomial4DEval(w,x,y,z, myPoly));
+        return(chebPolynomial4DEval(w,x,y,z, myPoly));
     } else {
-        psError(__func__, "Unknown polynomial type 0x%x\n", myPoly->type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+                PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE,
+                myPoly->type);
     }
     return(0.0);
@@ -924 +1383 @@
     return(tmp);
 }
-
-
-
 
 
@@ -1092 +1548 @@
 }
 
-static void dPolynomial1DFree(psDPolynomial1D* myPoly)
-{
-    psFree(myPoly->coeff);
-    psFree(myPoly->coeffErr);
-    psFree(myPoly->mask);
-}
-
-static void dPolynomial2DFree(psDPolynomial2D* myPoly)
-{
-    psS32 x = 0;
-
-    for (x = 0; x < myPoly->nX; x++) {
-        psFree(myPoly->coeff[x]);
-        psFree(myPoly->coeffErr[x]);
-        psFree(myPoly->mask[x]);
-    }
-    psFree(myPoly->coeff);
-    psFree(myPoly->coeffErr);
-    psFree(myPoly->mask);
-}
-
-static void dPolynomial3DFree(psDPolynomial3D* myPoly)
-{
-    psS32 x = 0;
-    psS32 y = 0;
-
-    for (x = 0; x < myPoly->nX; x++) {
-        for (y = 0; y < myPoly->nY; y++) {
-            psFree(myPoly->coeff[x][y]);
-            psFree(myPoly->coeffErr[x][y]);
-            psFree(myPoly->mask[x][y]);
-        }
-        psFree(myPoly->coeff[x]);
-        psFree(myPoly->coeffErr[x]);
-        psFree(myPoly->mask[x]);
-    }
-
-    psFree(myPoly->coeff);
-    psFree(myPoly->coeffErr);
-    psFree(myPoly->mask);
-}
-
-static void dPolynomial4DFree(psDPolynomial4D* myPoly)
-{
-    psS32 w = 0;
-    psS32 x = 0;
-    psS32 y = 0;
-
-    for (w = 0; w < myPoly->nW; w++) {
-        for (x = 0; x < myPoly->nX; x++) {
-            for (y = 0; y < myPoly->nY; y++) {
-                psFree(myPoly->coeff[w][x][y]);
-                psFree(myPoly->coeffErr[w][x][y]);
-                psFree(myPoly->mask[w][x][y]);
-            }
-            psFree(myPoly->coeff[w][x]);
-            psFree(myPoly->coeffErr[w][x]);
-            psFree(myPoly->mask[w][x]);
-        }
-        psFree(myPoly->coeff[w]);
-        psFree(myPoly->coeffErr[w]);
-        psFree(myPoly->mask[w]);
-    }
-
-    psFree(myPoly->coeff);
-    psFree(myPoly->coeffErr);
-    psFree(myPoly->mask);
-}
-
-/*****************************************************************************
-    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
- *****************************************************************************/
-double p_psDOrdPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
-{
-    psS32 loop_x = 0;
-    double polySum = 0.0;
-    double xSum = 1.0;
-
-    for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
-        if (myPoly->mask[loop_x] == 0) {
-            polySum += xSum * myPoly->coeff[loop_x];
-            xSum *= x;
-        }
-    }
-
-    return(polySum);
-}
-
-// XXX: You can do this without having to psAlloc() vector d.
-// XXX: How does the mask vector effect Crenshaw's formula?
-double p_psDChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
-{
-    psVector *d;
-    psS32 n;
-    psS32 i;
-    double tmp;
-
-    n = myPoly->n;
-    d = psVectorAlloc(n, PS_TYPE_F64);
-    d->data.F64[n-1] = myPoly->coeff[n-1];
-    d->data.F64[n-2] = (2.0 * x * d->data.F64[n-1]) + myPoly->coeff[n-2];
-    for (i=n-3;i>=1;i--) {
-        d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) -
-                         (d->data.F64[i+2]) +
-                         (myPoly->coeff[i]);
-    }
-
-    tmp = (x * d->data.F64[1]) -
-          (d->data.F64[2]) +
-          (0.5 * myPoly->coeff[0]);
-
-    psFree(d);
-    return(tmp);
-}
 
 double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
@@ -1212 +1554 @@
 
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
-        return(p_psDOrdPolynomial1DEval(x, myPoly));
+        return(dOrdPolynomial1DEval(x, myPoly));
     } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
-        return(p_psDChebPolynomial1DEval(x, myPoly));
+        return(dChebPolynomial1DEval(x, myPoly));
     } else {
-        psError(__func__, "Unknown polynomial type 0x%x\n", myPoly->type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+                PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE,
+                myPoly->type);
     }
     return(0.0);
@@ -1244 +1588 @@
 
 
-
-double p_psDOrdPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
-{
-    psS32 loop_x = 0;
-    psS32 loop_y = 0;
-    double polySum = 0.0;
-    double xSum = 1.0;
-    double ySum = 1.0;
-
-    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
-        ySum = xSum;
-        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
-            if (myPoly->mask[loop_x][loop_y] == 0) {
-                polySum += ySum * myPoly->coeff[loop_x][loop_y];
-                ySum *= y;
-            }
-        }
-        xSum *= x;
-    }
-
-    return(polySum);
-}
-
-double p_psDChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
-{
-    psS32 loop_x = 0;
-    psS32 loop_y = 0;
-    psS32 i = 0;
-    double polySum = 0.0;
-    psPolynomial1D* *chebPolys = NULL;
-    psS32 maxChebyPoly = 0;
-
-    // Determine how many Chebyshev polynomials
-    // are needed, then create them.
-    maxChebyPoly = myPoly->nX;
-    if (myPoly->nY > maxChebyPoly) {
-        maxChebyPoly = myPoly->nY;
-    }
-    chebPolys = CreateChebyshevPolys(maxChebyPoly);
-
-    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
-        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
-            if (myPoly->mask[loop_x][loop_y] == 0) {
-                polySum += myPoly->coeff[loop_x][loop_y] *
-                           psPolynomial1DEval(x, chebPolys[loop_x]) *
-                           psPolynomial1DEval(y, chebPolys[loop_y]);
-            }
-        }
-    }
-
-    for (i=0;i<maxChebyPoly;i++) {
-        psFree(chebPolys[i]);
-    }
-    psFree(chebPolys);
-    return(polySum);
-}
-
 double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
 {
@@ -1306 +1593 @@
 
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
-        return(p_psDOrdPolynomial2DEval(x, y, myPoly));
+        return(dOrdPolynomial2DEval(x, y, myPoly));
     } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
-        return(p_psDChebPolynomial2DEval(x, y, myPoly));
+        return(dChebPolynomial2DEval(x, y, myPoly));
     } else {
-        psError(__func__, "Unknown polynomial type 0x%x\n", myPoly->type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+                PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE,
+                myPoly->type);
     }
     return(0.0);
@@ -1353 +1642 @@
 
 
-
-double p_psDOrdPolynomial3DEval(double x, double y, double 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;
-
-    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
-        ySum = xSum;
-        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
-            zSum = ySum;
-            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
-                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
-                    polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z];
-                    zSum *= z;
-                }
-            }
-            ySum *= y;
-        }
-        xSum *= x;
-    }
-
-    return(polySum);
-}
-
-double p_psDChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
-{
-    psS32 loop_x = 0;
-    psS32 loop_y = 0;
-    psS32 loop_z = 0;
-    psS32 i = 0;
-    double polySum = 0.0;
-    psPolynomial1D* *chebPolys = NULL;
-    psS32 maxChebyPoly = 0;
-
-    // Determine how many Chebyshev polynomials
-    // are needed, then create them.
-    maxChebyPoly = myPoly->nX;
-    if (myPoly->nY > maxChebyPoly) {
-        maxChebyPoly = myPoly->nY;
-    }
-    if (myPoly->nZ > maxChebyPoly) {
-        maxChebyPoly = myPoly->nZ;
-    }
-    chebPolys = CreateChebyshevPolys(maxChebyPoly);
-
-    for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
-        for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
-            for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
-                if (myPoly->mask[loop_x][loop_y][loop_z] == 0) {
-                    polySum += myPoly->coeff[loop_x][loop_y][loop_z] *
-                               psPolynomial1DEval(x, chebPolys[loop_x]) *
-                               psPolynomial1DEval(y, chebPolys[loop_y]) *
-                               psPolynomial1DEval(z, chebPolys[loop_z]);
-                }
-            }
-        }
-    }
-
-    for (i=0;i<maxChebyPoly;i++) {
-        psFree(chebPolys[i]);
-    }
-    psFree(chebPolys);
-    return(polySum);
-}
-
 double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
 {
@@ -1428 +1647 @@
 
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
-        return(p_psDOrdPolynomial3DEval(x, y, z, myPoly));
+        return(dOrdPolynomial3DEval(x, y, z, myPoly));
     } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
-        return(p_psDChebPolynomial3DEval(x, y, z, myPoly));
+        return(dChebPolynomial3DEval(x, y, z, myPoly));
     } else {
-        psError(__func__, "Unknown polynomial type 0x%x\n", myPoly->type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+                PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE,
+                myPoly->type);
     }
     return(0.0);
@@ -1485 +1706 @@
 }
 
-
-
-
-
-
-
-
-double p_psDOrdPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
-{
-    psS32 loop_w = 0;
-    psS32 loop_x = 0;
-    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;
-
-    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
-        xSum = wSum;
-        for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
-            ySum = xSum;
-            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
-                zSum = ySum;
-                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
-                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
-                        polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
-                        zSum *= z;
-                    }
-                }
-                ySum *= y;
-            }
-            xSum *= x;
-        }
-        wSum *= w;
-    }
-
-    return(polySum);
-}
-
-double p_psDChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
-{
-    psS32 loop_w = 0;
-    psS32 loop_x = 0;
-    psS32 loop_y = 0;
-    psS32 loop_z = 0;
-    psS32 i = 0;
-    double polySum = 0.0;
-    psPolynomial1D* *chebPolys = NULL;
-    psS32 maxChebyPoly = 0;
-
-    // Determine how many Chebyshev polynomials
-    // are needed, then create them.
-    maxChebyPoly = myPoly->nW;
-    if (myPoly->nX > maxChebyPoly) {
-        maxChebyPoly = myPoly->nX;
-    }
-    if (myPoly->nY > maxChebyPoly) {
-        maxChebyPoly = myPoly->nY;
-    }
-    if (myPoly->nZ > maxChebyPoly) {
-        maxChebyPoly = myPoly->nZ;
-    }
-    chebPolys = CreateChebyshevPolys(maxChebyPoly);
-
-    for (loop_w = 0; loop_w < myPoly->nW; loop_w++) {
-        for (loop_x = 0; loop_x < myPoly->nX; loop_x++) {
-            for (loop_y = 0; loop_y < myPoly->nY; loop_y++) {
-                for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) {
-                    if (myPoly->mask[loop_w][loop_x][loop_y][loop_z] == 0) {
-                        polySum += myPoly->coeff[loop_w][loop_x][loop_y][loop_z] *
-                                   psPolynomial1DEval(w, chebPolys[loop_w]) *
-                                   psPolynomial1DEval(x, chebPolys[loop_x]) *
-                                   psPolynomial1DEval(y, chebPolys[loop_y]) *
-                                   psPolynomial1DEval(z, chebPolys[loop_z]);
-                    }
-                }
-            }
-        }
-    }
-
-    for (i=0;i<maxChebyPoly;i++) {
-        psFree(chebPolys[i]);
-    }
-    psFree(chebPolys);
-    return(polySum);
-}
-
 double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
 {
@@ -1579 +1711 @@
 
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
-        return(p_psDOrdPolynomial4DEval(w,x,y,z, myPoly));
+        return(dOrdPolynomial4DEval(w,x,y,z, myPoly));
     } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
-        return(p_psDChebPolynomial4DEval(w,x,y,z, myPoly));
+        return(dChebPolynomial4DEval(w,x,y,z, myPoly));
     } else {
-        psError(__func__, "Unknown polynomial type 0x%x\n", myPoly->type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+                PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE,
+                myPoly->type);
     }
     return(0.0);
@@ -1699 +1833 @@
     (tmp->domains)[numSplines] = max;
 
+    p_psMemSetDeallocator(tmp,(psFreeFcn)spline1DFree);
     return(tmp);
 }
 
-// XXX: Have Robert put the dealocator in the memory file.
-psS32 p_psSpline1DFree(psSpline1D *tmpSpline)
-{
-    psS32 i;
-
-    if (tmpSpline == NULL) {
-        return(0);
-    }
-
-    if (tmpSpline->spline != NULL) {
-        for (i=0;i<tmpSpline->n;i++) {
-            psFree((tmpSpline->spline)[i]);
-        }
-        psFree(tmpSpline->spline);
-    }
-
-    if (tmpSpline->p_psDeriv2 != NULL) {
-        psFree(tmpSpline->p_psDeriv2);
-    }
-    psFree(tmpSpline->domains);
-    psFree(tmpSpline);
-
-    return(0);
-}
 
 /*****************************************************************************
@@ -1765 +1876 @@
 
 /*****************************************************************************
-p_psVectorBinDisectF32(): This is a private function which takes as input a
+vectorBinDisectF32(): This is a private function which takes as input a
 vector of floating point data as well as a single floating point values.
 The input vector values are assumed to be non-decreasing (v[i-1] <= v[j] for
@@ -1776 +1887 @@
 XXX: name since we don't take psVectors as input.
  *****************************************************************************/
-psS32 p_psVectorBinDisectF32(float *bins,
-                             psS32 numBins,
-                             float x)
+static psS32 vectorBinDisectF32(float *bins,
+                                psS32 numBins,
+                                float x)
 {
     psS32 min;
@@ -1784 +1895 @@
     psS32 mid;
 
-    psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4,
-            "---- Calling p_psVectorBinDisectF32(%f)\n", x);
+    psTrace(".psLib.dataManip.psFunctions.vectorBinDisectF32", 4,
+            "---- Calling vectorBinDisectF32(%f)\n", x);
 
     if (x < bins[0]) {
         psLogMsg(__func__, PS_LOG_WARN,
-                 "p_psVectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).",
+                 "vectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).",
                  x, bins[0], bins[numBins-1]);
         return(-2);
@@ -1796 +1907 @@
     if (x > bins[numBins-1]) {
         psLogMsg(__func__, PS_LOG_WARN,
-                 "p_psVectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).",
+                 "vectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).",
                  x, bins[0], bins[numBins-1]);
         return(-1);
@@ -1806 +1917 @@
 
     while (min != max) {
-        psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4,
+        psTrace(".psLib.dataManip.psFunctions.vectorBinDisectF32", 4,
                 "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
                 min, mid, max, x, bins[mid]);
 
         if (x == bins[mid]) {
-            psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4,
-                    "---- Exiting p_psVectorBinDisectF32(): bin %d\n", mid);
+            psTrace(".psLib.dataManip.psFunctions.vectorBinDisectF32", 4,
+                    "---- Exiting vectorBinDisectF32(): bin %d\n", mid);
             return(mid);
         } else if (x < bins[mid]) {
@@ -1822 +1933 @@
     }
 
-    psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4,
-            "---- Exiting p_psVectorBinDisectF32(): bin %d\n", min);
+    psTrace(".psLib.dataManip.psFunctions.vectorBinDisectF32", 4,
+            "---- Exiting vectorBinDisectF32(): bin %d\n", min);
     return(min);
 }
 
 /*****************************************************************************
-p_psVectorBinDisectS32(): integer version of above.
+vectorBinDisectS32(): integer version of above.
  *****************************************************************************/
-psS32 p_psVectorBinDisectS32(psS32 *bins,
-                             psS32 numBins,
-                             psS32 x)
+static psS32 vectorBinDisectS32(psS32 *bins,
+                                psS32 numBins,
+                                psS32 x)
 {
     psS32 min;
@@ -1838 +1949 @@
     psS32 mid;
 
-    psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4,
-            "---- Calling p_psVectorBinDisectS32(%f)\n", x);
+    psTrace(".psLib.dataManip.psFunctions.vectorBinDisectS32", 4,
+            "---- Calling vectorBinDisectS32(%f)\n", x);
 
     if ((x < bins[0]) ||
             (x > bins[numBins-1])) {
         psLogMsg(__func__, PS_LOG_WARN,
-                 "p_psVectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).",
+                 "vectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).",
                  x, bins[0], bins[numBins-1]);
         return(-1);
@@ -1854 +1965 @@
 
     while (min != max) {
-        psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4,
+        psTrace(".psLib.dataManip.psFunctions.vectorBinDisectS32", 4,
                 "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
                 min, mid, max, x, bins[mid]);
 
         if (x == bins[mid]) {
-            psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4,
-                    "---- Exiting p_psVectorBinDisectS32(): bin %d\n", min);
+            psTrace(".psLib.dataManip.psFunctions.vectorBinDisectS32", 4,
+                    "---- Exiting vectorBinDisectS32(): bin %d\n", min);
             return(min);
         } else if (x < bins[mid]) {
@@ -1870 +1981 @@
     }
 
-    psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4,
-            "---- Exiting p_psVectorBinDisectS32(): bin %d\n", min);
+    psTrace(".psLib.dataManip.psFunctions.vectorBinDisectS32", 4,
+            "---- Exiting vectorBinDisectS32(): bin %d\n", min);
     return(min);
 }
@@ -1884 +1995 @@
 
     if (x->type.type == PS_TYPE_S32) {
-        return(p_psVectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));
+        return(vectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));
     } else if (x->type.type == PS_TYPE_F32) {
-        return(p_psVectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));
+        return(vectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));
     } else {
-        psError(__func__, "Unallowable data type.");
+        char* strType;
+        PS_TYPE_NAME(strType,x->type.type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE,
+                PS_ERRORTEXT_psFunctions_TYPE_NOT_SUPPORTED,
+                strType);
         return(-2);
     }
     return(-1);
-}
-
-/*****************************************************************************
-p_psInterpolate1D(): This routine will take as input n-element floating
-point arrays domain and range, and the x value, assumed to lie with the
-domain vector.  It produces as output the (n-1)-order LaGrange interpolated
-value of x.
- 
-XXX: do we error check for non-distinct domain values?
- *****************************************************************************/
-float p_ps1DFullInterpolateF32(float *domain,
-                               float *range,
-                               psS32 n,
-                               float x)
-{
-    PS_INT_CHECK_NON_NEGATIVE(n, NAN);
-    PS_PTR_CHECK_NULL(domain, NAN);
-    PS_PTR_CHECK_NULL(range, NAN);
-
-    psS32 i;
-    psS32 m;
-    static psVector *p = NULL;
-    p = psVectorRecycle(p, n, PS_TYPE_F32);
-    p_psMemSetPersistent(p, true);
-    p_psMemSetPersistent(p->data.F32, true);
-    /*
-        psVector *p = psVectorAlloc(n, PS_TYPE_F32);
-        float tmp;
-    */
-
-    psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
-            "---- p_ps1DFullInterpolateF32() begin (%d-order at x=%f) (%d data points)----\n", n-1, x, n);
-
-    for (i=0;i<n;i++) {
-        psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
-                "domain/range is (%f %f)\n", domain[i], range[i]);
-    }
-
-    for (i=0;i<n;i++) {
-        p->data.F32[i] = range[i];
-        psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
-                "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
-
-    }
-
-    // From NR, during each iteration of the m loop, we are computing the
-    // p_{i ... i+m} terms.
-    for (m=1;m<n;m++) {
-        for (i=0;i<n-m;i++) {
-            // From NR: we are computing P_{i ... i+m}
-            p->data.F32[i] = (((x-domain[i+m]) * p->data.F32[i]) +
-                              ((domain[i]-x) * p->data.F32[i+1])) /
-                             (domain[i] - domain[i+m]);
-            //printf("((%f-%f * %f) + (%f-%f * %f)) / (%f - %f)\n", x, domain[i+m], p->data.F32[i], domain[i], x, p->data.F32[i+1], domain[i], domain[i+m]);
-            psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
-                    "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
-        }
-    }
-    psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
-            "---- p_ps1DFullInterpolateF32() end ----\n");
-
-    /*
-        tmp = p->data.F32[0];
-        psFree(p);
-        return(tmp);
-    */
-    return(p->data.F32[0]);
-}
-
-
-/*****************************************************************************
-p_ps1DInterpolateF32(): this is the base 1-D flat memory routine to perform
-LaGrange interpolation.
- *****************************************************************************/
-float p_ps1DInterpolateF32(float *domain,
-                           float *range,
-                           psS32 n,
-                           psS32 order,
-                           float x)
-{
-    psS32 binNum;
-    psS32 numIntPoints = order+1;
-    psS32 origin;
-
-    psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4,
-            "---- p_ps1DInterpolateF32() begin ----\n");
-
-    binNum = p_psVectorBinDisectF32(domain, n, x);
-
-    if (0 == numIntPoints%2) {
-        origin = binNum - ((numIntPoints/2) - 1);
-    } else {
-        origin = binNum - (numIntPoints/2);
-        if ((x-domain[binNum]) > (domain[binNum+1]-x)) {
-            // x is closer to binNum+1.
-            origin = 1 + (binNum - (numIntPoints/2));
-        }
-    }
-    if (origin < 0) {
-        origin = 0;
-    }
-    if ((origin + numIntPoints) > n) {
-        origin = n - numIntPoints;
-    }
-
-    psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4,
-            "---- p_ps1DInterpolateF32() end ----\n");
-    return(p_ps1DFullInterpolateF32(&domain[origin], &range[origin], order+1, x));
 }
 
@@ -2035 +2042 @@
 
     if (order > (domain->n - 1)) {
-        psError(__func__, "not enough data points for %d-order interpolation.\n", order);
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
+                PS_ERRORTEXT_psFunctions_NOT_ENOUGH_DATAPOINTS,
+                order);
         return(NULL);
     }
@@ -2042 +2051 @@
         psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
                 "---- p_psVectorInterpolate() end ----\n");
-        return(psScalarAlloc(p_ps1DInterpolateF32(domain->data.F32,
-                             range->data.F32,
-                             domain->n,
-                             order,
-                             x->data.F32), PS_TYPE_F32));
+        return(psScalarAlloc(interpolate1DF32(domain->data.F32,
+                                              range->data.F32,
+                                              domain->n,
+                                              order,
+                                              x->data.F32), PS_TYPE_F32));
     } else if (x->type.type == PS_TYPE_F64) {
         // XXX: use recycled vectors here.
@@ -2053 +2062 @@
 
         psScalar *tmpScalar = psScalarAlloc((double)
-                                            p_ps1DInterpolateF32(domain32->data.F32,
-                                                                 range32->data.F32,
-                                                                 domain32->n,
-                                                                 order,
-                                                                 (float) x->data.F64), PS_TYPE_F64);
+                                            interpolate1DF32(domain32->data.F32,
+                                                             range32->data.F32,
+                                                             domain32->n,
+                                                             order,
+                                                             (float) x->data.F64), PS_TYPE_F64);
         psFree(range32);
         psFree(domain32);
@@ -2067 +2076 @@
 
     } else {
-        // XXX psError: type not supported
-        psError(__func__, "type %d not supported\n", x->type.type);
+        char* strType;
+        PS_TYPE_NAME(strType,x->type.type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE,
+                PS_ERRORTEXT_psFunctions_TYPE_NOT_SUPPORTED,
+                strType);
     }
 
@@ -2084 +2096 @@
 and an independent x value.  Each determines which spline that x corresponds
 to by doing a bracket disection on the domains of the spline data structure
-(p_psVectorBinDisectF32()).  Then it evaluates the spline at that x location
+(vectorBinDisectF32()).  Then it evaluates the spline at that x location
 by a call to the 1D polynomial functions.
  
@@ -2100 +2112 @@
 
     n = spline->n;
-    binNum = p_psVectorBinDisectF32(spline->domains, (spline->n)+1, x);
+    binNum = vectorBinDisectF32(spline->domains, (spline->n)+1, x);
     if (binNum < 0) {
         psLogMsg(__func__, PS_LOG_WARN,
@@ -2139 +2151 @@
         }
     } else {
-        psError(__func__, "Unknown data type.\n");
+        char* strType;
+        PS_TYPE_NAME(strType,x->type.type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE,
+                PS_ERRORTEXT_psFunctions_TYPE_NOT_SUPPORTED,
+                strType);
         return(NULL);
     }