Index: /trunk/psLib/src/dataManip/psConstants.h
===================================================================
--- /trunk/psLib/src/dataManip/psConstants.h	(revision 2216)
+++ /trunk/psLib/src/dataManip/psConstants.h	(revision 2217)
@@ -6,6 +6,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-27 20:32:49 $
+ *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-27 21:06:30 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -28,4 +28,15 @@
 
 /*****************************************************************************
+ 
+*****************************************************************************/
+#define PS_INT_CHECK_NON_NEGATIVE(NAME, RVAL) \
+if (NAME < 0) { \
+    psError(__func__,"%s is less than 0.", #NAME); \
+    return(RVAL); \
+}
+
+
+
+/*****************************************************************************
 Macros which take a generic psLib type and determine if it is NULL, or has
 the wrong type.
@@ -48,7 +59,19 @@
     PS_VECTOR macros:
  *****************************************************************************/
-#define PS_VECTOR_CHECK_TYPE_F32_OR_F64(VEC, RVAL) \
-if ((VEC->type.type != PS_TYPE_F32) && (VEC->type.type != PS_TYPE_F64)) { \
-    psAbort(__func__, "Bad type for VEC (%d)", VEC->type.type); \
+#define PS_VECTOR_CHECK_NULL(NAME, RVAL) \
+if (NAME == NULL || NAME->data.V == NULL) { \
+    psError(__func__,"Unallowable operation: psVector %s or its data is NULL.", #NAME); \
+    return(RVAL); \
+} \
+
+#define PS_VECTOR_CHECK_EMPTY(NAME, RVAL) \
+if (NAME->n < 1) { \
+    psError(__func__,"Unallowable operation: psVector %s has no elements.", #NAME); \
+    return(RVAL); \
+} \
+
+#define PS_VECTOR_CHECK_TYPE_F32_OR_F64(NAME, RVAL) \
+if ((NAME->type.type != PS_TYPE_F32) && (NAME->type.type != PS_TYPE_F64)) { \
+    psAbort(__func__, "psVector %s: bad type(%d)", #NAME, NAME->type.type); \
     return(RVAL); \
 } \
@@ -56,23 +79,11 @@
 #define PS_VECTOR_CHECK_TYPE(NAME, TYPE, RVAL) \
 if (NAME->type.type != TYPE) { \
-    psError(__func__,"Unallowable operation: %s has incorrect type.", #NAME); \
-    return(RVAL); \
-}
-
-#define PS_VECTOR_CHECK_NULL(NAME, RVAL) \
-if (NAME == NULL || NAME->data.V == NULL) { \
-    psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME); \
-    return(RVAL); \
-} \
-
-#define PS_VECTOR_CHECK_EMPTY(NAME, RVAL) \
-if (NAME->n < 1) { \
-    psError(__func__,"Unallowable operation: %s has zero n value.", #NAME); \
-    return(RVAL); \
-} \
+    psError(__func__,"Unallowable operation: psVector %s has incorrect type.", #NAME); \
+    return(RVAL); \
+}
 
 #define PS_VECTOR_CHECK_SIZE_EQUAL(VEC1, VEC2, RVAL) \
 if (VEC1->n != VEC2->n) { \
-    psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
+    psError(__func__,"psVector %s has size %d, psVector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
     return(RVAL); \
 }
@@ -166,8 +177,8 @@
     PS_POLY macros:
 *****************************************************************************/
-#define PS_POLY_CHECK_NULL(NAME, NULL) \
+#define PS_POLY_CHECK_NULL(NAME, RVAL) \
 if (NAME == NULL || NAME->coeff == NULL) { \
     psError(__func__,"Unallowable operation: polynomial %s or its coeffs is NULL.", #NAME); \
-    return(NULL); \
+    return(RVAL); \
 } \
 
Index: /trunk/psLib/src/dataManip/psFunctions.c
===================================================================
--- /trunk/psLib/src/dataManip/psFunctions.c	(revision 2216)
+++ /trunk/psLib/src/dataManip/psFunctions.c	(revision 2217)
@@ -7,8 +7,12 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.54 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-27 00:57:31 $
+ *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-27 21:06:30 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *
+ *  XXX: What happens if the polyEVal functions are called with dat of the wrong
+ *       type?
+ *
  */
 /*****************************************************************************/
@@ -30,53 +34,9 @@
 #include "psAbort.h"
 #include "psLogMsg.h"
-
 #include "psFunctions.h"
+#include "psConstants.h"
 /*****************************************************************************/
 /* DEFINE STATEMENTS                                                         */
 /*****************************************************************************/
-/*
- * The following macros take the vector IN as input.  If IN is not of
- * the correct type, then a new vector OUT is created of the correct
- * type, and the elements of IN and converted to OUT.
- *
- * If vector IN is the correct type, then OUT is simply set to IN.
- */
-#define PS_CONVERT_VECTOR_F32(IN, OUT) \
-if (IN->type.type != PS_TYPE_F32) { \
-    psLogMsg(__func__, PS_LOG_WARN, \
-             "Input vector has incorrect type (%d).\n", \
-             IN->type); \
-    OUT = psVectorAlloc(IN->n, PS_TYPE_F32); \
-    \
-    if (IN->type.type == PS_TYPE_F64) { \
-        for (psS32 i=0;i<IN->n;i++) { \
-            OUT->data.F32[i] = (float) (IN->data.F64[i]); \
-        } \
-    }\
-    else { \
-        psAbort(__func__, "Wrong type.\n"); \
-    } \
-} else { \
-    OUT = (psVector *) IN; \
-}
-
-#define PS_CONVERT_VECTOR_F64(IN, OUT) \
-if (IN->type.type != PS_TYPE_F64) { \
-    psLogMsg(__func__, PS_LOG_WARN, \
-             "Input vector has incorrect type (%d).\n", \
-             IN->type); \
-    OUT = psVectorAlloc(IN->n, PS_TYPE_F64); \
-    \
-    if (IN->type.type == PS_TYPE_F32) { \
-        for (psS32 i=0;i<IN->n;i++) { \
-            OUT->data.F64[i] = (float) (IN->data.F32[i]); \
-        } \
-    }\
-    else { \
-        psAbort(__func__, "Wrong type.\n"); \
-    } \
-} else { \
-    OUT = (psVector *) IN; \
-}
 
 /*****************************************************************************/
@@ -119,4 +79,6 @@
 static psPolynomial1D **CreateChebyshevPolys(psS32 maxChebyPoly)
 {
+    PS_INT_CHECK_NON_NEGATIVE(maxChebyPoly, NULL);
+
     psPolynomial1D **chebPolys = NULL;
     psS32 i = 0;
@@ -187,8 +149,11 @@
 psVector* psGaussianDev(float mean, float sigma, psS32 Npts)
 {
+    PS_INT_CHECK_NON_NEGATIVE(Npts, NULL);
+
     psVector* gauss = NULL;
     const gsl_rng_type *T = NULL;
     gsl_rng *r = NULL;
     psS32 i = 0;
+
 
     gauss = psVectorAlloc(Npts, PS_TYPE_F32);
@@ -212,4 +177,6 @@
                                     psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(n, NULL);
+
     psS32 i = 0;
     psPolynomial1D* newPoly = NULL;
@@ -235,4 +202,7 @@
                                     psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+
     psS32 x = 0;
     psS32 y = 0;
@@ -268,4 +238,8 @@
                                     psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
+
     psS32 x = 0;
     psS32 y = 0;
@@ -310,4 +284,9 @@
                                     psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nW, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
+
     psS32 w = 0;
     psS32 x = 0;
@@ -522,4 +501,6 @@
 float psPolynomial1DEval(float x, const psPolynomial1D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psOrdPolynomial1DEval(x, myPoly));
@@ -535,10 +516,12 @@
                                    const psPolynomial1D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+
     psVector *tmp;
     psVector *myX;
     psS32 i;
 
-    PS_CONVERT_VECTOR_F32(x, myX);
-
     tmp = psVectorAlloc(x->n, PS_TYPE_F32);
     for (i=0;i<x->n;i++) {
@@ -555,4 +538,6 @@
 float p_psOrdPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     psS32 loop_x = 0;
     psS32 loop_y = 0;
@@ -577,4 +562,6 @@
 float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     psS32 loop_x = 0;
     psS32 loop_y = 0;
@@ -610,4 +597,6 @@
 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));
@@ -625,4 +614,10 @@
                                    const psPolynomial2D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
+
     psVector *tmp;
     psVector *myX;
@@ -631,6 +626,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F32(x, myX);
-    PS_CONVERT_VECTOR_F32(y, myY);
     vecLen=x->n;
     if (y->n < vecLen) {
@@ -727,4 +720,6 @@
 float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psOrdPolynomial3DEval(x, y, z, myPoly));
@@ -742,4 +737,12 @@
                                    const psPolynomial3D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(z, NULL);
+    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F32, NULL);
+
     psVector *tmp;
     psVector *myX;
@@ -749,7 +752,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F32(x, myX);
-    PS_CONVERT_VECTOR_F32(y, myY);
-    PS_CONVERT_VECTOR_F32(z, myZ);
     vecLen=x->n;
     if (y->n < vecLen) {
@@ -869,4 +869,6 @@
 float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psOrdPolynomial4DEval(w,x,y,z, myPoly));
@@ -885,4 +887,14 @@
                                    const psPolynomial4D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(w, NULL);
+    PS_VECTOR_CHECK_TYPE(w, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(z, NULL);
+    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F32, NULL);
+
     psVector *tmp;
     psVector *myW;
@@ -893,8 +905,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F32(w, myW);
-    PS_CONVERT_VECTOR_F32(x, myX);
-    PS_CONVERT_VECTOR_F32(y, myY);
-    PS_CONVERT_VECTOR_F32(z, myZ);
     vecLen=w->n;
     if (y->n < vecLen) {
@@ -939,4 +947,6 @@
                                       psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(n, NULL);
+
     psS32 i = 0;
     psDPolynomial1D* newPoly = NULL;
@@ -962,4 +972,7 @@
                                       psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+
     psS32 x = 0;
     psS32 y = 0;
@@ -995,4 +1008,8 @@
                                       psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
+
     psS32 x = 0;
     psS32 y = 0;
@@ -1037,4 +1054,9 @@
                                       psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nW, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
+
     psS32 w = 0;
     psS32 x = 0;
@@ -1207,4 +1229,6 @@
 double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psDOrdPolynomial1DEval(x, myPoly));
@@ -1220,9 +1244,11 @@
                                     const psDPolynomial1D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+
     psVector *tmp;
     psVector *myX;
     psS32 i;
-
-    PS_CONVERT_VECTOR_F64(x, myX);
 
     tmp = psVectorAlloc(x->n, PS_TYPE_F64);
@@ -1297,4 +1323,6 @@
 double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psDOrdPolynomial2DEval(x, y, myPoly));
@@ -1311,4 +1339,10 @@
                                     const psDPolynomial2D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
+
     psVector *tmp;
     psVector *myX;
@@ -1317,6 +1351,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F64(x, myX);
-    PS_CONVERT_VECTOR_F64(y, myY);
     vecLen=x->n;
     if (y->n < vecLen) {
@@ -1413,4 +1445,6 @@
 double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psDOrdPolynomial3DEval(x, y, z, myPoly));
@@ -1428,4 +1462,12 @@
                                     const psDPolynomial3D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(z, NULL);
+    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F64, NULL);
+
     psVector *tmp;
     psVector *myX;
@@ -1435,7 +1477,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F64(x, myX);
-    PS_CONVERT_VECTOR_F64(y, myY);
-    PS_CONVERT_VECTOR_F64(z, myZ);
     vecLen=x->n;
     if (y->n < vecLen) {
@@ -1557,4 +1596,6 @@
 double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psDOrdPolynomial4DEval(w,x,y,z, myPoly));
@@ -1573,4 +1614,14 @@
                                     const psDPolynomial4D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(w, NULL);
+    PS_VECTOR_CHECK_TYPE(w, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(z, NULL);
+    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F64, NULL);
+
     psVector *tmp;
     psVector *myW;
@@ -1581,8 +1632,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F64(w, myW);
-    PS_CONVERT_VECTOR_F64(x, myX);
-    PS_CONVERT_VECTOR_F64(y, myY);
-    PS_CONVERT_VECTOR_F64(z, myZ);
     vecLen=w->n;
     if (y->n < vecLen) {
Index: /trunk/psLib/src/dataManip/psStats.c
===================================================================
--- /trunk/psLib/src/dataManip/psStats.c	(revision 2216)
+++ /trunk/psLib/src/dataManip/psStats.c	(revision 2217)
@@ -9,7 +9,7 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.73 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-27 20:11:47 $
- *
+ *  @version $Revision: 1.74 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-27 21:06:30 $
+n *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -39,19 +39,14 @@
 /* DEFINE STATEMENTS                                                         */
 /*****************************************************************************/
-// will use robust statistical methods.
-#define GAUSS_WIDTH 5       // The width of the Gaussian or boxcar smoothing.
-#define CLIPPED_NUM_ITER_LB 1
-#define CLIPPED_NUM_ITER_UB 10
-#define CLIPPED_SIGMA_LB 1.0
-#define CLIPPED_SIGMA_UB 10.0
-#define true 1
-#define false 0
-#define MY_MAX_FLOAT HUGE
-#define MAX_ITERATIONS 10
-
-void p_psVectorRobustStats(const psVector* restrict myVector,
-                           const psVector* restrict maskVector, psU32 maskVal, psStats* stats);
-
-
+#define PS_GAUSS_WIDTH 5       // The width of the Gaussian or boxcar smoothing.
+#define PS_CLIPPED_NUM_ITER_LB 1
+#define PS_CLIPPED_NUM_ITER_UB 10
+#define PS_CLIPPED_SIGMA_LB 1.0
+#define PS_CLIPPED_SIGMA_UB 10.0
+#define PS_MAX_FLOAT HUGE
+#define PS_POLY_MEDIAN_MAX_ITERATIONS 10
+
+#define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
+(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))
 /*****************************************************************************/
 /* TYPE DEFINITIONS                                                          */
@@ -59,4 +54,8 @@
 psVector* p_psConvertToF32(psVector* in);
 
+void p_psVectorRobustStats(const psVector* restrict myVector,
+                           const psVector* restrict maskVector,
+                           psU32 maskVal,
+                           psStats* stats);
 /*****************************************************************************/
 /* GLOBAL VARIABLES                                                          */
@@ -245,5 +244,5 @@
 {
     psS32 i = 0;                  // Loop index variable
-    float max = -MY_MAX_FLOAT;  // The calculated maximum
+    float max = -PS_MAX_FLOAT;  // The calculated maximum
     float rangeMin = 0.0;       // Exclude data below this
     float rangeMax = 0.0;       // Exclude date above this
@@ -306,5 +305,5 @@
 {
     psS32 i = 0;                  // Loop index variable
-    float min = MY_MAX_FLOAT;   // The calculated maximum
+    float min = PS_MAX_FLOAT;   // The calculated maximum
     float rangeMin = 0.0;       // Exclude data below this
     float rangeMax = 0.0;       // Exclude date above this
@@ -493,5 +492,4 @@
 robustHistogram with a Gaussian of width sigma.
  *****************************************************************************/
-#define GAUSS_WIDTH_OTHER 5.0
 psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram,
                                        float sigma)
@@ -521,5 +519,5 @@
         // warnings?
 
-        x.data.F32 = iMid - (GAUSS_WIDTH * sigma);
+        x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma);
         if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
             jMin = p_psVectorBinDisect(robustHistogram->bounds, &x);
@@ -531,5 +529,5 @@
         }
 
-        x.data.F32 = iMid + (GAUSS_WIDTH * sigma);
+        x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma);
         if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
             jMax = p_psVectorBinDisect(robustHistogram->bounds, &x);
@@ -744,9 +742,9 @@
 
     // Endure that stats->clipIter is within the proper range.
-    if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter) && (stats->clipIter <= CLIPPED_NUM_ITER_UB))) {
+    if (!((PS_CLIPPED_NUM_ITER_LB <= stats->clipIter) && (stats->clipIter <= PS_CLIPPED_NUM_ITER_UB))) {
         psAbort(__func__, "Unallowed value for clipIter (%d).\n", stats->clipIter);
     }
     // Endure that stats->clipSigma is within the proper range.
-    if (!((CLIPPED_SIGMA_LB <= stats->clipSigma) && (stats->clipSigma <= CLIPPED_SIGMA_UB))) {
+    if (!((PS_CLIPPED_SIGMA_LB <= stats->clipSigma) && (stats->clipSigma <= PS_CLIPPED_SIGMA_UB))) {
         psAbort(__func__, "Unallowed value for clipSigma (%f).\n", stats->clipSigma);
     }
@@ -996,5 +994,5 @@
     // printf("p_ps1DPolyMedian(%f, %f, %f) \n", rangeLow, rangeHigh, getThisValue);
 
-    while (numIterations < MAX_ITERATIONS) {
+    while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) {
         midpoint = (rangeHigh + rangeLow) / 2.0;
         if (fabs(midpoint - oldMidpoint) <= FLT_EPSILON) {
@@ -1099,7 +1097,4 @@
     return(tmpFloat);
 }
-
-#define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
-(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))
 
 /******************************************************************************
Index: /trunk/psLib/src/math/psConstants.h
===================================================================
--- /trunk/psLib/src/math/psConstants.h	(revision 2216)
+++ /trunk/psLib/src/math/psConstants.h	(revision 2217)
@@ -6,6 +6,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-27 20:32:49 $
+ *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-27 21:06:30 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -28,4 +28,15 @@
 
 /*****************************************************************************
+ 
+*****************************************************************************/
+#define PS_INT_CHECK_NON_NEGATIVE(NAME, RVAL) \
+if (NAME < 0) { \
+    psError(__func__,"%s is less than 0.", #NAME); \
+    return(RVAL); \
+}
+
+
+
+/*****************************************************************************
 Macros which take a generic psLib type and determine if it is NULL, or has
 the wrong type.
@@ -48,7 +59,19 @@
     PS_VECTOR macros:
  *****************************************************************************/
-#define PS_VECTOR_CHECK_TYPE_F32_OR_F64(VEC, RVAL) \
-if ((VEC->type.type != PS_TYPE_F32) && (VEC->type.type != PS_TYPE_F64)) { \
-    psAbort(__func__, "Bad type for VEC (%d)", VEC->type.type); \
+#define PS_VECTOR_CHECK_NULL(NAME, RVAL) \
+if (NAME == NULL || NAME->data.V == NULL) { \
+    psError(__func__,"Unallowable operation: psVector %s or its data is NULL.", #NAME); \
+    return(RVAL); \
+} \
+
+#define PS_VECTOR_CHECK_EMPTY(NAME, RVAL) \
+if (NAME->n < 1) { \
+    psError(__func__,"Unallowable operation: psVector %s has no elements.", #NAME); \
+    return(RVAL); \
+} \
+
+#define PS_VECTOR_CHECK_TYPE_F32_OR_F64(NAME, RVAL) \
+if ((NAME->type.type != PS_TYPE_F32) && (NAME->type.type != PS_TYPE_F64)) { \
+    psAbort(__func__, "psVector %s: bad type(%d)", #NAME, NAME->type.type); \
     return(RVAL); \
 } \
@@ -56,23 +79,11 @@
 #define PS_VECTOR_CHECK_TYPE(NAME, TYPE, RVAL) \
 if (NAME->type.type != TYPE) { \
-    psError(__func__,"Unallowable operation: %s has incorrect type.", #NAME); \
-    return(RVAL); \
-}
-
-#define PS_VECTOR_CHECK_NULL(NAME, RVAL) \
-if (NAME == NULL || NAME->data.V == NULL) { \
-    psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME); \
-    return(RVAL); \
-} \
-
-#define PS_VECTOR_CHECK_EMPTY(NAME, RVAL) \
-if (NAME->n < 1) { \
-    psError(__func__,"Unallowable operation: %s has zero n value.", #NAME); \
-    return(RVAL); \
-} \
+    psError(__func__,"Unallowable operation: psVector %s has incorrect type.", #NAME); \
+    return(RVAL); \
+}
 
 #define PS_VECTOR_CHECK_SIZE_EQUAL(VEC1, VEC2, RVAL) \
 if (VEC1->n != VEC2->n) { \
-    psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
+    psError(__func__,"psVector %s has size %d, psVector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
     return(RVAL); \
 }
@@ -166,8 +177,8 @@
     PS_POLY macros:
 *****************************************************************************/
-#define PS_POLY_CHECK_NULL(NAME, NULL) \
+#define PS_POLY_CHECK_NULL(NAME, RVAL) \
 if (NAME == NULL || NAME->coeff == NULL) { \
     psError(__func__,"Unallowable operation: polynomial %s or its coeffs is NULL.", #NAME); \
-    return(NULL); \
+    return(RVAL); \
 } \
 
Index: /trunk/psLib/src/math/psPolynomial.c
===================================================================
--- /trunk/psLib/src/math/psPolynomial.c	(revision 2216)
+++ /trunk/psLib/src/math/psPolynomial.c	(revision 2217)
@@ -7,8 +7,12 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.54 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-27 00:57:31 $
+ *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-27 21:06:30 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *
+ *  XXX: What happens if the polyEVal functions are called with dat of the wrong
+ *       type?
+ *
  */
 /*****************************************************************************/
@@ -30,53 +34,9 @@
 #include "psAbort.h"
 #include "psLogMsg.h"
-
 #include "psFunctions.h"
+#include "psConstants.h"
 /*****************************************************************************/
 /* DEFINE STATEMENTS                                                         */
 /*****************************************************************************/
-/*
- * The following macros take the vector IN as input.  If IN is not of
- * the correct type, then a new vector OUT is created of the correct
- * type, and the elements of IN and converted to OUT.
- *
- * If vector IN is the correct type, then OUT is simply set to IN.
- */
-#define PS_CONVERT_VECTOR_F32(IN, OUT) \
-if (IN->type.type != PS_TYPE_F32) { \
-    psLogMsg(__func__, PS_LOG_WARN, \
-             "Input vector has incorrect type (%d).\n", \
-             IN->type); \
-    OUT = psVectorAlloc(IN->n, PS_TYPE_F32); \
-    \
-    if (IN->type.type == PS_TYPE_F64) { \
-        for (psS32 i=0;i<IN->n;i++) { \
-            OUT->data.F32[i] = (float) (IN->data.F64[i]); \
-        } \
-    }\
-    else { \
-        psAbort(__func__, "Wrong type.\n"); \
-    } \
-} else { \
-    OUT = (psVector *) IN; \
-}
-
-#define PS_CONVERT_VECTOR_F64(IN, OUT) \
-if (IN->type.type != PS_TYPE_F64) { \
-    psLogMsg(__func__, PS_LOG_WARN, \
-             "Input vector has incorrect type (%d).\n", \
-             IN->type); \
-    OUT = psVectorAlloc(IN->n, PS_TYPE_F64); \
-    \
-    if (IN->type.type == PS_TYPE_F32) { \
-        for (psS32 i=0;i<IN->n;i++) { \
-            OUT->data.F64[i] = (float) (IN->data.F32[i]); \
-        } \
-    }\
-    else { \
-        psAbort(__func__, "Wrong type.\n"); \
-    } \
-} else { \
-    OUT = (psVector *) IN; \
-}
 
 /*****************************************************************************/
@@ -119,4 +79,6 @@
 static psPolynomial1D **CreateChebyshevPolys(psS32 maxChebyPoly)
 {
+    PS_INT_CHECK_NON_NEGATIVE(maxChebyPoly, NULL);
+
     psPolynomial1D **chebPolys = NULL;
     psS32 i = 0;
@@ -187,8 +149,11 @@
 psVector* psGaussianDev(float mean, float sigma, psS32 Npts)
 {
+    PS_INT_CHECK_NON_NEGATIVE(Npts, NULL);
+
     psVector* gauss = NULL;
     const gsl_rng_type *T = NULL;
     gsl_rng *r = NULL;
     psS32 i = 0;
+
 
     gauss = psVectorAlloc(Npts, PS_TYPE_F32);
@@ -212,4 +177,6 @@
                                     psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(n, NULL);
+
     psS32 i = 0;
     psPolynomial1D* newPoly = NULL;
@@ -235,4 +202,7 @@
                                     psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+
     psS32 x = 0;
     psS32 y = 0;
@@ -268,4 +238,8 @@
                                     psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
+
     psS32 x = 0;
     psS32 y = 0;
@@ -310,4 +284,9 @@
                                     psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nW, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
+
     psS32 w = 0;
     psS32 x = 0;
@@ -522,4 +501,6 @@
 float psPolynomial1DEval(float x, const psPolynomial1D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psOrdPolynomial1DEval(x, myPoly));
@@ -535,10 +516,12 @@
                                    const psPolynomial1D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+
     psVector *tmp;
     psVector *myX;
     psS32 i;
 
-    PS_CONVERT_VECTOR_F32(x, myX);
-
     tmp = psVectorAlloc(x->n, PS_TYPE_F32);
     for (i=0;i<x->n;i++) {
@@ -555,4 +538,6 @@
 float p_psOrdPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     psS32 loop_x = 0;
     psS32 loop_y = 0;
@@ -577,4 +562,6 @@
 float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     psS32 loop_x = 0;
     psS32 loop_y = 0;
@@ -610,4 +597,6 @@
 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));
@@ -625,4 +614,10 @@
                                    const psPolynomial2D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
+
     psVector *tmp;
     psVector *myX;
@@ -631,6 +626,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F32(x, myX);
-    PS_CONVERT_VECTOR_F32(y, myY);
     vecLen=x->n;
     if (y->n < vecLen) {
@@ -727,4 +720,6 @@
 float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psOrdPolynomial3DEval(x, y, z, myPoly));
@@ -742,4 +737,12 @@
                                    const psPolynomial3D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(z, NULL);
+    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F32, NULL);
+
     psVector *tmp;
     psVector *myX;
@@ -749,7 +752,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F32(x, myX);
-    PS_CONVERT_VECTOR_F32(y, myY);
-    PS_CONVERT_VECTOR_F32(z, myZ);
     vecLen=x->n;
     if (y->n < vecLen) {
@@ -869,4 +869,6 @@
 float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psOrdPolynomial4DEval(w,x,y,z, myPoly));
@@ -885,4 +887,14 @@
                                    const psPolynomial4D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(w, NULL);
+    PS_VECTOR_CHECK_TYPE(w, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(z, NULL);
+    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F32, NULL);
+
     psVector *tmp;
     psVector *myW;
@@ -893,8 +905,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F32(w, myW);
-    PS_CONVERT_VECTOR_F32(x, myX);
-    PS_CONVERT_VECTOR_F32(y, myY);
-    PS_CONVERT_VECTOR_F32(z, myZ);
     vecLen=w->n;
     if (y->n < vecLen) {
@@ -939,4 +947,6 @@
                                       psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(n, NULL);
+
     psS32 i = 0;
     psDPolynomial1D* newPoly = NULL;
@@ -962,4 +972,7 @@
                                       psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+
     psS32 x = 0;
     psS32 y = 0;
@@ -995,4 +1008,8 @@
                                       psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
+
     psS32 x = 0;
     psS32 y = 0;
@@ -1037,4 +1054,9 @@
                                       psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nW, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
+
     psS32 w = 0;
     psS32 x = 0;
@@ -1207,4 +1229,6 @@
 double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psDOrdPolynomial1DEval(x, myPoly));
@@ -1220,9 +1244,11 @@
                                     const psDPolynomial1D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+
     psVector *tmp;
     psVector *myX;
     psS32 i;
-
-    PS_CONVERT_VECTOR_F64(x, myX);
 
     tmp = psVectorAlloc(x->n, PS_TYPE_F64);
@@ -1297,4 +1323,6 @@
 double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psDOrdPolynomial2DEval(x, y, myPoly));
@@ -1311,4 +1339,10 @@
                                     const psDPolynomial2D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
+
     psVector *tmp;
     psVector *myX;
@@ -1317,6 +1351,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F64(x, myX);
-    PS_CONVERT_VECTOR_F64(y, myY);
     vecLen=x->n;
     if (y->n < vecLen) {
@@ -1413,4 +1445,6 @@
 double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psDOrdPolynomial3DEval(x, y, z, myPoly));
@@ -1428,4 +1462,12 @@
                                     const psDPolynomial3D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(z, NULL);
+    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F64, NULL);
+
     psVector *tmp;
     psVector *myX;
@@ -1435,7 +1477,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F64(x, myX);
-    PS_CONVERT_VECTOR_F64(y, myY);
-    PS_CONVERT_VECTOR_F64(z, myZ);
     vecLen=x->n;
     if (y->n < vecLen) {
@@ -1557,4 +1596,6 @@
 double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psDOrdPolynomial4DEval(w,x,y,z, myPoly));
@@ -1573,4 +1614,14 @@
                                     const psDPolynomial4D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(w, NULL);
+    PS_VECTOR_CHECK_TYPE(w, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(z, NULL);
+    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F64, NULL);
+
     psVector *tmp;
     psVector *myW;
@@ -1581,8 +1632,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F64(w, myW);
-    PS_CONVERT_VECTOR_F64(x, myX);
-    PS_CONVERT_VECTOR_F64(y, myY);
-    PS_CONVERT_VECTOR_F64(z, myZ);
     vecLen=w->n;
     if (y->n < vecLen) {
Index: /trunk/psLib/src/math/psSpline.c
===================================================================
--- /trunk/psLib/src/math/psSpline.c	(revision 2216)
+++ /trunk/psLib/src/math/psSpline.c	(revision 2217)
@@ -7,8 +7,12 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.54 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-27 00:57:31 $
+ *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-27 21:06:30 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *
+ *  XXX: What happens if the polyEVal functions are called with dat of the wrong
+ *       type?
+ *
  */
 /*****************************************************************************/
@@ -30,53 +34,9 @@
 #include "psAbort.h"
 #include "psLogMsg.h"
-
 #include "psFunctions.h"
+#include "psConstants.h"
 /*****************************************************************************/
 /* DEFINE STATEMENTS                                                         */
 /*****************************************************************************/
-/*
- * The following macros take the vector IN as input.  If IN is not of
- * the correct type, then a new vector OUT is created of the correct
- * type, and the elements of IN and converted to OUT.
- *
- * If vector IN is the correct type, then OUT is simply set to IN.
- */
-#define PS_CONVERT_VECTOR_F32(IN, OUT) \
-if (IN->type.type != PS_TYPE_F32) { \
-    psLogMsg(__func__, PS_LOG_WARN, \
-             "Input vector has incorrect type (%d).\n", \
-             IN->type); \
-    OUT = psVectorAlloc(IN->n, PS_TYPE_F32); \
-    \
-    if (IN->type.type == PS_TYPE_F64) { \
-        for (psS32 i=0;i<IN->n;i++) { \
-            OUT->data.F32[i] = (float) (IN->data.F64[i]); \
-        } \
-    }\
-    else { \
-        psAbort(__func__, "Wrong type.\n"); \
-    } \
-} else { \
-    OUT = (psVector *) IN; \
-}
-
-#define PS_CONVERT_VECTOR_F64(IN, OUT) \
-if (IN->type.type != PS_TYPE_F64) { \
-    psLogMsg(__func__, PS_LOG_WARN, \
-             "Input vector has incorrect type (%d).\n", \
-             IN->type); \
-    OUT = psVectorAlloc(IN->n, PS_TYPE_F64); \
-    \
-    if (IN->type.type == PS_TYPE_F32) { \
-        for (psS32 i=0;i<IN->n;i++) { \
-            OUT->data.F64[i] = (float) (IN->data.F32[i]); \
-        } \
-    }\
-    else { \
-        psAbort(__func__, "Wrong type.\n"); \
-    } \
-} else { \
-    OUT = (psVector *) IN; \
-}
 
 /*****************************************************************************/
@@ -119,4 +79,6 @@
 static psPolynomial1D **CreateChebyshevPolys(psS32 maxChebyPoly)
 {
+    PS_INT_CHECK_NON_NEGATIVE(maxChebyPoly, NULL);
+
     psPolynomial1D **chebPolys = NULL;
     psS32 i = 0;
@@ -187,8 +149,11 @@
 psVector* psGaussianDev(float mean, float sigma, psS32 Npts)
 {
+    PS_INT_CHECK_NON_NEGATIVE(Npts, NULL);
+
     psVector* gauss = NULL;
     const gsl_rng_type *T = NULL;
     gsl_rng *r = NULL;
     psS32 i = 0;
+
 
     gauss = psVectorAlloc(Npts, PS_TYPE_F32);
@@ -212,4 +177,6 @@
                                     psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(n, NULL);
+
     psS32 i = 0;
     psPolynomial1D* newPoly = NULL;
@@ -235,4 +202,7 @@
                                     psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+
     psS32 x = 0;
     psS32 y = 0;
@@ -268,4 +238,8 @@
                                     psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
+
     psS32 x = 0;
     psS32 y = 0;
@@ -310,4 +284,9 @@
                                     psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nW, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
+
     psS32 w = 0;
     psS32 x = 0;
@@ -522,4 +501,6 @@
 float psPolynomial1DEval(float x, const psPolynomial1D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psOrdPolynomial1DEval(x, myPoly));
@@ -535,10 +516,12 @@
                                    const psPolynomial1D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+
     psVector *tmp;
     psVector *myX;
     psS32 i;
 
-    PS_CONVERT_VECTOR_F32(x, myX);
-
     tmp = psVectorAlloc(x->n, PS_TYPE_F32);
     for (i=0;i<x->n;i++) {
@@ -555,4 +538,6 @@
 float p_psOrdPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     psS32 loop_x = 0;
     psS32 loop_y = 0;
@@ -577,4 +562,6 @@
 float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     psS32 loop_x = 0;
     psS32 loop_y = 0;
@@ -610,4 +597,6 @@
 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));
@@ -625,4 +614,10 @@
                                    const psPolynomial2D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
+
     psVector *tmp;
     psVector *myX;
@@ -631,6 +626,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F32(x, myX);
-    PS_CONVERT_VECTOR_F32(y, myY);
     vecLen=x->n;
     if (y->n < vecLen) {
@@ -727,4 +720,6 @@
 float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psOrdPolynomial3DEval(x, y, z, myPoly));
@@ -742,4 +737,12 @@
                                    const psPolynomial3D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(z, NULL);
+    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F32, NULL);
+
     psVector *tmp;
     psVector *myX;
@@ -749,7 +752,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F32(x, myX);
-    PS_CONVERT_VECTOR_F32(y, myY);
-    PS_CONVERT_VECTOR_F32(z, myZ);
     vecLen=x->n;
     if (y->n < vecLen) {
@@ -869,4 +869,6 @@
 float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psOrdPolynomial4DEval(w,x,y,z, myPoly));
@@ -885,4 +887,14 @@
                                    const psPolynomial4D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(w, NULL);
+    PS_VECTOR_CHECK_TYPE(w, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NULL);
+    PS_VECTOR_CHECK_NULL(z, NULL);
+    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F32, NULL);
+
     psVector *tmp;
     psVector *myW;
@@ -893,8 +905,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F32(w, myW);
-    PS_CONVERT_VECTOR_F32(x, myX);
-    PS_CONVERT_VECTOR_F32(y, myY);
-    PS_CONVERT_VECTOR_F32(z, myZ);
     vecLen=w->n;
     if (y->n < vecLen) {
@@ -939,4 +947,6 @@
                                       psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(n, NULL);
+
     psS32 i = 0;
     psDPolynomial1D* newPoly = NULL;
@@ -962,4 +972,7 @@
                                       psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+
     psS32 x = 0;
     psS32 y = 0;
@@ -995,4 +1008,8 @@
                                       psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
+
     psS32 x = 0;
     psS32 y = 0;
@@ -1037,4 +1054,9 @@
                                       psPolynomialType type)
 {
+    PS_INT_CHECK_NON_NEGATIVE(nW, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nX, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nY, NULL);
+    PS_INT_CHECK_NON_NEGATIVE(nZ, NULL);
+
     psS32 w = 0;
     psS32 x = 0;
@@ -1207,4 +1229,6 @@
 double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psDOrdPolynomial1DEval(x, myPoly));
@@ -1220,9 +1244,11 @@
                                     const psDPolynomial1D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+
     psVector *tmp;
     psVector *myX;
     psS32 i;
-
-    PS_CONVERT_VECTOR_F64(x, myX);
 
     tmp = psVectorAlloc(x->n, PS_TYPE_F64);
@@ -1297,4 +1323,6 @@
 double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psDOrdPolynomial2DEval(x, y, myPoly));
@@ -1311,4 +1339,10 @@
                                     const psDPolynomial2D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
+
     psVector *tmp;
     psVector *myX;
@@ -1317,6 +1351,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F64(x, myX);
-    PS_CONVERT_VECTOR_F64(y, myY);
     vecLen=x->n;
     if (y->n < vecLen) {
@@ -1413,4 +1445,6 @@
 double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psDOrdPolynomial3DEval(x, y, z, myPoly));
@@ -1428,4 +1462,12 @@
                                     const psDPolynomial3D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(z, NULL);
+    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F64, NULL);
+
     psVector *tmp;
     psVector *myX;
@@ -1435,7 +1477,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F64(x, myX);
-    PS_CONVERT_VECTOR_F64(y, myY);
-    PS_CONVERT_VECTOR_F64(z, myZ);
     vecLen=x->n;
     if (y->n < vecLen) {
@@ -1557,4 +1596,6 @@
 double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NAN);
+
     if (myPoly->type == PS_POLYNOMIAL_ORD) {
         return(p_psDOrdPolynomial4DEval(w,x,y,z, myPoly));
@@ -1573,4 +1614,14 @@
                                     const psDPolynomial4D *myPoly)
 {
+    PS_POLY_CHECK_NULL(myPoly, NULL);
+    PS_VECTOR_CHECK_NULL(w, NULL);
+    PS_VECTOR_CHECK_TYPE(w, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F64, NULL);
+    PS_VECTOR_CHECK_NULL(z, NULL);
+    PS_VECTOR_CHECK_TYPE(z, PS_TYPE_F64, NULL);
+
     psVector *tmp;
     psVector *myW;
@@ -1581,8 +1632,4 @@
     psS32 vecLen=x->n;
 
-    PS_CONVERT_VECTOR_F64(w, myW);
-    PS_CONVERT_VECTOR_F64(x, myX);
-    PS_CONVERT_VECTOR_F64(y, myY);
-    PS_CONVERT_VECTOR_F64(z, myZ);
     vecLen=w->n;
     if (y->n < vecLen) {
Index: /trunk/psLib/src/math/psStats.c
===================================================================
--- /trunk/psLib/src/math/psStats.c	(revision 2216)
+++ /trunk/psLib/src/math/psStats.c	(revision 2217)
@@ -9,7 +9,7 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.73 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-27 20:11:47 $
- *
+ *  @version $Revision: 1.74 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-27 21:06:30 $
+n *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -39,19 +39,14 @@
 /* DEFINE STATEMENTS                                                         */
 /*****************************************************************************/
-// will use robust statistical methods.
-#define GAUSS_WIDTH 5       // The width of the Gaussian or boxcar smoothing.
-#define CLIPPED_NUM_ITER_LB 1
-#define CLIPPED_NUM_ITER_UB 10
-#define CLIPPED_SIGMA_LB 1.0
-#define CLIPPED_SIGMA_UB 10.0
-#define true 1
-#define false 0
-#define MY_MAX_FLOAT HUGE
-#define MAX_ITERATIONS 10
-
-void p_psVectorRobustStats(const psVector* restrict myVector,
-                           const psVector* restrict maskVector, psU32 maskVal, psStats* stats);
-
-
+#define PS_GAUSS_WIDTH 5       // The width of the Gaussian or boxcar smoothing.
+#define PS_CLIPPED_NUM_ITER_LB 1
+#define PS_CLIPPED_NUM_ITER_UB 10
+#define PS_CLIPPED_SIGMA_LB 1.0
+#define PS_CLIPPED_SIGMA_UB 10.0
+#define PS_MAX_FLOAT HUGE
+#define PS_POLY_MEDIAN_MAX_ITERATIONS 10
+
+#define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
+(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))
 /*****************************************************************************/
 /* TYPE DEFINITIONS                                                          */
@@ -59,4 +54,8 @@
 psVector* p_psConvertToF32(psVector* in);
 
+void p_psVectorRobustStats(const psVector* restrict myVector,
+                           const psVector* restrict maskVector,
+                           psU32 maskVal,
+                           psStats* stats);
 /*****************************************************************************/
 /* GLOBAL VARIABLES                                                          */
@@ -245,5 +244,5 @@
 {
     psS32 i = 0;                  // Loop index variable
-    float max = -MY_MAX_FLOAT;  // The calculated maximum
+    float max = -PS_MAX_FLOAT;  // The calculated maximum
     float rangeMin = 0.0;       // Exclude data below this
     float rangeMax = 0.0;       // Exclude date above this
@@ -306,5 +305,5 @@
 {
     psS32 i = 0;                  // Loop index variable
-    float min = MY_MAX_FLOAT;   // The calculated maximum
+    float min = PS_MAX_FLOAT;   // The calculated maximum
     float rangeMin = 0.0;       // Exclude data below this
     float rangeMax = 0.0;       // Exclude date above this
@@ -493,5 +492,4 @@
 robustHistogram with a Gaussian of width sigma.
  *****************************************************************************/
-#define GAUSS_WIDTH_OTHER 5.0
 psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram,
                                        float sigma)
@@ -521,5 +519,5 @@
         // warnings?
 
-        x.data.F32 = iMid - (GAUSS_WIDTH * sigma);
+        x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma);
         if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
             jMin = p_psVectorBinDisect(robustHistogram->bounds, &x);
@@ -531,5 +529,5 @@
         }
 
-        x.data.F32 = iMid + (GAUSS_WIDTH * sigma);
+        x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma);
         if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
             jMax = p_psVectorBinDisect(robustHistogram->bounds, &x);
@@ -744,9 +742,9 @@
 
     // Endure that stats->clipIter is within the proper range.
-    if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter) && (stats->clipIter <= CLIPPED_NUM_ITER_UB))) {
+    if (!((PS_CLIPPED_NUM_ITER_LB <= stats->clipIter) && (stats->clipIter <= PS_CLIPPED_NUM_ITER_UB))) {
         psAbort(__func__, "Unallowed value for clipIter (%d).\n", stats->clipIter);
     }
     // Endure that stats->clipSigma is within the proper range.
-    if (!((CLIPPED_SIGMA_LB <= stats->clipSigma) && (stats->clipSigma <= CLIPPED_SIGMA_UB))) {
+    if (!((PS_CLIPPED_SIGMA_LB <= stats->clipSigma) && (stats->clipSigma <= PS_CLIPPED_SIGMA_UB))) {
         psAbort(__func__, "Unallowed value for clipSigma (%f).\n", stats->clipSigma);
     }
@@ -996,5 +994,5 @@
     // printf("p_ps1DPolyMedian(%f, %f, %f) \n", rangeLow, rangeHigh, getThisValue);
 
-    while (numIterations < MAX_ITERATIONS) {
+    while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) {
         midpoint = (rangeHigh + rangeLow) / 2.0;
         if (fabs(midpoint - oldMidpoint) <= FLT_EPSILON) {
@@ -1099,7 +1097,4 @@
     return(tmpFloat);
 }
-
-#define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
-(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))
 
 /******************************************************************************