Index: trunk/psLib/src/math/psPolynomial.c
===================================================================
--- trunk/psLib/src/math/psPolynomial.c	(revision 737)
+++ trunk/psLib/src/math/psPolynomial.c	(revision 750)
@@ -24,7 +24,24 @@
 psGaussian(float x,
            float mean,
-           float stddev)
-{
-    return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
+           float stddev,
+           int normal)
+{
+    float tmp = 0.0;
+
+    if (nomal == true) {
+        tmp = 1.0 / sqrtf(2.0 * PI * (stddev * stddev));
+        return(tmp * exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
+    } else {
+        return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
+    }
+}
+
+psVector *psGaussianDev(float mean,
+                        float sigma,
+                        int Npts)
+{
+    psVector *gauss = NULL;
+
+    gauss = psVectorAlloc(Npts, PS_TYPE_FLOAT);
 }
 
@@ -204,5 +221,5 @@
 
 /*****************************************************************************
-    Polynomial coefficients will be accessed in (x,y,z) fashion.
+    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
  *****************************************************************************/
 float
@@ -212,7 +229,9 @@
     int loop_x = 0;
     float polySum = 0.0;
+    float xSum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->n;loop_x++) {
-        polySum+= x + myPoly->coeff[loop_x];
+        polySum+= xSum * myPoly->coeff[loop_x];
+        xSum*=x;
     }
 
@@ -229,9 +248,14 @@
     int 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++) {
+        xSum = xSum;
         for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
-            polySum+= x + myPoly->coeff[loop_x][loop_y];
-        }
+            polySum+= ySum * myPoly->coeff[loop_x][loop_y];
+            ySum*=y;
+        }
+        xSum*=x;
     }
 
@@ -249,11 +273,19 @@
     int 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++) {
-                polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z];
-            }
-        }
+                polySum+= zSum * myPoly->coeff[loop_x][loop_y][loop_z];
+                zSum*=z;
+            }
+            ySum*=y;
+        }
+        xSum*=x;
     }
 
@@ -273,13 +305,25 @@
     int 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++) {
-                    polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                    polySum+= zSum *
+                              myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                    zSum*=z;
                 }
-            }
-        }
+                ySum*=y;
+            }
+            xSum*=x;
+        }
+        wSum*=w;
     }
 
@@ -461,13 +505,18 @@
 }
 
+/*****************************************************************************
+    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
+ *****************************************************************************/
 double
-psevalDPolynomial1D(double x,
-                    const psPolynomial1D *myPoly)
+psDEvalPolynomial1D(double x,
+                    const psDPolynomial1D *myPoly)
 {
     int loop_x = 0;
     double polySum = 0.0;
+    double xSum = 1.0;
 
     for (loop_x=0;loop_x<myPoly->n;loop_x++) {
-        polySum+= x + myPoly->coeff[loop_x];
+        polySum+= xSum * myPoly->coeff[loop_x];
+        xSum*=x;
     }
 
@@ -477,16 +526,21 @@
 
 double
-psevalDPolynomial2D(double x,
+psDEvalPolynomial2D(double x,
                     double y,
-                    const psPolynomial2D *myPoly)
+                    const psDPolynomial2D *myPoly)
 {
     int loop_x = 0;
     int 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++) {
+        xSum = xSum;
         for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
-            polySum+= x + myPoly->coeff[loop_x][loop_y];
-        }
+            polySum+= ySum * myPoly->coeff[loop_x][loop_y];
+            ySum*=y;
+        }
+        xSum*=x;
     }
 
@@ -495,8 +549,8 @@
 
 double
-psevalDPolynomial3D(double x,
+psDEvalPolynomial3D(double x,
                     double y,
                     double z,
-                    const psPolynomial3D *myPoly)
+                    const psDPolynomial3D *myPoly)
 {
     int loop_x = 0;
@@ -504,11 +558,19 @@
     int 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++) {
-                polySum+= x + myPoly->coeff[loop_x][loop_y][loop_z];
-            }
-        }
+                polySum+= zSum * myPoly->coeff[loop_x][loop_y][loop_z];
+                zSum*=z;
+            }
+            ySum*=y;
+        }
+        xSum*=x;
     }
 
@@ -517,9 +579,9 @@
 
 double
-psevalDPolynomial4D(double w,
+psDEvalPolynomial4D(double w,
                     double x,
                     double y,
                     double z,
-                    const psPolynomial4D *myPoly)
+                    const psDPolynomial4D *myPoly)
 {
     int loop_w = 0;
@@ -528,15 +590,27 @@
     int 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++) {
-                    polySum+= x + myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                    polySum+= zSum *
+                              myPoly->coeff[loop_w][loop_x][loop_y][loop_z];
+                    zSum*=z;
                 }
-            }
-        }
-    }
-
-    return(polySum);
-}
+                ySum*=y;
+            }
+            xSum*=x;
+        }
+        wSum*=w;
+    }
+
+    return(polySum);
+}