Index: /trunk/psLib/src/dataManip/Makefile
===================================================================
--- /trunk/psLib/src/dataManip/Makefile	(revision 735)
+++ /trunk/psLib/src/dataManip/Makefile	(revision 736)
@@ -6,5 +6,5 @@
 include ../Makefile.Globals
 CFLAGS := $(CFLAGS_RELOC) -I. -I../sysUtils -I../collections -I..
-SRC_OBJS = psStats.o
+SRC_OBJS = psStats.o psFunctions.o
 
 all: $(TARGET_STATIC)
Index: /trunk/psLib/src/dataManip/psFunctions.c
===================================================================
--- /trunk/psLib/src/dataManip/psFunctions.c	(revision 736)
+++ /trunk/psLib/src/dataManip/psFunctions.c	(revision 736)
@@ -0,0 +1,542 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdarg.h>
+#include <float.h>
+#include <math.h>
+
+#include "psMemory.h"
+#include "psVector.h"
+#include "psTrace.h"
+#include "psError.h"
+#include "psAbort.h"
+#include "psFunctions.h"
+#include "psSort.h"
+
+#include "float.h"
+#include <math.h>
+/*****************************************************************************
+    Evaluate a non-normalized Gaussian with the given mean and sigma at the
+    given coordianate.  Note that this is not a Gaussian deviate.  The
+    evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f]
+ *****************************************************************************/
+float
+psGaussian(float x,
+           float mean,
+           float stddev)
+{
+    return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
+}
+
+/*****************************************************************************
+    This routine must allocate memory for the polynomial structures.
+ *****************************************************************************/
+psPolynomial1D *psPolynomial1DAlloc(int n)
+{
+    psPolynomial1D *newPoly = NULL;
+
+    newPoly = (psPolynomial1D *) psAlloc(sizeof(psPolynomial1D));
+    newPoly->n = n;
+    newPoly->coeff    = (float *) psAlloc(n * sizeof(float));
+    newPoly->coeffErr = (float *) psAlloc(n * sizeof(float));
+    newPoly->mask     = (char *)  psAlloc(n * sizeof(char));
+
+    return(newPoly);
+}
+
+psPolynomial2D *psPolynomial2DAlloc(int nX, int nY)
+{
+    int x = 0;
+    psPolynomial2D *newPoly = NULL;
+
+    newPoly = (psPolynomial2D *) psAlloc(sizeof(psPolynomial2D));
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+
+    newPoly->coeff    = (float **) psAlloc(nX * sizeof(float));
+    newPoly->coeffErr = (float **) psAlloc(nX * sizeof(float));
+    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char));
+    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));
+    }
+
+    return(newPoly);
+}
+
+psPolynomial3D *psPolynomial3DAlloc(int nX, int nY, int nZ)
+{
+    int x = 0;
+    int y = 0;
+    psPolynomial3D *newPoly = NULL;
+
+    newPoly = (psPolynomial3D *) psAlloc(sizeof(psPolynomial3D));
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+    newPoly->nZ = nZ;
+
+    newPoly->coeff    = (float ***) psAlloc(nX * sizeof(float));
+    newPoly->coeffErr = (float ***) psAlloc(nX * sizeof(float));
+    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char));
+    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));
+        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));
+        }
+    }
+
+    return(newPoly);
+}
+
+psPolynomial4D *psPolynomial4DAlloc(int nW, int nX, int nY, int nZ)
+{
+    int w = 0;
+    int x = 0;
+    int y = 0;
+    psPolynomial4D *newPoly = NULL;
+
+    newPoly = (psPolynomial4D *) psAlloc(sizeof(psPolynomial4D));
+    newPoly->nW = nW;
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+    newPoly->nZ = nZ;
+
+    newPoly->coeff    = (float ****) psAlloc(nW * sizeof(float));
+    newPoly->coeffErr = (float ****) psAlloc(nW * sizeof(float));
+    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char));
+    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));
+        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));
+            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));
+            }
+        }
+    }
+
+    return(newPoly);
+}
+
+void psPolynomial1DFree(psPolynomial1D *myPoly)
+{
+    psFree(myPoly->coeff);
+    psFree(myPoly->coeffErr);
+    psFree(myPoly->mask);
+    psFree(myPoly);
+}
+
+void psPolynomial2DFree(psPolynomial2D *myPoly)
+{
+    int 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);
+
+    psFree(myPoly);
+}
+
+void psPolynomial3DFree(psPolynomial3D *myPoly)
+{
+    int x = 0;
+    int 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);
+    psFree(myPoly);
+}
+
+void psPolynomial4DFree(psPolynomial4D *myPoly)
+{
+    int w = 0;
+    int x = 0;
+    int 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);
+    psFree(myPoly);
+}
+
+/*****************************************************************************
+    Polynomial coefficients will be accessed in (x,y,z) fashion.
+ *****************************************************************************/
+float
+psEvalPolynomial1D(float x,
+                   const psPolynomial1D *myPoly)
+{
+    int loop_x = 0;
+    float polySum = 0.0;
+
+    for (loop_x=0;loop_x<myPoly->n;loop_x++) {
+        polySum+= x + myPoly->coeff[loop_x];
+    }
+
+    return(polySum);
+}
+
+
+float
+psEvalPolynomial2D(float x,
+                   float y,
+                   const psPolynomial2D *myPoly)
+{
+    int loop_x = 0;
+    int loop_y = 0;
+    float polySum = 0.0;
+
+    for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+            polySum+= x + myPoly->coeff[loop_x][loop_y];
+        }
+    }
+
+    return(polySum);
+}
+
+float
+psEvalPolynomial3D(float x,
+                   float y,
+                   float z,
+                   const psPolynomial3D *myPoly)
+{
+    int loop_x = 0;
+    int loop_y = 0;
+    int loop_z = 0;
+    float polySum = 0.0;
+
+    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++) {
+                polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z];
+            }
+        }
+    }
+
+    return(polySum);
+}
+
+float
+psEvalPolynomial4D(float w,
+                   float x,
+                   float y,
+                   float z,
+                   const psPolynomial4D *myPoly)
+{
+    int loop_w = 0;
+    int loop_x = 0;
+    int loop_y = 0;
+    int loop_z = 0;
+    float polySum = 0.0;
+
+    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++) {
+                    polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                }
+            }
+        }
+    }
+
+    return(polySum);
+}
+
+
+
+
+psDPolynomial1D *psDPolynomial1DAlloc(int n)
+{
+    psDPolynomial1D *newPoly = NULL;
+
+    newPoly = (psDPolynomial1D *) psAlloc(sizeof(psDPolynomial1D));
+    newPoly->n = n;
+    newPoly->coeff    = (double *) psAlloc(n * sizeof(double));
+    newPoly->coeffErr = (double *) psAlloc(n * sizeof(double));
+    newPoly->mask     = (char *)  psAlloc(n * sizeof(char));
+
+    return(newPoly);
+}
+
+psDPolynomial2D *psDPolynomial2DAlloc(int nX, int nY)
+{
+    int x = 0;
+    psDPolynomial2D *newPoly = NULL;
+
+    newPoly = (psDPolynomial2D *) psAlloc(sizeof(psDPolynomial2D));
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+
+    newPoly->coeff    = (double **) psAlloc(nX * sizeof(double));
+    newPoly->coeffErr = (double **) psAlloc(nX * sizeof(double));
+    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char));
+    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));
+    }
+
+    return(newPoly);
+}
+
+psDPolynomial3D *psDPolynomial3DAlloc(int nX, int nY, int nZ)
+{
+    int x = 0;
+    int y = 0;
+    psDPolynomial3D *newPoly = NULL;
+
+    newPoly = (psDPolynomial3D *) psAlloc(sizeof(psDPolynomial3D));
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+    newPoly->nZ = nZ;
+
+    newPoly->coeff    = (double ***) psAlloc(nX * sizeof(double));
+    newPoly->coeffErr = (double ***) psAlloc(nX * sizeof(double));
+    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char));
+    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));
+        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));
+        }
+    }
+
+    return(newPoly);
+}
+
+psDPolynomial4D *psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ)
+{
+    int w = 0;
+    int x = 0;
+    int y = 0;
+    psDPolynomial4D *newPoly = NULL;
+
+    newPoly = (psDPolynomial4D *) psAlloc(sizeof(psDPolynomial4D));
+    newPoly->nW = nW;
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+    newPoly->nZ = nZ;
+
+    newPoly->coeff    = (double ****) psAlloc(nW * sizeof(double));
+    newPoly->coeffErr = (double ****) psAlloc(nW * sizeof(double));
+    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char));
+    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));
+        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));
+            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));
+            }
+        }
+    }
+
+    return(newPoly);
+}
+
+void psDPolynomial1DFree(psDPolynomial1D *myPoly)
+{
+    psFree(myPoly->coeff);
+    psFree(myPoly->coeffErr);
+    psFree(myPoly->mask);
+    psFree(myPoly);
+}
+
+void psDPolynomial2DFree(psDPolynomial2D *myPoly)
+{
+    int 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);
+
+    psFree(myPoly);
+}
+
+void psDPolynomial3DFree(psDPolynomial3D *myPoly)
+{
+    int x = 0;
+    int 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);
+    psFree(myPoly);
+}
+
+void psDPolynomial4DFree(psDPolynomial4D *myPoly)
+{
+    int w = 0;
+    int x = 0;
+    int 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);
+    psFree(myPoly);
+}
+
+double
+psevalDPolynomial1D(double x,
+                    const psPolynomial1D *myPoly)
+{
+    int loop_x = 0;
+    double polySum = 0.0;
+
+    for (loop_x=0;loop_x<myPoly->n;loop_x++) {
+        polySum+= x + myPoly->coeff[loop_x];
+    }
+
+    return(polySum);
+}
+
+
+double
+psevalDPolynomial2D(double x,
+                    double y,
+                    const psPolynomial2D *myPoly)
+{
+    int loop_x = 0;
+    int loop_y = 0;
+    double polySum = 0.0;
+
+    for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+            polySum+= x + myPoly->coeff[loop_x][loop_y];
+        }
+    }
+
+    return(polySum);
+}
+
+double
+psevalDPolynomial3D(double x,
+                    double y,
+                    double z,
+                    const psPolynomial3D *myPoly)
+{
+    int loop_x = 0;
+    int loop_y = 0;
+    int loop_z = 0;
+    double polySum = 0.0;
+
+    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++) {
+                polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z];
+            }
+        }
+    }
+
+    return(polySum);
+}
+
+double
+psevalDPolynomial4D(double w,
+                    double x,
+                    double y,
+                    double z,
+                    const psPolynomial4D *myPoly)
+{
+    int loop_w = 0;
+    int loop_x = 0;
+    int loop_y = 0;
+    int loop_z = 0;
+    double polySum = 0.0;
+
+    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++) {
+                    polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                }
+            }
+        }
+    }
+
+    return(polySum);
+}
Index: /trunk/psLib/src/dataManip/psFunctions.h
===================================================================
--- /trunk/psLib/src/dataManip/psFunctions.h	(revision 736)
+++ /trunk/psLib/src/dataManip/psFunctions.h	(revision 736)
@@ -0,0 +1,235 @@
+#if !defined(PS_FUNCTIONS_H)
+#define PS_FUNCTIONS_H
+
+/** \file psFunctions.h
+ *  \brief Standard Mathematical Functions. 
+ *  \ingroup MathGroup
+ */
+
+/** Evaluate a non-normalized Gaussian with the given mean and sigma at the
+    given coordianate.  Note that this is not a Gaussian deviate.  The
+    evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] */
+
+float
+psGaussian(float x,   ///< Value at which to evaluate
+           float mean,   ///< Mean for the Gaussian
+           float stddev   ///< Standard deviation for the Gaussian
+          );
+
+/*****************************************************************************/
+
+/** One-dimensional polynomial */
+typedef struct
+{
+    int n;    ///< Number of terms
+    float *coeff;   ///< Coefficients
+    float *coeffErr;   ///< Error in coefficients
+    char *mask;    ///< Coefficient mask
+}
+psPolynomial1D;
+
+/** Two-dimensional polynomial */
+typedef struct
+{
+    int nX, nY;    ///< Number of terms in x and y
+    float **coeff;   ///< Coefficients
+    float **coeffErr;   ///< Error in coefficients
+    char **mask;   ///< Coefficients mask
+}
+psPolynomial2D;
+
+/** Three-dimensional polynomial */
+typedef struct
+{
+    int nX, nY, nZ;   ///< Number of terms in x, y and z
+    float ***coeff;    ///< Coefficients
+    float ***coeffErr;    ///< Error in coefficients
+    char ***mask;    ///< Coefficients mask
+}
+psPolynomial3D;
+
+/** Four-dimensional polynomial */
+typedef struct
+{
+    int nW, nX, nY, nZ;   ///< Number of terms in w, x, y and z
+    float ****coeff;    ///< Coefficients
+    float ****coeffErr;   ///< Error in coefficients
+    char ****mask;    ///< Coefficients mask
+}
+psPolynomial4D;
+
+
+/** Functions **************************************************************/
+/** \addtogroup MathGroup Math Utilities
+ *  \{
+ */
+
+/** Constructor */
+psPolynomial1D *psPolynomial1DAlloc(int n ///< Number of terms
+                                   );
+/** Constructor */
+psPolynomial2D *psPolynomial2DAlloc(int nX, int nY ///< Number of terms in x and y
+                                   );
+/** Constructor */
+psPolynomial3D *psPolynomial3DAlloc(int nX, int nY, int nZ ///< Number of terms in x, y and z
+                                   );
+/** Constructor */
+psPolynomial4D *psPolynomial4DAlloc(int nW, int nX, int nY, int nZ ///< Number of terms in w, x, y and z
+                                   );
+
+/** Destructor */
+void psPolynomial1DFree(psPolynomial1D *myPoly ///< Polynomial to destroy
+                       );
+
+/** Destructor */
+void psPolynomial2DFree(psPolynomial2D *myPoly ///< Polynomial to destroy
+                       );
+/** Destructor */
+void psPolynomial3DFree(psPolynomial3D *myPoly ///< Polynomial to destroy
+                       );
+/** Destructor */
+void psPolynomial4DFree(psPolynomial4D *myPoly ///< Polynomial to destroy
+                       );
+
+/** Evaluate 1D polynomial */
+float
+psEvalPolynomial1D(float x,  ///< Value at which to evaluate
+                   const psPolynomial1D *myPoly ///< Coefficients for the polynomial
+                  );
+
+/** Evaluate 2D polynomial */
+float
+psEvalPolynomial2D(float x,  ///< Value x at which to evaluate
+                   float y,  ///< Value y at which to evaluate
+                   const psPolynomial2D *myPoly ///< Coefficients for the polynomial
+                  );
+
+/** Evaluate 3D polynomial */
+float
+psEvalPolynomial3D(float x,  ///< Value x at which to evaluate
+                   float y,  ///< Value y at which to evaluate
+                   float z,  ///< Value z at which to evaluate
+                   const psPolynomial3D *myPoly ///< Coefficients for the polynomial
+                  );
+
+/** Evaluate 4D polynomial */
+float
+psEvalPolynomial4D(float w,  ///< Value w at which to evaluate
+                   float x,  ///< Value x at which to evaluate
+                   float y,  ///< Value y at which to evaluate
+                   float z,  ///< Value z at which to evaluate
+                   const psPolynomial4D *myPoly ///< Coefficients for the polynomial
+                  );
+
+/* \} */ // End of MathGroup Functions
+
+/*****************************************************************************/
+
+/* Double-precision polynomials, mainly for use in astrometry */
+
+/** Double-precision one-dimensional polynomial */
+typedef struct
+{
+    int n;    ///< Number of terms
+    double *coeff;   ///< Coefficients
+    double *coeffErr;   ///< Error in coefficients
+    char *mask;    ///< Coefficient mask
+}
+psDPolynomial1D;
+
+/** Double-precision two-dimensional polynomial */
+typedef struct
+{
+    int nX, nY;    ///< Number of terms in x and y
+    double **coeff;   ///< Coefficients
+    double **coeffErr;    ///< Error in coefficients
+    char **mask;   ///< Coefficients mask
+}
+psDPolynomial2D;
+
+/** Double-precision three-dimensional polynomial */
+typedef struct
+{
+    int nX, nY, nZ;   ///< Number of terms in x, y and z
+    double ***coeff;   ///< Coefficients
+    double ***coeffErr;   ///< Error in coefficients
+    char ***mask;    ///< Coefficient mask
+}
+psDPolynomial3D;
+
+/** Double-precision four-dimensional polynomial */
+typedef struct
+{
+    int nW, nX, nY, nZ;   ///< Number of terms in w, x, y and z
+    double ****coeff;    ///< Coefficients
+    double ****coeffErr;   ///< Error in coefficients
+    char ****mask;    ///< Coefficients mask
+}
+psDPolynomial4D;
+
+/** Functions **************************************************************/
+/** \addtogroup MathGroup Math Utilities
+ *  \{
+ */
+
+/** Constructor */
+psDPolynomial1D *psDPolynomial1DAlloc(int n ///< Number of terms
+                                     );
+/** Constructor */
+psDPolynomial2D *psDPolynomial2DAlloc(int nX, int nY ///< Number of terms in x and y
+                                     );
+/** Constructor */
+psDPolynomial3D *psDPolynomial3DAlloc(int nX, int nY, int nZ ///< Number of terms in x, y and z
+                                     );
+/** Constructor */
+psDPolynomial4D *psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ ///< Number of terms in w, x, y and z
+                                     );
+
+
+/** Destructor */
+void psDPolynomial1DFree(psDPolynomial1D *myPoly ///< Polynomial to destroy
+                        );
+/** Destructor */
+void psDPolynomial2DFree(psDPolynomial2D *myPoly ///< Polynomial to destroy
+                        );
+/** Destructor */
+void psDPolynomial3DFree(psDPolynomial3D *myPoly ///< Polynomial to destroy
+                        );
+/** Destructor */
+void psDPolynomial4DFree(psDPolynomial4D *myPoly ///< Polynomial to destroy
+                        );
+
+
+/** Evaluate 1D polynomial (double precision) */
+double
+psEvalDPolynomial1D(double x,  ///< Value at which to evaluate
+                    const psDPolynomial1D *myPoly ///< Coefficients for the polynomial
+                   );
+
+/** Evaluate 2D polynomial (double precision) */
+double
+psEvalDPolynomial2D(double x,  ///< Value x at which to evaluate
+                    double y,  ///< Value y at which to evaluate
+                    const psDPolynomial2D *myPoly ///< Coefficients for the polynomial
+                   );
+
+/** Evaluate 3D polynomial (double precision) */
+double
+psEvalDPolynomial3D(double x,  ///< Value x at which to evaluate
+                    double y,  ///< Value y at which to evaluate
+                    double z,  ///< Value z at which to evaluate
+                    const psDPolynomial3D *myPoly ///< Coefficients for the polynomial
+                   );
+
+/** Evaluate 4D polynomial (double precision) */
+double
+psEvalDPolynomial4D(double w,  ///< Value w at which to evaluate
+                    double x,  ///< Value x at which to evaluate
+                    double y,  ///< Value y at which to evaluate
+                    double z,  ///< Value z at which to evaluate
+                    const psDPolynomial4D *myPoly ///< Coefficients for the polynomial
+                   );
+
+/* \} */ // End of MathGroup Functions
+
+#endif
Index: /trunk/psLib/src/math/psPolynomial.c
===================================================================
--- /trunk/psLib/src/math/psPolynomial.c	(revision 736)
+++ /trunk/psLib/src/math/psPolynomial.c	(revision 736)
@@ -0,0 +1,542 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdarg.h>
+#include <float.h>
+#include <math.h>
+
+#include "psMemory.h"
+#include "psVector.h"
+#include "psTrace.h"
+#include "psError.h"
+#include "psAbort.h"
+#include "psFunctions.h"
+#include "psSort.h"
+
+#include "float.h"
+#include <math.h>
+/*****************************************************************************
+    Evaluate a non-normalized Gaussian with the given mean and sigma at the
+    given coordianate.  Note that this is not a Gaussian deviate.  The
+    evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f]
+ *****************************************************************************/
+float
+psGaussian(float x,
+           float mean,
+           float stddev)
+{
+    return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
+}
+
+/*****************************************************************************
+    This routine must allocate memory for the polynomial structures.
+ *****************************************************************************/
+psPolynomial1D *psPolynomial1DAlloc(int n)
+{
+    psPolynomial1D *newPoly = NULL;
+
+    newPoly = (psPolynomial1D *) psAlloc(sizeof(psPolynomial1D));
+    newPoly->n = n;
+    newPoly->coeff    = (float *) psAlloc(n * sizeof(float));
+    newPoly->coeffErr = (float *) psAlloc(n * sizeof(float));
+    newPoly->mask     = (char *)  psAlloc(n * sizeof(char));
+
+    return(newPoly);
+}
+
+psPolynomial2D *psPolynomial2DAlloc(int nX, int nY)
+{
+    int x = 0;
+    psPolynomial2D *newPoly = NULL;
+
+    newPoly = (psPolynomial2D *) psAlloc(sizeof(psPolynomial2D));
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+
+    newPoly->coeff    = (float **) psAlloc(nX * sizeof(float));
+    newPoly->coeffErr = (float **) psAlloc(nX * sizeof(float));
+    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char));
+    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));
+    }
+
+    return(newPoly);
+}
+
+psPolynomial3D *psPolynomial3DAlloc(int nX, int nY, int nZ)
+{
+    int x = 0;
+    int y = 0;
+    psPolynomial3D *newPoly = NULL;
+
+    newPoly = (psPolynomial3D *) psAlloc(sizeof(psPolynomial3D));
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+    newPoly->nZ = nZ;
+
+    newPoly->coeff    = (float ***) psAlloc(nX * sizeof(float));
+    newPoly->coeffErr = (float ***) psAlloc(nX * sizeof(float));
+    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char));
+    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));
+        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));
+        }
+    }
+
+    return(newPoly);
+}
+
+psPolynomial4D *psPolynomial4DAlloc(int nW, int nX, int nY, int nZ)
+{
+    int w = 0;
+    int x = 0;
+    int y = 0;
+    psPolynomial4D *newPoly = NULL;
+
+    newPoly = (psPolynomial4D *) psAlloc(sizeof(psPolynomial4D));
+    newPoly->nW = nW;
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+    newPoly->nZ = nZ;
+
+    newPoly->coeff    = (float ****) psAlloc(nW * sizeof(float));
+    newPoly->coeffErr = (float ****) psAlloc(nW * sizeof(float));
+    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char));
+    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));
+        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));
+            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));
+            }
+        }
+    }
+
+    return(newPoly);
+}
+
+void psPolynomial1DFree(psPolynomial1D *myPoly)
+{
+    psFree(myPoly->coeff);
+    psFree(myPoly->coeffErr);
+    psFree(myPoly->mask);
+    psFree(myPoly);
+}
+
+void psPolynomial2DFree(psPolynomial2D *myPoly)
+{
+    int 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);
+
+    psFree(myPoly);
+}
+
+void psPolynomial3DFree(psPolynomial3D *myPoly)
+{
+    int x = 0;
+    int 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);
+    psFree(myPoly);
+}
+
+void psPolynomial4DFree(psPolynomial4D *myPoly)
+{
+    int w = 0;
+    int x = 0;
+    int 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);
+    psFree(myPoly);
+}
+
+/*****************************************************************************
+    Polynomial coefficients will be accessed in (x,y,z) fashion.
+ *****************************************************************************/
+float
+psEvalPolynomial1D(float x,
+                   const psPolynomial1D *myPoly)
+{
+    int loop_x = 0;
+    float polySum = 0.0;
+
+    for (loop_x=0;loop_x<myPoly->n;loop_x++) {
+        polySum+= x + myPoly->coeff[loop_x];
+    }
+
+    return(polySum);
+}
+
+
+float
+psEvalPolynomial2D(float x,
+                   float y,
+                   const psPolynomial2D *myPoly)
+{
+    int loop_x = 0;
+    int loop_y = 0;
+    float polySum = 0.0;
+
+    for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+            polySum+= x + myPoly->coeff[loop_x][loop_y];
+        }
+    }
+
+    return(polySum);
+}
+
+float
+psEvalPolynomial3D(float x,
+                   float y,
+                   float z,
+                   const psPolynomial3D *myPoly)
+{
+    int loop_x = 0;
+    int loop_y = 0;
+    int loop_z = 0;
+    float polySum = 0.0;
+
+    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++) {
+                polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z];
+            }
+        }
+    }
+
+    return(polySum);
+}
+
+float
+psEvalPolynomial4D(float w,
+                   float x,
+                   float y,
+                   float z,
+                   const psPolynomial4D *myPoly)
+{
+    int loop_w = 0;
+    int loop_x = 0;
+    int loop_y = 0;
+    int loop_z = 0;
+    float polySum = 0.0;
+
+    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++) {
+                    polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                }
+            }
+        }
+    }
+
+    return(polySum);
+}
+
+
+
+
+psDPolynomial1D *psDPolynomial1DAlloc(int n)
+{
+    psDPolynomial1D *newPoly = NULL;
+
+    newPoly = (psDPolynomial1D *) psAlloc(sizeof(psDPolynomial1D));
+    newPoly->n = n;
+    newPoly->coeff    = (double *) psAlloc(n * sizeof(double));
+    newPoly->coeffErr = (double *) psAlloc(n * sizeof(double));
+    newPoly->mask     = (char *)  psAlloc(n * sizeof(char));
+
+    return(newPoly);
+}
+
+psDPolynomial2D *psDPolynomial2DAlloc(int nX, int nY)
+{
+    int x = 0;
+    psDPolynomial2D *newPoly = NULL;
+
+    newPoly = (psDPolynomial2D *) psAlloc(sizeof(psDPolynomial2D));
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+
+    newPoly->coeff    = (double **) psAlloc(nX * sizeof(double));
+    newPoly->coeffErr = (double **) psAlloc(nX * sizeof(double));
+    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char));
+    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));
+    }
+
+    return(newPoly);
+}
+
+psDPolynomial3D *psDPolynomial3DAlloc(int nX, int nY, int nZ)
+{
+    int x = 0;
+    int y = 0;
+    psDPolynomial3D *newPoly = NULL;
+
+    newPoly = (psDPolynomial3D *) psAlloc(sizeof(psDPolynomial3D));
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+    newPoly->nZ = nZ;
+
+    newPoly->coeff    = (double ***) psAlloc(nX * sizeof(double));
+    newPoly->coeffErr = (double ***) psAlloc(nX * sizeof(double));
+    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char));
+    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));
+        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));
+        }
+    }
+
+    return(newPoly);
+}
+
+psDPolynomial4D *psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ)
+{
+    int w = 0;
+    int x = 0;
+    int y = 0;
+    psDPolynomial4D *newPoly = NULL;
+
+    newPoly = (psDPolynomial4D *) psAlloc(sizeof(psDPolynomial4D));
+    newPoly->nW = nW;
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+    newPoly->nZ = nZ;
+
+    newPoly->coeff    = (double ****) psAlloc(nW * sizeof(double));
+    newPoly->coeffErr = (double ****) psAlloc(nW * sizeof(double));
+    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char));
+    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));
+        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));
+            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));
+            }
+        }
+    }
+
+    return(newPoly);
+}
+
+void psDPolynomial1DFree(psDPolynomial1D *myPoly)
+{
+    psFree(myPoly->coeff);
+    psFree(myPoly->coeffErr);
+    psFree(myPoly->mask);
+    psFree(myPoly);
+}
+
+void psDPolynomial2DFree(psDPolynomial2D *myPoly)
+{
+    int 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);
+
+    psFree(myPoly);
+}
+
+void psDPolynomial3DFree(psDPolynomial3D *myPoly)
+{
+    int x = 0;
+    int 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);
+    psFree(myPoly);
+}
+
+void psDPolynomial4DFree(psDPolynomial4D *myPoly)
+{
+    int w = 0;
+    int x = 0;
+    int 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);
+    psFree(myPoly);
+}
+
+double
+psevalDPolynomial1D(double x,
+                    const psPolynomial1D *myPoly)
+{
+    int loop_x = 0;
+    double polySum = 0.0;
+
+    for (loop_x=0;loop_x<myPoly->n;loop_x++) {
+        polySum+= x + myPoly->coeff[loop_x];
+    }
+
+    return(polySum);
+}
+
+
+double
+psevalDPolynomial2D(double x,
+                    double y,
+                    const psPolynomial2D *myPoly)
+{
+    int loop_x = 0;
+    int loop_y = 0;
+    double polySum = 0.0;
+
+    for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+            polySum+= x + myPoly->coeff[loop_x][loop_y];
+        }
+    }
+
+    return(polySum);
+}
+
+double
+psevalDPolynomial3D(double x,
+                    double y,
+                    double z,
+                    const psPolynomial3D *myPoly)
+{
+    int loop_x = 0;
+    int loop_y = 0;
+    int loop_z = 0;
+    double polySum = 0.0;
+
+    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++) {
+                polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z];
+            }
+        }
+    }
+
+    return(polySum);
+}
+
+double
+psevalDPolynomial4D(double w,
+                    double x,
+                    double y,
+                    double z,
+                    const psPolynomial4D *myPoly)
+{
+    int loop_w = 0;
+    int loop_x = 0;
+    int loop_y = 0;
+    int loop_z = 0;
+    double polySum = 0.0;
+
+    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++) {
+                    polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                }
+            }
+        }
+    }
+
+    return(polySum);
+}
Index: /trunk/psLib/src/math/psPolynomial.h
===================================================================
--- /trunk/psLib/src/math/psPolynomial.h	(revision 736)
+++ /trunk/psLib/src/math/psPolynomial.h	(revision 736)
@@ -0,0 +1,235 @@
+#if !defined(PS_FUNCTIONS_H)
+#define PS_FUNCTIONS_H
+
+/** \file psFunctions.h
+ *  \brief Standard Mathematical Functions. 
+ *  \ingroup MathGroup
+ */
+
+/** Evaluate a non-normalized Gaussian with the given mean and sigma at the
+    given coordianate.  Note that this is not a Gaussian deviate.  The
+    evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] */
+
+float
+psGaussian(float x,   ///< Value at which to evaluate
+           float mean,   ///< Mean for the Gaussian
+           float stddev   ///< Standard deviation for the Gaussian
+          );
+
+/*****************************************************************************/
+
+/** One-dimensional polynomial */
+typedef struct
+{
+    int n;    ///< Number of terms
+    float *coeff;   ///< Coefficients
+    float *coeffErr;   ///< Error in coefficients
+    char *mask;    ///< Coefficient mask
+}
+psPolynomial1D;
+
+/** Two-dimensional polynomial */
+typedef struct
+{
+    int nX, nY;    ///< Number of terms in x and y
+    float **coeff;   ///< Coefficients
+    float **coeffErr;   ///< Error in coefficients
+    char **mask;   ///< Coefficients mask
+}
+psPolynomial2D;
+
+/** Three-dimensional polynomial */
+typedef struct
+{
+    int nX, nY, nZ;   ///< Number of terms in x, y and z
+    float ***coeff;    ///< Coefficients
+    float ***coeffErr;    ///< Error in coefficients
+    char ***mask;    ///< Coefficients mask
+}
+psPolynomial3D;
+
+/** Four-dimensional polynomial */
+typedef struct
+{
+    int nW, nX, nY, nZ;   ///< Number of terms in w, x, y and z
+    float ****coeff;    ///< Coefficients
+    float ****coeffErr;   ///< Error in coefficients
+    char ****mask;    ///< Coefficients mask
+}
+psPolynomial4D;
+
+
+/** Functions **************************************************************/
+/** \addtogroup MathGroup Math Utilities
+ *  \{
+ */
+
+/** Constructor */
+psPolynomial1D *psPolynomial1DAlloc(int n ///< Number of terms
+                                   );
+/** Constructor */
+psPolynomial2D *psPolynomial2DAlloc(int nX, int nY ///< Number of terms in x and y
+                                   );
+/** Constructor */
+psPolynomial3D *psPolynomial3DAlloc(int nX, int nY, int nZ ///< Number of terms in x, y and z
+                                   );
+/** Constructor */
+psPolynomial4D *psPolynomial4DAlloc(int nW, int nX, int nY, int nZ ///< Number of terms in w, x, y and z
+                                   );
+
+/** Destructor */
+void psPolynomial1DFree(psPolynomial1D *myPoly ///< Polynomial to destroy
+                       );
+
+/** Destructor */
+void psPolynomial2DFree(psPolynomial2D *myPoly ///< Polynomial to destroy
+                       );
+/** Destructor */
+void psPolynomial3DFree(psPolynomial3D *myPoly ///< Polynomial to destroy
+                       );
+/** Destructor */
+void psPolynomial4DFree(psPolynomial4D *myPoly ///< Polynomial to destroy
+                       );
+
+/** Evaluate 1D polynomial */
+float
+psEvalPolynomial1D(float x,  ///< Value at which to evaluate
+                   const psPolynomial1D *myPoly ///< Coefficients for the polynomial
+                  );
+
+/** Evaluate 2D polynomial */
+float
+psEvalPolynomial2D(float x,  ///< Value x at which to evaluate
+                   float y,  ///< Value y at which to evaluate
+                   const psPolynomial2D *myPoly ///< Coefficients for the polynomial
+                  );
+
+/** Evaluate 3D polynomial */
+float
+psEvalPolynomial3D(float x,  ///< Value x at which to evaluate
+                   float y,  ///< Value y at which to evaluate
+                   float z,  ///< Value z at which to evaluate
+                   const psPolynomial3D *myPoly ///< Coefficients for the polynomial
+                  );
+
+/** Evaluate 4D polynomial */
+float
+psEvalPolynomial4D(float w,  ///< Value w at which to evaluate
+                   float x,  ///< Value x at which to evaluate
+                   float y,  ///< Value y at which to evaluate
+                   float z,  ///< Value z at which to evaluate
+                   const psPolynomial4D *myPoly ///< Coefficients for the polynomial
+                  );
+
+/* \} */ // End of MathGroup Functions
+
+/*****************************************************************************/
+
+/* Double-precision polynomials, mainly for use in astrometry */
+
+/** Double-precision one-dimensional polynomial */
+typedef struct
+{
+    int n;    ///< Number of terms
+    double *coeff;   ///< Coefficients
+    double *coeffErr;   ///< Error in coefficients
+    char *mask;    ///< Coefficient mask
+}
+psDPolynomial1D;
+
+/** Double-precision two-dimensional polynomial */
+typedef struct
+{
+    int nX, nY;    ///< Number of terms in x and y
+    double **coeff;   ///< Coefficients
+    double **coeffErr;    ///< Error in coefficients
+    char **mask;   ///< Coefficients mask
+}
+psDPolynomial2D;
+
+/** Double-precision three-dimensional polynomial */
+typedef struct
+{
+    int nX, nY, nZ;   ///< Number of terms in x, y and z
+    double ***coeff;   ///< Coefficients
+    double ***coeffErr;   ///< Error in coefficients
+    char ***mask;    ///< Coefficient mask
+}
+psDPolynomial3D;
+
+/** Double-precision four-dimensional polynomial */
+typedef struct
+{
+    int nW, nX, nY, nZ;   ///< Number of terms in w, x, y and z
+    double ****coeff;    ///< Coefficients
+    double ****coeffErr;   ///< Error in coefficients
+    char ****mask;    ///< Coefficients mask
+}
+psDPolynomial4D;
+
+/** Functions **************************************************************/
+/** \addtogroup MathGroup Math Utilities
+ *  \{
+ */
+
+/** Constructor */
+psDPolynomial1D *psDPolynomial1DAlloc(int n ///< Number of terms
+                                     );
+/** Constructor */
+psDPolynomial2D *psDPolynomial2DAlloc(int nX, int nY ///< Number of terms in x and y
+                                     );
+/** Constructor */
+psDPolynomial3D *psDPolynomial3DAlloc(int nX, int nY, int nZ ///< Number of terms in x, y and z
+                                     );
+/** Constructor */
+psDPolynomial4D *psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ ///< Number of terms in w, x, y and z
+                                     );
+
+
+/** Destructor */
+void psDPolynomial1DFree(psDPolynomial1D *myPoly ///< Polynomial to destroy
+                        );
+/** Destructor */
+void psDPolynomial2DFree(psDPolynomial2D *myPoly ///< Polynomial to destroy
+                        );
+/** Destructor */
+void psDPolynomial3DFree(psDPolynomial3D *myPoly ///< Polynomial to destroy
+                        );
+/** Destructor */
+void psDPolynomial4DFree(psDPolynomial4D *myPoly ///< Polynomial to destroy
+                        );
+
+
+/** Evaluate 1D polynomial (double precision) */
+double
+psEvalDPolynomial1D(double x,  ///< Value at which to evaluate
+                    const psDPolynomial1D *myPoly ///< Coefficients for the polynomial
+                   );
+
+/** Evaluate 2D polynomial (double precision) */
+double
+psEvalDPolynomial2D(double x,  ///< Value x at which to evaluate
+                    double y,  ///< Value y at which to evaluate
+                    const psDPolynomial2D *myPoly ///< Coefficients for the polynomial
+                   );
+
+/** Evaluate 3D polynomial (double precision) */
+double
+psEvalDPolynomial3D(double x,  ///< Value x at which to evaluate
+                    double y,  ///< Value y at which to evaluate
+                    double z,  ///< Value z at which to evaluate
+                    const psDPolynomial3D *myPoly ///< Coefficients for the polynomial
+                   );
+
+/** Evaluate 4D polynomial (double precision) */
+double
+psEvalDPolynomial4D(double w,  ///< Value w at which to evaluate
+                    double x,  ///< Value x at which to evaluate
+                    double y,  ///< Value y at which to evaluate
+                    double z,  ///< Value z at which to evaluate
+                    const psDPolynomial4D *myPoly ///< Coefficients for the polynomial
+                   );
+
+/* \} */ // End of MathGroup Functions
+
+#endif
Index: /trunk/psLib/src/math/psSpline.c
===================================================================
--- /trunk/psLib/src/math/psSpline.c	(revision 736)
+++ /trunk/psLib/src/math/psSpline.c	(revision 736)
@@ -0,0 +1,542 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <stdarg.h>
+#include <float.h>
+#include <math.h>
+
+#include "psMemory.h"
+#include "psVector.h"
+#include "psTrace.h"
+#include "psError.h"
+#include "psAbort.h"
+#include "psFunctions.h"
+#include "psSort.h"
+
+#include "float.h"
+#include <math.h>
+/*****************************************************************************
+    Evaluate a non-normalized Gaussian with the given mean and sigma at the
+    given coordianate.  Note that this is not a Gaussian deviate.  The
+    evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f]
+ *****************************************************************************/
+float
+psGaussian(float x,
+           float mean,
+           float stddev)
+{
+    return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
+}
+
+/*****************************************************************************
+    This routine must allocate memory for the polynomial structures.
+ *****************************************************************************/
+psPolynomial1D *psPolynomial1DAlloc(int n)
+{
+    psPolynomial1D *newPoly = NULL;
+
+    newPoly = (psPolynomial1D *) psAlloc(sizeof(psPolynomial1D));
+    newPoly->n = n;
+    newPoly->coeff    = (float *) psAlloc(n * sizeof(float));
+    newPoly->coeffErr = (float *) psAlloc(n * sizeof(float));
+    newPoly->mask     = (char *)  psAlloc(n * sizeof(char));
+
+    return(newPoly);
+}
+
+psPolynomial2D *psPolynomial2DAlloc(int nX, int nY)
+{
+    int x = 0;
+    psPolynomial2D *newPoly = NULL;
+
+    newPoly = (psPolynomial2D *) psAlloc(sizeof(psPolynomial2D));
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+
+    newPoly->coeff    = (float **) psAlloc(nX * sizeof(float));
+    newPoly->coeffErr = (float **) psAlloc(nX * sizeof(float));
+    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char));
+    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));
+    }
+
+    return(newPoly);
+}
+
+psPolynomial3D *psPolynomial3DAlloc(int nX, int nY, int nZ)
+{
+    int x = 0;
+    int y = 0;
+    psPolynomial3D *newPoly = NULL;
+
+    newPoly = (psPolynomial3D *) psAlloc(sizeof(psPolynomial3D));
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+    newPoly->nZ = nZ;
+
+    newPoly->coeff    = (float ***) psAlloc(nX * sizeof(float));
+    newPoly->coeffErr = (float ***) psAlloc(nX * sizeof(float));
+    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char));
+    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));
+        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));
+        }
+    }
+
+    return(newPoly);
+}
+
+psPolynomial4D *psPolynomial4DAlloc(int nW, int nX, int nY, int nZ)
+{
+    int w = 0;
+    int x = 0;
+    int y = 0;
+    psPolynomial4D *newPoly = NULL;
+
+    newPoly = (psPolynomial4D *) psAlloc(sizeof(psPolynomial4D));
+    newPoly->nW = nW;
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+    newPoly->nZ = nZ;
+
+    newPoly->coeff    = (float ****) psAlloc(nW * sizeof(float));
+    newPoly->coeffErr = (float ****) psAlloc(nW * sizeof(float));
+    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char));
+    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));
+        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));
+            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));
+            }
+        }
+    }
+
+    return(newPoly);
+}
+
+void psPolynomial1DFree(psPolynomial1D *myPoly)
+{
+    psFree(myPoly->coeff);
+    psFree(myPoly->coeffErr);
+    psFree(myPoly->mask);
+    psFree(myPoly);
+}
+
+void psPolynomial2DFree(psPolynomial2D *myPoly)
+{
+    int 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);
+
+    psFree(myPoly);
+}
+
+void psPolynomial3DFree(psPolynomial3D *myPoly)
+{
+    int x = 0;
+    int 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);
+    psFree(myPoly);
+}
+
+void psPolynomial4DFree(psPolynomial4D *myPoly)
+{
+    int w = 0;
+    int x = 0;
+    int 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);
+    psFree(myPoly);
+}
+
+/*****************************************************************************
+    Polynomial coefficients will be accessed in (x,y,z) fashion.
+ *****************************************************************************/
+float
+psEvalPolynomial1D(float x,
+                   const psPolynomial1D *myPoly)
+{
+    int loop_x = 0;
+    float polySum = 0.0;
+
+    for (loop_x=0;loop_x<myPoly->n;loop_x++) {
+        polySum+= x + myPoly->coeff[loop_x];
+    }
+
+    return(polySum);
+}
+
+
+float
+psEvalPolynomial2D(float x,
+                   float y,
+                   const psPolynomial2D *myPoly)
+{
+    int loop_x = 0;
+    int loop_y = 0;
+    float polySum = 0.0;
+
+    for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+            polySum+= x + myPoly->coeff[loop_x][loop_y];
+        }
+    }
+
+    return(polySum);
+}
+
+float
+psEvalPolynomial3D(float x,
+                   float y,
+                   float z,
+                   const psPolynomial3D *myPoly)
+{
+    int loop_x = 0;
+    int loop_y = 0;
+    int loop_z = 0;
+    float polySum = 0.0;
+
+    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++) {
+                polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z];
+            }
+        }
+    }
+
+    return(polySum);
+}
+
+float
+psEvalPolynomial4D(float w,
+                   float x,
+                   float y,
+                   float z,
+                   const psPolynomial4D *myPoly)
+{
+    int loop_w = 0;
+    int loop_x = 0;
+    int loop_y = 0;
+    int loop_z = 0;
+    float polySum = 0.0;
+
+    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++) {
+                    polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                }
+            }
+        }
+    }
+
+    return(polySum);
+}
+
+
+
+
+psDPolynomial1D *psDPolynomial1DAlloc(int n)
+{
+    psDPolynomial1D *newPoly = NULL;
+
+    newPoly = (psDPolynomial1D *) psAlloc(sizeof(psDPolynomial1D));
+    newPoly->n = n;
+    newPoly->coeff    = (double *) psAlloc(n * sizeof(double));
+    newPoly->coeffErr = (double *) psAlloc(n * sizeof(double));
+    newPoly->mask     = (char *)  psAlloc(n * sizeof(char));
+
+    return(newPoly);
+}
+
+psDPolynomial2D *psDPolynomial2DAlloc(int nX, int nY)
+{
+    int x = 0;
+    psDPolynomial2D *newPoly = NULL;
+
+    newPoly = (psDPolynomial2D *) psAlloc(sizeof(psDPolynomial2D));
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+
+    newPoly->coeff    = (double **) psAlloc(nX * sizeof(double));
+    newPoly->coeffErr = (double **) psAlloc(nX * sizeof(double));
+    newPoly->mask     = (char **)  psAlloc(nX * sizeof(char));
+    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));
+    }
+
+    return(newPoly);
+}
+
+psDPolynomial3D *psDPolynomial3DAlloc(int nX, int nY, int nZ)
+{
+    int x = 0;
+    int y = 0;
+    psDPolynomial3D *newPoly = NULL;
+
+    newPoly = (psDPolynomial3D *) psAlloc(sizeof(psDPolynomial3D));
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+    newPoly->nZ = nZ;
+
+    newPoly->coeff    = (double ***) psAlloc(nX * sizeof(double));
+    newPoly->coeffErr = (double ***) psAlloc(nX * sizeof(double));
+    newPoly->mask     = (char ***)  psAlloc(nX * sizeof(char));
+    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));
+        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));
+        }
+    }
+
+    return(newPoly);
+}
+
+psDPolynomial4D *psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ)
+{
+    int w = 0;
+    int x = 0;
+    int y = 0;
+    psDPolynomial4D *newPoly = NULL;
+
+    newPoly = (psDPolynomial4D *) psAlloc(sizeof(psDPolynomial4D));
+    newPoly->nW = nW;
+    newPoly->nX = nX;
+    newPoly->nY = nY;
+    newPoly->nZ = nZ;
+
+    newPoly->coeff    = (double ****) psAlloc(nW * sizeof(double));
+    newPoly->coeffErr = (double ****) psAlloc(nW * sizeof(double));
+    newPoly->mask     = (char ****)  psAlloc(nW * sizeof(char));
+    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));
+        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));
+            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));
+            }
+        }
+    }
+
+    return(newPoly);
+}
+
+void psDPolynomial1DFree(psDPolynomial1D *myPoly)
+{
+    psFree(myPoly->coeff);
+    psFree(myPoly->coeffErr);
+    psFree(myPoly->mask);
+    psFree(myPoly);
+}
+
+void psDPolynomial2DFree(psDPolynomial2D *myPoly)
+{
+    int 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);
+
+    psFree(myPoly);
+}
+
+void psDPolynomial3DFree(psDPolynomial3D *myPoly)
+{
+    int x = 0;
+    int 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);
+    psFree(myPoly);
+}
+
+void psDPolynomial4DFree(psDPolynomial4D *myPoly)
+{
+    int w = 0;
+    int x = 0;
+    int 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);
+    psFree(myPoly);
+}
+
+double
+psevalDPolynomial1D(double x,
+                    const psPolynomial1D *myPoly)
+{
+    int loop_x = 0;
+    double polySum = 0.0;
+
+    for (loop_x=0;loop_x<myPoly->n;loop_x++) {
+        polySum+= x + myPoly->coeff[loop_x];
+    }
+
+    return(polySum);
+}
+
+
+double
+psevalDPolynomial2D(double x,
+                    double y,
+                    const psPolynomial2D *myPoly)
+{
+    int loop_x = 0;
+    int loop_y = 0;
+    double polySum = 0.0;
+
+    for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
+        for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
+            polySum+= x + myPoly->coeff[loop_x][loop_y];
+        }
+    }
+
+    return(polySum);
+}
+
+double
+psevalDPolynomial3D(double x,
+                    double y,
+                    double z,
+                    const psPolynomial3D *myPoly)
+{
+    int loop_x = 0;
+    int loop_y = 0;
+    int loop_z = 0;
+    double polySum = 0.0;
+
+    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++) {
+                polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z];
+            }
+        }
+    }
+
+    return(polySum);
+}
+
+double
+psevalDPolynomial4D(double w,
+                    double x,
+                    double y,
+                    double z,
+                    const psPolynomial4D *myPoly)
+{
+    int loop_w = 0;
+    int loop_x = 0;
+    int loop_y = 0;
+    int loop_z = 0;
+    double polySum = 0.0;
+
+    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++) {
+                    polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                }
+            }
+        }
+    }
+
+    return(polySum);
+}
Index: /trunk/psLib/src/math/psSpline.h
===================================================================
--- /trunk/psLib/src/math/psSpline.h	(revision 736)
+++ /trunk/psLib/src/math/psSpline.h	(revision 736)
@@ -0,0 +1,235 @@
+#if !defined(PS_FUNCTIONS_H)
+#define PS_FUNCTIONS_H
+
+/** \file psFunctions.h
+ *  \brief Standard Mathematical Functions. 
+ *  \ingroup MathGroup
+ */
+
+/** Evaluate a non-normalized Gaussian with the given mean and sigma at the
+    given coordianate.  Note that this is not a Gaussian deviate.  The
+    evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] */
+
+float
+psGaussian(float x,   ///< Value at which to evaluate
+           float mean,   ///< Mean for the Gaussian
+           float stddev   ///< Standard deviation for the Gaussian
+          );
+
+/*****************************************************************************/
+
+/** One-dimensional polynomial */
+typedef struct
+{
+    int n;    ///< Number of terms
+    float *coeff;   ///< Coefficients
+    float *coeffErr;   ///< Error in coefficients
+    char *mask;    ///< Coefficient mask
+}
+psPolynomial1D;
+
+/** Two-dimensional polynomial */
+typedef struct
+{
+    int nX, nY;    ///< Number of terms in x and y
+    float **coeff;   ///< Coefficients
+    float **coeffErr;   ///< Error in coefficients
+    char **mask;   ///< Coefficients mask
+}
+psPolynomial2D;
+
+/** Three-dimensional polynomial */
+typedef struct
+{
+    int nX, nY, nZ;   ///< Number of terms in x, y and z
+    float ***coeff;    ///< Coefficients
+    float ***coeffErr;    ///< Error in coefficients
+    char ***mask;    ///< Coefficients mask
+}
+psPolynomial3D;
+
+/** Four-dimensional polynomial */
+typedef struct
+{
+    int nW, nX, nY, nZ;   ///< Number of terms in w, x, y and z
+    float ****coeff;    ///< Coefficients
+    float ****coeffErr;   ///< Error in coefficients
+    char ****mask;    ///< Coefficients mask
+}
+psPolynomial4D;
+
+
+/** Functions **************************************************************/
+/** \addtogroup MathGroup Math Utilities
+ *  \{
+ */
+
+/** Constructor */
+psPolynomial1D *psPolynomial1DAlloc(int n ///< Number of terms
+                                   );
+/** Constructor */
+psPolynomial2D *psPolynomial2DAlloc(int nX, int nY ///< Number of terms in x and y
+                                   );
+/** Constructor */
+psPolynomial3D *psPolynomial3DAlloc(int nX, int nY, int nZ ///< Number of terms in x, y and z
+                                   );
+/** Constructor */
+psPolynomial4D *psPolynomial4DAlloc(int nW, int nX, int nY, int nZ ///< Number of terms in w, x, y and z
+                                   );
+
+/** Destructor */
+void psPolynomial1DFree(psPolynomial1D *myPoly ///< Polynomial to destroy
+                       );
+
+/** Destructor */
+void psPolynomial2DFree(psPolynomial2D *myPoly ///< Polynomial to destroy
+                       );
+/** Destructor */
+void psPolynomial3DFree(psPolynomial3D *myPoly ///< Polynomial to destroy
+                       );
+/** Destructor */
+void psPolynomial4DFree(psPolynomial4D *myPoly ///< Polynomial to destroy
+                       );
+
+/** Evaluate 1D polynomial */
+float
+psEvalPolynomial1D(float x,  ///< Value at which to evaluate
+                   const psPolynomial1D *myPoly ///< Coefficients for the polynomial
+                  );
+
+/** Evaluate 2D polynomial */
+float
+psEvalPolynomial2D(float x,  ///< Value x at which to evaluate
+                   float y,  ///< Value y at which to evaluate
+                   const psPolynomial2D *myPoly ///< Coefficients for the polynomial
+                  );
+
+/** Evaluate 3D polynomial */
+float
+psEvalPolynomial3D(float x,  ///< Value x at which to evaluate
+                   float y,  ///< Value y at which to evaluate
+                   float z,  ///< Value z at which to evaluate
+                   const psPolynomial3D *myPoly ///< Coefficients for the polynomial
+                  );
+
+/** Evaluate 4D polynomial */
+float
+psEvalPolynomial4D(float w,  ///< Value w at which to evaluate
+                   float x,  ///< Value x at which to evaluate
+                   float y,  ///< Value y at which to evaluate
+                   float z,  ///< Value z at which to evaluate
+                   const psPolynomial4D *myPoly ///< Coefficients for the polynomial
+                  );
+
+/* \} */ // End of MathGroup Functions
+
+/*****************************************************************************/
+
+/* Double-precision polynomials, mainly for use in astrometry */
+
+/** Double-precision one-dimensional polynomial */
+typedef struct
+{
+    int n;    ///< Number of terms
+    double *coeff;   ///< Coefficients
+    double *coeffErr;   ///< Error in coefficients
+    char *mask;    ///< Coefficient mask
+}
+psDPolynomial1D;
+
+/** Double-precision two-dimensional polynomial */
+typedef struct
+{
+    int nX, nY;    ///< Number of terms in x and y
+    double **coeff;   ///< Coefficients
+    double **coeffErr;    ///< Error in coefficients
+    char **mask;   ///< Coefficients mask
+}
+psDPolynomial2D;
+
+/** Double-precision three-dimensional polynomial */
+typedef struct
+{
+    int nX, nY, nZ;   ///< Number of terms in x, y and z
+    double ***coeff;   ///< Coefficients
+    double ***coeffErr;   ///< Error in coefficients
+    char ***mask;    ///< Coefficient mask
+}
+psDPolynomial3D;
+
+/** Double-precision four-dimensional polynomial */
+typedef struct
+{
+    int nW, nX, nY, nZ;   ///< Number of terms in w, x, y and z
+    double ****coeff;    ///< Coefficients
+    double ****coeffErr;   ///< Error in coefficients
+    char ****mask;    ///< Coefficients mask
+}
+psDPolynomial4D;
+
+/** Functions **************************************************************/
+/** \addtogroup MathGroup Math Utilities
+ *  \{
+ */
+
+/** Constructor */
+psDPolynomial1D *psDPolynomial1DAlloc(int n ///< Number of terms
+                                     );
+/** Constructor */
+psDPolynomial2D *psDPolynomial2DAlloc(int nX, int nY ///< Number of terms in x and y
+                                     );
+/** Constructor */
+psDPolynomial3D *psDPolynomial3DAlloc(int nX, int nY, int nZ ///< Number of terms in x, y and z
+                                     );
+/** Constructor */
+psDPolynomial4D *psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ ///< Number of terms in w, x, y and z
+                                     );
+
+
+/** Destructor */
+void psDPolynomial1DFree(psDPolynomial1D *myPoly ///< Polynomial to destroy
+                        );
+/** Destructor */
+void psDPolynomial2DFree(psDPolynomial2D *myPoly ///< Polynomial to destroy
+                        );
+/** Destructor */
+void psDPolynomial3DFree(psDPolynomial3D *myPoly ///< Polynomial to destroy
+                        );
+/** Destructor */
+void psDPolynomial4DFree(psDPolynomial4D *myPoly ///< Polynomial to destroy
+                        );
+
+
+/** Evaluate 1D polynomial (double precision) */
+double
+psEvalDPolynomial1D(double x,  ///< Value at which to evaluate
+                    const psDPolynomial1D *myPoly ///< Coefficients for the polynomial
+                   );
+
+/** Evaluate 2D polynomial (double precision) */
+double
+psEvalDPolynomial2D(double x,  ///< Value x at which to evaluate
+                    double y,  ///< Value y at which to evaluate
+                    const psDPolynomial2D *myPoly ///< Coefficients for the polynomial
+                   );
+
+/** Evaluate 3D polynomial (double precision) */
+double
+psEvalDPolynomial3D(double x,  ///< Value x at which to evaluate
+                    double y,  ///< Value y at which to evaluate
+                    double z,  ///< Value z at which to evaluate
+                    const psDPolynomial3D *myPoly ///< Coefficients for the polynomial
+                   );
+
+/** Evaluate 4D polynomial (double precision) */
+double
+psEvalDPolynomial4D(double w,  ///< Value w at which to evaluate
+                    double x,  ///< Value x at which to evaluate
+                    double y,  ///< Value y at which to evaluate
+                    double z,  ///< Value z at which to evaluate
+                    const psDPolynomial4D *myPoly ///< Coefficients for the polynomial
+                   );
+
+/* \} */ // End of MathGroup Functions
+
+#endif