Index: /trunk/psLib/src/math/Makefile.am
===================================================================
--- /trunk/psLib/src/math/Makefile.am	(revision 6304)
+++ /trunk/psLib/src/math/Makefile.am	(revision 6305)
@@ -15,5 +15,6 @@
 	psPolynomial.c \
 	psSpline.c \
-	psStats.c
+	psStats.c \
+	psMathUtils.c
 
 EXTRA_DIST = math.i
@@ -32,3 +33,4 @@
 	psPolynomial.h \
 	psSpline.h \
-	psStats.h
+	psStats.h \
+	psMathUtils.h
Index: /trunk/psLib/src/math/psMathUtils.c
===================================================================
--- /trunk/psLib/src/math/psMathUtils.c	(revision 6305)
+++ /trunk/psLib/src/math/psMathUtils.c	(revision 6305)
@@ -0,0 +1,405 @@
+/** @file psMathUtils.c
+ *
+ *  This file contains standard math rotines.
+ *
+ *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-02 21:09:07 $
+ *
+ *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
+ *
+ */
+/*****************************************************************************/
+/*  INCLUDE FILES                                                            */
+/*****************************************************************************/
+#include <stdio.h>
+#include <stdbool.h>
+#include <float.h>
+#include <math.h>
+#include "psMemory.h"
+#include "psVector.h"
+#include "psScalar.h"
+#include "psTrace.h"
+#include "psError.h"
+#include "psLogMsg.h"
+#include "psPolynomial.h"
+#include "psMathUtils.h"
+#include "psConstants.h"
+#include "psErrorText.h"
+
+/*****************************************************************************/
+/* DEFINE STATEMENTS                                                         */
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* TYPE DEFINITIONS                                                          */
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* GLOBAL VARIABLES                                                          */
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* FILE STATIC VARIABLES                                                     */
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* FUNCTION IMPLEMENTATION - LOCAL                                           */
+/*****************************************************************************/
+
+/*****************************************************************************
+vectorBinDisectTYPE(): This is a macro for a private function which takes as
+input an array of data as well as a single value for that data.  The input
+vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all i).
+This routine does a binary disection of the vector and returns "i" such that
+(v[i] <= x <= v[i+1).  If x lies outside the range of v[], then this routine
+prints a warning message and returns (-2 or -1).
+ *****************************************************************************/
+#define FUNC_MACRO_VECTOR_BIN_DISECT(TYPE) \
+static psS32 vectorBinDisect##TYPE(ps##TYPE *bins, \
+                                   psS32 numBins, \
+                                   ps##TYPE x) \
+{ \
+    psS32 min; \
+    psS32 max; \
+    psS32 mid; \
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__); \
+    /* psTrace(__func__, 6, "Determining the bin for: %f\n", x); */\
+    if (x < bins[0]) { \
+        psLogMsg(__func__, PS_LOG_WARN, \
+                 "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \
+                 #TYPE, x, bins[0], bins[numBins-1]); \
+        return(-2); \
+    } \
+    if (x > bins[numBins-1]) { \
+        psLogMsg(__func__, PS_LOG_WARN, \
+                 "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \
+                 #TYPE, x, bins[0], bins[numBins-1]); \
+        return(-1); \
+    } \
+    \
+    min = 0; \
+    max = numBins-2; \
+    mid = ((max+1)-min)/2; \
+    while (min != max) { \
+        /*psTrace(__func__, 6, "(min, mid, max) is (%u, %u, %u): (x, bins[mid]) is (%f, %f)\n", min, mid, max, x, bins[mid]); */\
+        \
+        if (x == bins[mid]) { \
+            /*psTrace(__func__, 6, "%.2f should be in this range (%.2f - %.2f)\n", x, bins[mid], bins[mid+1]); */\
+            psTrace(__func__, 4, "---- %s(%d) end (1) ----\n", __func__, mid); \
+            return(mid); \
+        } else if (x < bins[mid]) { \
+            max = mid-1; \
+        } else { \
+            min = mid; \
+        } \
+        mid = ((max+1)+min)/2; \
+    } \
+    /*psTrace(__func__, 6, "%.2f should be in this range (%.2f - %.2f)\n", x, bins[min], bins[min+1]);*/ \
+    psTrace(__func__, 4, "---- %s(%d) end (2) ----\n", __func__, min); \
+    return(min); \
+} \
+
+FUNC_MACRO_VECTOR_BIN_DISECT(S8)
+FUNC_MACRO_VECTOR_BIN_DISECT(S16)
+FUNC_MACRO_VECTOR_BIN_DISECT(S32)
+FUNC_MACRO_VECTOR_BIN_DISECT(S64)
+FUNC_MACRO_VECTOR_BIN_DISECT(U8)
+FUNC_MACRO_VECTOR_BIN_DISECT(U16)
+FUNC_MACRO_VECTOR_BIN_DISECT(U32)
+FUNC_MACRO_VECTOR_BIN_DISECT(U64)
+FUNC_MACRO_VECTOR_BIN_DISECT(F32)
+FUNC_MACRO_VECTOR_BIN_DISECT(F64)
+
+/*****************************************************************************
+p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect().
+  *****************************************************************************/
+psS32 p_psVectorBinDisect(
+    psVector *bins,
+    psScalar *x)
+{
+    PS_ASSERT_VECTOR_NON_NULL(bins, -4);
+    PS_ASSERT_VECTOR_NON_EMPTY(bins, -4);
+    PS_ASSERT_PTR_NON_NULL(x, -6);
+    PS_ASSERT_PTR_TYPE_EQUAL(x, bins, -3);
+    char* strType;
+
+    switch (x->type.type) {
+    case PS_TYPE_U8:
+        return(vectorBinDisectU8(bins->data.U8, bins->n, x->data.U8));
+    case PS_TYPE_U16:
+        return(vectorBinDisectU16(bins->data.U16, bins->n, x->data.U16));
+    case PS_TYPE_U32:
+        return(vectorBinDisectU32(bins->data.U32, bins->n, x->data.U32));
+    case PS_TYPE_U64:
+        return(vectorBinDisectU64(bins->data.U64, bins->n, x->data.U64));
+    case PS_TYPE_S8:
+        return(vectorBinDisectS8(bins->data.S8, bins->n, x->data.S8));
+    case PS_TYPE_S16:
+        return(vectorBinDisectS16(bins->data.S16, bins->n, x->data.S16));
+    case PS_TYPE_S32:
+        return(vectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));
+    case PS_TYPE_S64:
+        return(vectorBinDisectS64(bins->data.S64, bins->n, x->data.S64));
+    case PS_TYPE_F32:
+        return(vectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));
+    case PS_TYPE_F64:
+        return(vectorBinDisectF64(bins->data.F64, bins->n, x->data.F64));
+    default:
+        PS_TYPE_NAME(strType, x->type.type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType);
+        return -1;
+    }
+    return(-3);
+}
+
+
+
+/*****************************************************************************
+fullInterpolateTYPE(): This routine will take as input n-element psVectors
+domain and range, and the x value, assumed to lie with the domain vector.  It
+produces as output the (order-1)-order LaGrange interpolated value of x.
+ 
+XXX: do we error check for non-distinct domain values?
+ 
+This routine will only be called inside this file.  So, we do not have to
+do PS_ASSERTS.
+ *****************************************************************************/
+#define FUNC_MACRO_FULL_INTERPOLATE(TYPE) \
+static psScalar *vectorInterpolate##TYPE( \
+        psScalar *out, \
+        psVector *domain, \
+        psVector *range, \
+        psS32 order, \
+        psScalar *x) \
+{ \
+    psTrace(__func__, 4, "---- %s() begin %u-order.) (%d data points) ----\n", __func__, order, order+1); \
+    if (x->data.TYPE < domain->data.TYPE[0]) { \
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: x is outside the domain of input data.\n"); \
+        out->data.TYPE = range->data.TYPE[0]; \
+        return(out); \
+    } \
+    if (x->data.TYPE > domain->data.TYPE[domain->n-1]) { \
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: x is outside the domain of input data.\n"); \
+        out->data.TYPE = range->data.TYPE[domain->n-1]; \
+        return(out); \
+    } \
+    psVector *p = psVectorCopy(NULL, range, PS_TYPE_##TYPE); \
+    psS32 binNum = p_psVectorBinDisect(domain, x); \
+    \
+    psS32 numIntPoints = order+1; \
+    psS32 origin; \
+    if (0 == numIntPoints%2) { \
+        origin = binNum - ((numIntPoints/2) - 1); \
+    } else { \
+        origin = binNum - (numIntPoints/2); \
+        if ((x->data.TYPE-domain->data.TYPE[binNum]) > (domain->data.TYPE[binNum+1]-x->data.TYPE)) { \
+            /* x is closer to binNum+1. */\
+            origin = 1 + (binNum - (numIntPoints/2)); \
+        } \
+    } \
+    origin = PS_MAX(origin, 0); \
+    origin = PS_MIN(origin, (domain->n - numIntPoints)); \
+    \
+    /* From NR, during each iteration of the m loop, we are computing the p_{i ... i+m} terms. */ \
+    for (psU32 m = 1 ; m < numIntPoints ; m++) { \
+        for (psU32 i = origin ; i < (numIntPoints+origin-m) ; i++) { \
+            p->data.TYPE[i] = (((x->data.TYPE - domain->data.TYPE[i+m]) * p->data.TYPE[i]) + \
+                               ((domain->data.TYPE[i] - x->data.TYPE) * p->data.TYPE[i+1])) / \
+                              (domain->data.TYPE[i] - domain->data.TYPE[i+m]); \
+        } \
+    } \
+    out->data.TYPE = p->data.TYPE[origin]; \
+    psFree(p); \
+    psTrace(__func__, 4, "---- %s(....) end ----\n", __func__); \
+    return(out); \
+} \
+
+FUNC_MACRO_FULL_INTERPOLATE(U8)
+FUNC_MACRO_FULL_INTERPOLATE(U16)
+FUNC_MACRO_FULL_INTERPOLATE(U32)
+FUNC_MACRO_FULL_INTERPOLATE(U64)
+FUNC_MACRO_FULL_INTERPOLATE(S8)
+FUNC_MACRO_FULL_INTERPOLATE(S16)
+FUNC_MACRO_FULL_INTERPOLATE(S32)
+FUNC_MACRO_FULL_INTERPOLATE(S64)
+FUNC_MACRO_FULL_INTERPOLATE(F32)
+FUNC_MACRO_FULL_INTERPOLATE(F64)
+
+/*****************************************************************************
+p_psVectorInterpolate(): This routine will take as input psVectors domain and
+range, and the x value, assumed to lie with the domain vector.  It produces
+as output the LaGrange interpolated value of a polynomial of the specified
+order around the point x.
+ 
+XXX: This stuff does not currently work with a mask.
+ *****************************************************************************/
+psScalar *p_psVectorInterpolate(
+    psScalar *out,
+    psVector *domain,
+    psVector *range,
+    psS32 order,
+    psScalar *x)
+{
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_VECTOR_NON_NULL(domain, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(range, NULL);
+    PS_ASSERT_PTR_NON_NULL(x, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(order, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(domain, range, NULL);
+    PS_ASSERT_PTR_TYPE_EQUAL(domain, range, NULL);
+    PS_ASSERT_PTR_TYPE_EQUAL(domain, x, NULL);
+    if (out == NULL) {
+        out = psScalarAlloc(0, x->type.type);
+    } else {
+        PS_ASSERT_PTR_TYPE_EQUAL(domain, out, NULL);
+    }
+    char* strType;
+
+    switch (x->type.type) {
+    case PS_TYPE_U8:
+        return(vectorInterpolateU8(out, domain, range, order, x));
+    case PS_TYPE_U16:
+        return(vectorInterpolateU16(out, domain, range, order, x));
+    case PS_TYPE_U32:
+        return(vectorInterpolateU32(out, domain, range, order, x));
+    case PS_TYPE_U64:
+        return(vectorInterpolateU64(out, domain, range, order, x));
+    case PS_TYPE_S8:
+        return(vectorInterpolateS8(out, domain, range, order, x));
+    case PS_TYPE_S16:
+        return(vectorInterpolateS16(out, domain, range, order, x));
+    case PS_TYPE_S32:
+        return(vectorInterpolateS32(out, domain, range, order, x));
+    case PS_TYPE_S64:
+        return(vectorInterpolateS64(out, domain, range, order, x));
+    case PS_TYPE_F32:
+        return(vectorInterpolateF32(out, domain, range, order, x));
+    case PS_TYPE_F64:
+        return(vectorInterpolateF64(out, domain, range, order, x));
+    default:
+        PS_TYPE_NAME(strType, x->type.type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType);
+        return NULL;
+    }
+
+    return NULL;
+}
+
+/*****************************************************************************
+These macros and functions define the following functions:
+ 
+    p_psNormalizeVectorRange(myData, low, high)
+ 
+That assumes that the low/high arguments are PS_TYPE_F64; the vector myData
+can be of any type.  Arguments low/high will be converted to the appropriate
+type and one of the type-specific functions below will be called:
+ 
+    p_psNormalizeVectorRangeU8(myData, low, high)
+    p_psNormalizeVectorRangeU16(myData, low, high)
+    p_psNormalizeVectorRangeU32(myData, low, high)
+    p_psNormalizeVectorRangeU64(myData, low, high)
+    p_psNormalizeVectorRangeS8(myData, low, high)
+    p_psNormalizeVectorRangeS16(myData, low, high)
+    p_psNormalizeVectorRangeS32(myData, low, high)
+    p_psNormalizeVectorRangeS64(myData, low, high)
+    p_psNormalizeVectorRangeF32(myData, low, high)
+    p_psNormalizeVectorRangeF64(myData, low, high)
+ *****************************************************************************/
+#define PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(TYPE) \
+static void p_psNormalizeVectorRange##TYPE( \
+        psVector* myData, \
+        ps##TYPE outLow, \
+        ps##TYPE outHigh) \
+{ \
+    ps##TYPE min = (ps##TYPE) PS_MAX_##TYPE; \
+    ps##TYPE max = (ps##TYPE) -PS_MAX_##TYPE; \
+    psS32 i = 0; \
+    \
+    for (i = 0; i < myData->n; i++) { \
+        if (myData->data.TYPE[i] < min) { \
+            min = myData->data.TYPE[i]; \
+        } \
+        if (myData->data.TYPE[i] > max) { \
+            max = myData->data.TYPE[i]; \
+        } \
+    } \
+    \
+    /* Ensure that max!=min before we divide by (max-min) */ \
+    if (max != min) { \
+        for (i = 0; i < myData->n; i++) { \
+            myData->data.TYPE[i] = (outLow + (myData->data.TYPE[i] - min) * \
+                                    (outHigh - outLow) / (max - min)); \
+        } \
+    } else { \
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: (max==min).  Setting all elements to min.\n"); \
+        for (i = 0; i < myData->n; i++) \
+        { \
+            \
+            myData->data.TYPE[i] = outLow; \
+            \
+        } \
+    } \
+} \
+
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U8)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U16)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U32)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U64)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S8)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S16)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S32)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S64)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F32)
+PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F64)
+
+psS32 p_psNormalizeVectorRange(psVector* myData,
+                               psF64 outLow,
+                               psF64 outHigh)
+{
+    PS_ASSERT_VECTOR_NON_NULL(myData, -1);
+    char* strType;
+
+    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
+    switch (myData->type.type) {
+    case PS_TYPE_U8:
+        p_psNormalizeVectorRangeU8(myData, (psU8) outLow, (psU8) outHigh);
+        break;
+    case PS_TYPE_U16:
+        p_psNormalizeVectorRangeU16(myData, (psU16) outLow, (psU16) outHigh);
+        break;
+    case PS_TYPE_U32:
+        p_psNormalizeVectorRangeU32(myData, (psU32) outLow, (psU32) outHigh);
+        break;
+    case PS_TYPE_U64:
+        p_psNormalizeVectorRangeU64(myData, (psU64) outLow, (psU64) outHigh);
+        break;
+    case PS_TYPE_S8:
+        p_psNormalizeVectorRangeS8(myData, (psS8) outLow, (psS8) outHigh);
+        break;
+    case PS_TYPE_S16:
+        p_psNormalizeVectorRangeS16(myData, (psS16) outLow, (psS16) outHigh);
+        break;
+    case PS_TYPE_S32:
+        p_psNormalizeVectorRangeS32(myData, (psS32) outLow, (psS32) outHigh);
+        break;
+    case PS_TYPE_S64:
+        p_psNormalizeVectorRangeS64(myData, (psS64) outLow, (psS64) outHigh);
+        break;
+    case PS_TYPE_F32:
+        p_psNormalizeVectorRangeF32(myData, (psF32) outLow, (psF32) outHigh);
+        break;
+    case PS_TYPE_F64:
+        p_psNormalizeVectorRangeF64(myData, (psF64) outLow, (psF64) outHigh);
+        break;
+    case PS_TYPE_C32:
+    case PS_TYPE_C64:
+    default:
+        PS_TYPE_NAME(strType, myData->type.type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE, "psLib type %s is not supported.", strType);
+        break;
+    }
+    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
+    return(0);
+}
+
+
Index: /trunk/psLib/src/math/psMathUtils.h
===================================================================
--- /trunk/psLib/src/math/psMathUtils.h	(revision 6305)
+++ /trunk/psLib/src/math/psMathUtils.h	(revision 6305)
@@ -0,0 +1,63 @@
+/** @file psMathUtils.h
+ *  @brief Standard Mathematical Functions.
+ *  @ingroup GROUP00
+ *
+ *  This file contains standard math rotines.
+ *
+ *  @author GLG, MHPCC
+ *
+ *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-02 21:09:07 $
+ *
+ *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
+ */
+
+#ifndef PS_MATHUTILS_H
+#define PS_MATHUTILS_H
+
+#include <stdio.h>
+#include <stdbool.h>
+#include <float.h>
+#include <math.h>
+#include "psVector.h"
+#include "psScalar.h"
+#include "psPolynomial.h"
+
+/** \addtogroup GROUP00
+ *  \{
+ */
+
+/** Performs a binary disection on a monotonically non-decreasing vector.
+ *  Searches through an array of data for a specified value.
+ *
+ *  @return psS32    corresponding index number of specified value
+ */
+psS32 p_psVectorBinDisect(
+    psVector *bins,                    ///< Array of non-decreasing values
+    psScalar *x                        ///< Target value to find
+);
+
+/** Interpolates a series of data points for evaluation at a specific coordinate.  Uses a
+ *  Lagrange interpolation method.
+ *
+ *  @return psScalar*    Lagrange interpolation value at given location
+ */
+psScalar *p_psVectorInterpolate(
+    psScalar *out,
+    psVector *domain,                  ///< Domain (x coords) for interpolation
+    psVector *range,                   ///< Range (y coords) for interpolation
+    psS32 order,                       ///< Order of interpolation function
+    psScalar *x                        ///< Location at which to evaluate
+);
+
+// XXX: Create a single, generic, version of the vector normalize function.
+// XXX: Ask IfA for a public psLib function.
+
+psS32 p_psNormalizeVectorRange(psVector* myData,
+                               psF64 outLow,
+                               psF64 outHigh);
+
+/** \} */ // End of MathGroup Functions
+
+#endif // #ifndef PS_MATHUTILS_H
+
Index: /trunk/psLib/src/math/psMinimizePolyFit.c
===================================================================
--- /trunk/psLib/src/math/psMinimizePolyFit.c	(revision 6304)
+++ /trunk/psLib/src/math/psMinimizePolyFit.c	(revision 6305)
@@ -10,6 +10,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-27 20:08:58 $
+ *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-02 21:09:07 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -34,4 +34,5 @@
 #include "psBinaryOp.h"
 #include "psLogMsg.h"
+#include "psMathUtils.h"
 /*****************************************************************************/
 /* DEFINE STATEMENTS                                                         */
@@ -684,6 +685,5 @@
             f->data.F64[i] = y->data.F64[x->n-1];
         } else {
-            fScalar = p_psVectorInterpolate((psVector *) x, (psVector *) y,
-                                            3, &tmpScalar);
+            fScalar = p_psVectorInterpolate(NULL, (psVector *) x, (psVector *) y, 3, &tmpScalar);
             if (fScalar == NULL) {
                 psError(PS_ERR_UNKNOWN, false, "Could not perform vector interpolation.  Returning NULL.\n");
Index: /trunk/psLib/src/math/psPolynomial.c
===================================================================
--- /trunk/psLib/src/math/psPolynomial.c	(revision 6304)
+++ /trunk/psLib/src/math/psPolynomial.c	(revision 6305)
@@ -7,6 +7,6 @@
 *  polynomials.  It also contains a Gaussian functions.
 *
-*  @version $Revision: 1.140 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-01-26 21:10:22 $
+*  @version $Revision: 1.141 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-02-02 21:09:07 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -286,6 +286,5 @@
 
         for (i=nTerms-3;i>=1;i--) {
-            d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) -
-                             (d->data.F64[i+2]);
+            d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) - (d->data.F64[i+2]);
             if(poly->mask[i] == 0) {
                 d->data.F64[i] += poly->coeff[i];
@@ -293,6 +292,5 @@
         }
 
-        tmp = (x * d->data.F64[1]) -
-              (d->data.F64[2]);
+        tmp = (x * d->data.F64[1]) - (d->data.F64[2]);
         if(poly->mask[0] == 0) {
             tmp += (0.5 * poly->coeff[0]);
Index: /trunk/psLib/src/math/psSpline.c
===================================================================
--- /trunk/psLib/src/math/psSpline.c	(revision 6304)
+++ /trunk/psLib/src/math/psSpline.c	(revision 6305)
@@ -4,9 +4,8 @@
 *         and evaluation routines.
 *
-*  This file will hold the functions for allocated, freeing, and evaluating
-*  splines.
+*  This file contains the routines that allocate, free, and evaluate splines.
 *
-*  @version $Revision: 1.134 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-01-26 21:10:22 $
+*  @version $Revision: 1.135 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-02-02 21:09:08 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -31,4 +30,5 @@
 #include "psConstants.h"
 #include "psErrorText.h"
+#include "psMathUtils.h"
 
 /*****************************************************************************/
@@ -72,4 +72,25 @@
 
     return;
+}
+
+
+static void PS_POLY1D_PRINT(
+    psPolynomial1D *poly)
+{
+    printf("-------------- PS_POLY1D_PRINT() --------------\n");
+    printf("poly->nX is %d\n", poly->nX);
+    for (psS32 i = 0 ; i < (1 + poly->nX) ; i++) {
+        printf("poly->coeff[%d] is %f\n", i, poly->coeff[i]);
+    }
+}
+
+static void PS_PRINT_SPLINE2(psSpline1D *mySpline)
+{
+    printf("-------------- PS_PRINT_SPLINE2() --------------\n");
+    printf("mySpline->n is %d\n", mySpline->n);
+    for (psS32 i = 0 ; i < mySpline->n ; i++) {
+        PS_POLY1D_PRINT(mySpline->spline[i]);
+    }
+    PS_VECTOR_PRINT_F32(mySpline->knots);
 }
 
@@ -139,310 +160,4 @@
     psTrace(__func__, 4, "---- %s() end ----\n", __func__);
     return(derivs2);
-}
-
-
-/*****************************************************************************
-vectorBinDisectTYPE(): This is a macro for a private function which takes as
-input a vector an array of data as well as a single value for that data.  The
-input vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all
-i).  This routine does a binary disection of the vector and returns "i" such
-that (v[i] <= x <= v[i+1).  If x lies outside the range of v[], then this
-routine prints a warning message and returns (-2 or -1).
- *****************************************************************************/
-#define FUNC_MACRO_VECTOR_BIN_DISECT(TYPE) \
-static psS32 vectorBinDisect##TYPE(ps##TYPE *bins, \
-                                   psS32 numBins, \
-                                   ps##TYPE x) \
-{ \
-    psS32 min; \
-    psS32 max; \
-    psS32 mid; \
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__); \
-    if (x < bins[0]) { \
-        psLogMsg(__func__, PS_LOG_WARN, \
-                 "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \
-                 #TYPE, x, bins[0], bins[numBins-1]); \
-        return(-2); \
-    } \
-    if (x > bins[numBins-1]) { \
-        psLogMsg(__func__, PS_LOG_WARN, \
-                 "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \
-                 #TYPE, x, bins[0], bins[numBins-1]); \
-        return(-1); \
-    } \
-    \
-    min = 0; \
-    max = numBins-2; \
-    mid = ((max+1)-min)/2; \
-    while (min != max) { \
-        psTrace(__func__, 6, "(min, mid, max) is (%u, %u, %u): (x, bins) is (%f, %f)\n", min, mid, max, x, bins[mid]); \
-        \
-        if (x == bins[mid]) { \
-            psTrace(__func__, 4, "---- %s(%d) end ----\n", __func__, mid); \
-            return(mid); \
-        } else if (x < bins[mid]) { \
-            max = mid-1; \
-        } else { \
-            min = mid; \
-        } \
-        mid = ((max+1)+min)/2; \
-    } \
-    psTrace(__func__, 4, "---- %s(%d) end ----\n", __func__, min); \
-    return(min); \
-} \
-
-FUNC_MACRO_VECTOR_BIN_DISECT(S8)
-FUNC_MACRO_VECTOR_BIN_DISECT(S16)
-FUNC_MACRO_VECTOR_BIN_DISECT(S32)
-FUNC_MACRO_VECTOR_BIN_DISECT(S64)
-FUNC_MACRO_VECTOR_BIN_DISECT(U8)
-FUNC_MACRO_VECTOR_BIN_DISECT(U16)
-FUNC_MACRO_VECTOR_BIN_DISECT(U32)
-FUNC_MACRO_VECTOR_BIN_DISECT(U64)
-FUNC_MACRO_VECTOR_BIN_DISECT(F32)
-FUNC_MACRO_VECTOR_BIN_DISECT(F64)
-
-/*****************************************************************************
-p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect().
-  *****************************************************************************/
-psS32 p_psVectorBinDisect(
-    psVector *bins,
-    psScalar *x)
-{
-    PS_ASSERT_VECTOR_NON_NULL(bins, -4);
-    PS_ASSERT_VECTOR_NON_EMPTY(bins, -4);
-    PS_ASSERT_PTR_NON_NULL(x, -6);
-    PS_ASSERT_PTR_TYPE_EQUAL(x, bins, -3);
-    char* strType;
-
-    switch (x->type.type) {
-    case PS_TYPE_U8:
-        return(vectorBinDisectU8(bins->data.U8, bins->n, x->data.U8));
-    case PS_TYPE_U16:
-        return(vectorBinDisectU16(bins->data.U16, bins->n, x->data.U16));
-    case PS_TYPE_U32:
-        return(vectorBinDisectU32(bins->data.U32, bins->n, x->data.U32));
-    case PS_TYPE_U64:
-        return(vectorBinDisectU64(bins->data.U64, bins->n, x->data.U64));
-    case PS_TYPE_S8:
-        return(vectorBinDisectS8(bins->data.S8, bins->n, x->data.S8));
-    case PS_TYPE_S16:
-        return(vectorBinDisectS16(bins->data.S16, bins->n, x->data.S16));
-    case PS_TYPE_S32:
-        return(vectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));
-    case PS_TYPE_S64:
-        return(vectorBinDisectS64(bins->data.S64, bins->n, x->data.S64));
-    case PS_TYPE_F32:
-        return(vectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));
-    case PS_TYPE_F64:
-        return(vectorBinDisectF64(bins->data.F64, bins->n, x->data.F64));
-    case PS_TYPE_C32:
-        PS_TYPE_NAME(strType,x->type.type);
-        psError(PS_ERR_BAD_PARAMETER_TYPE,
-                PS_ERRORTEXT_psSpline_TYPE_NOT_SUPPORTED,
-                strType);
-        return 0;
-    case PS_TYPE_C64:
-        PS_TYPE_NAME(strType,x->type.type);
-        psError(PS_ERR_BAD_PARAMETER_TYPE,
-                PS_ERRORTEXT_psSpline_TYPE_NOT_SUPPORTED,
-                strType);
-        return 0;
-    case PS_TYPE_BOOL:
-        PS_TYPE_NAME(strType,x->type.type);
-        psError(PS_ERR_BAD_PARAMETER_TYPE,
-                PS_ERRORTEXT_psSpline_TYPE_NOT_SUPPORTED,
-                strType);
-        return 0;
-    }
-    return(-3);
-}
-
-/*****************************************************************************
-fullInterpolate1DF32(): This routine will take as input n-element floating
-point arrays domain and range, and the x value, assumed to lie with the
-domain vector.  It produces as output the (n-1)-order LaGrange interpolated
-value of x.
- 
-XXX: do we error check for non-distinct domain values?
- 
-XXX: This is a somewhat general function.  There's no reason to put it in this
-file.
- 
-XXX: Sort all these interpolation functions out.  Why are there so many?
- 
-XXX: The fullInterpolate1D macros are not used anywhere.
- *****************************************************************************/
-#define FUNC_MACRO_FULL_INTERPOLATE_1D(TYPE) \
-static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \
-                                     ps##TYPE *range, \
-                                     unsigned int n, \
-                                     ps##TYPE x) \
-{ \
-    \
-    unsigned int i; \
-    unsigned int m; \
-    static psVector *p = NULL; \
-    p = psVectorRecycle(p, n, PS_TYPE_##TYPE); \
-    p_psMemSetPersistent(p, true); \
-    p_psMemSetPersistent(p->data.TYPE, true); \
-    \
-    psTrace(__func__, 4, "---- %s() begin %u-order at x=%f) (%d data points) ----\n", __func__, n-1, x, n); \
-    for (i=0;i<n;i++) { \
-        psTrace(__func__, 6, "domain/range is (%f %f)\n", domain[i], range[i]); \
-    } \
-    \
-    for (i=0;i<n;i++) { \
-        p->data.TYPE[i] = range[i]; \
-        psTrace(__func__, 6, "p->data.TYPE[%u] is %f\n", i, p->data.TYPE[i]); \
-        \
-    } \
-    \
-    /* From NR, during each iteration of the m loop, we are computing the \
-       p_{i ... i+m} terms. \
-    */ \
-    for (m=1;m<n;m++) { \
-        for (i=0;i<n-m;i++) { \
-            /* From NR: we are computing P_{i ... i+m} \
-             */ \
-            p->data.TYPE[i] = (((x-domain[i+m]) * p->data.TYPE[i]) + \
-                               ((domain[i]-x) * p->data.TYPE[i+1])) / \
-                              (domain[i] - domain[i+m]); \
-            psTrace(__func__, 6, "p->data.TYPE[%u] is %f\n", i, p->data.TYPE[i]); \
-        } \
-    } \
-    psTrace(__func__, 4, "---- %s(....) end ----\n", __func__); \
-    return(p->data.TYPE[0]); \
-} \
-
-/* XXX: Do this correctly.
-FUNC_MACRO_FULL_INTERPOLATE_1D(U8)
-FUNC_MACRO_FULL_INTERPOLATE_1D(U16)
-FUNC_MACRO_FULL_INTERPOLATE_1D(U32)
-FUNC_MACRO_FULL_INTERPOLATE_1D(U64)
-FUNC_MACRO_FULL_INTERPOLATE_1D(S8)
-FUNC_MACRO_FULL_INTERPOLATE_1D(S16)
-FUNC_MACRO_FULL_INTERPOLATE_1D(S32)
-FUNC_MACRO_FULL_INTERPOLATE_1D(S64)
-FUNC_MACRO_FULL_INTERPOLATE_1D(F64)
-*/
-FUNC_MACRO_FULL_INTERPOLATE_1D(F32)
-
-
-/*****************************************************************************
-interpolate1DF32(): this is the base 1-D flat memory routine to perform
-LaGrange interpolation.
- *****************************************************************************/
-static psF32 interpolate1DF32(
-    psF32 *domain,
-    psF32 *range,
-    unsigned int n,
-    unsigned int order,
-    psF32 x)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_PTR_NON_NULL(domain, NAN)
-    PS_ASSERT_PTR_NON_NULL(range, NAN)
-    // XXX: Check valid values for n, order, and x?
-
-    psS32 binNum;
-    psS32 numIntPoints = order+1;
-    psS32 origin;
-
-    binNum = vectorBinDisectF32(domain, n, x);
-
-    if (0 == numIntPoints%2) {
-        origin = binNum - ((numIntPoints/2) - 1);
-    } else {
-        origin = binNum - (numIntPoints/2);
-        if ((x-domain[binNum]) > (domain[binNum+1]-x)) {
-            // x is closer to binNum+1.
-            origin = 1 + (binNum - (numIntPoints/2));
-        }
-    }
-    if (origin < 0) {
-        origin = 0;
-    }
-    if ((origin + numIntPoints) > n) {
-        origin = n - numIntPoints;
-    }
-
-    psTrace(__func__, 4, "---- %s(....) end ----\n", __func__);
-    return(fullInterpolate1DF32(&domain[origin], &range[origin], order+1, x));
-}
-
-
-/*****************************************************************************
-p_psVectorInterpolate(): This routine will take as input psVectors domain and
-range, and the x value, assumed to lie with the domain vector.  It produces
-as output the LaGrange interpolated value of a polynomial of the specified
-order around the point x.
- 
-XXX: This stuff does not currently work with a mask.
- 
-XXX: add another psScalar argument for the result (so you don't have to
-     allocate it each time).
- 
-XXX: The VectorCopy routines seg fault when I declare range32 as static.
- *****************************************************************************/
-psScalar *p_psVectorInterpolate(
-    psVector *domain,
-    psVector *range,
-    psS32 order,
-    psScalar *x)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_VECTOR_NON_NULL(domain, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(range, NULL);
-    PS_ASSERT_PTR_NON_NULL(x, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(order, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(domain, range, NULL);
-    PS_ASSERT_PTR_TYPE_EQUAL(domain, range, NULL);
-    PS_ASSERT_PTR_TYPE_EQUAL(domain, x, NULL);
-
-    psVector *range32 = NULL;
-    psVector *domain32 = NULL;
-    if (order > (domain->n - 1)) {
-        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
-                PS_ERRORTEXT_psSpline_NOT_ENOUGH_DATAPOINTS,
-                order);
-        return(NULL);
-    }
-
-    if (x->type.type == PS_TYPE_F32) {
-        psTrace(__func__, 4, "---- %s(NULL) end ----\n", __func__);
-        return(psScalarAlloc(interpolate1DF32(domain->data.F32,
-                                              range->data.F32,
-                                              domain->n,
-                                              order,
-                                              x->data.F32), PS_TYPE_F32));
-    } else if (x->type.type == PS_TYPE_F64) {
-        // XXX: use recycled vectors here.
-        range32 = psVectorCopy(range32, range, PS_TYPE_F32);
-        domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32);
-
-        psScalar *tmpScalar = psScalarAlloc((psF64)
-                                            interpolate1DF32(domain32->data.F32,
-                                                             range32->data.F32,
-                                                             domain32->n,
-                                                             order,
-                                                             (psF32) x->data.F64), PS_TYPE_F64);
-        psFree(range32);
-        psFree(domain32);
-
-        psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-        // XXX: Convert data type to F64?
-        return(tmpScalar);
-
-    } else {
-        char* strType;
-        PS_TYPE_NAME(strType,x->type.type);
-        psError(PS_ERR_BAD_PARAMETER_TYPE,
-                PS_ERRORTEXT_psSpline_TYPE_NOT_SUPPORTED,
-                strType);
-    }
-
-    psTrace(__func__, 4, "---- %s(NULL) end ----\n", __func__);
-    return(NULL);
 }
 
@@ -627,59 +342,4 @@
 }
 
-void PS_POLY1D_PRINT(
-    psPolynomial1D *poly)
-{
-    printf("-------------- PS_POLY1D_PRINT() --------------\n");
-    printf("poly->nX is %d\n", poly->nX);
-    for (psS32 i = 0 ; i < (1 + poly->nX) ; i++) {
-        printf("poly->coeff[%d] is %f\n", i, poly->coeff[i]);
-    }
-}
-
-void PS_PRINT_SPLINE2(psSpline1D *mySpline)
-{
-    printf("-------------- PS_PRINT_SPLINE2() --------------\n");
-    printf("mySpline->n is %d\n", mySpline->n);
-    for (psS32 i = 0 ; i < mySpline->n ; i++) {
-        PS_POLY1D_PRINT(mySpline->spline[i]);
-    }
-    PS_VECTOR_PRINT_F32(mySpline->knots);
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
 
 /*****************************************************************************
@@ -706,5 +366,9 @@
 
     psS32 n = spline->n;
-    psS32 binNum = vectorBinDisectF32(spline->knots->data.F32, (spline->n)+1, x);
+    psScalar tmpScalar;
+    tmpScalar.type.type = PS_TYPE_F32;
+    tmpScalar.data.F32 = x;
+    psS32 binNum = p_psVectorBinDisect(spline->knots, &tmpScalar);
+
     if (binNum < 0) {
         //
@@ -746,5 +410,5 @@
     if (psTraceGetLevel(__func__) >= 6) {
         PS_VECTOR_PRINT_F32(x);
-        // XXX: print spline
+        PS_PRINT_SPLINE2((psSpline1D *) spline);
     }
 
Index: /trunk/psLib/src/math/psSpline.h
===================================================================
--- /trunk/psLib/src/math/psSpline.h	(revision 6304)
+++ /trunk/psLib/src/math/psSpline.h	(revision 6305)
@@ -12,6 +12,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.58 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-11-15 20:10:32 $
+ *  @version $Revision: 1.59 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-02 21:09:08 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -92,26 +92,4 @@
 );
 
-/** Performs a binary disection on a given vector.
- *  Searches through an array of data for a specified value.
- *
- *  @return psS32    corresponding index number of specified value
- */
-psS32 p_psVectorBinDisect(
-    psVector *bins,                    ///< Array of non-decreasing values
-    psScalar *x                        ///< Target value to find
-);
-
-/** Interpolates a series of data points for evaluation at a specific coordinate.  Uses a
- *  Lagrange interpolation method.
- *
- *  @return psScalar*    Lagrange interpolation value at given location
- */
-psScalar *p_psVectorInterpolate(
-    psVector *domain,                  ///< Domain (x coords) for interpolation
-    psVector *range,                   ///< Range (y coords) for interpolation
-    psS32 order,                       ///< Order of interpolation function
-    psScalar *x                        ///< Location at which to evaluate
-);
-
 /** \} */ // End of MathGroup Functions
 
Index: /trunk/psLib/src/math/psStats.c
===================================================================
--- /trunk/psLib/src/math/psStats.c	(revision 6304)
+++ /trunk/psLib/src/math/psStats.c	(revision 6305)
@@ -3,5 +3,5 @@
  *  @ingroup Stat
  *
- *  This file will hold the definition of the histogram and stats data
+ *  This file will hold the definitions of the histogram and stats data
  *  structures.  It also contains prototypes for procedures which operate
  *  on those data structures.
@@ -10,18 +10,12 @@
  *
  *  XXX: The following stats members are never used, or set in this code.
- *      stats->robustN50
  *      stats->clippedNvalues
- *      stats->binsize
  *
  * XXX: Must do
  * nSubsample points
  * use ->min and ->max (PS_STAT_USE_RANGE)
- * use ->binsize (PS_STAT_USE_BINSIZE)
  *
- *
- *
- *
- *  @version $Revision: 1.164 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-01 02:07:24 $
+ *  @version $Revision: 1.165 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-02 21:09:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -49,4 +43,5 @@
 #include "psPolynomial.h"
 #include "psConstants.h"
+#include "psMathUtils.h"
 
 #include "psErrorText.h"
@@ -61,5 +56,5 @@
 #define PS_CLIPPED_SIGMA_LB 1.0
 #define PS_CLIPPED_SIGMA_UB 10.0
-#define PS_POLY_MEDIAN_MAX_ITERATIONS 10
+#define PS_POLY_MEDIAN_MAX_ITERATIONS 30
 
 #define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
@@ -847,6 +842,4 @@
 
     // Calculate the quartile points exactly.
-    // XXX: We should probably do a bin-midpoint if the quartile is not an
-    // integer.
     stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)];
     stats->sampleLQ = sortedVector->data.F32[nValues / 4];
@@ -859,102 +852,4 @@
 }
 
-/******************************************************************************
-p_psVectorSampleStdevOLD(myVector, maskVector, maskVal, stats): calculates the
-stdev of the input vector.
-Inputs
-    myVector
-    maskVector
-    maskVal
-    stats
-Returns
-    NULL
- 
-XXX: remove this
- *****************************************************************************/
-void p_psVectorSampleStdevOLD(const psVector* myVector,
-                              const psVector* errors,
-                              const psVector* maskVector,
-                              psU32 maskVal,
-                              psStats* stats)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    psS32 i = 0;                  // Loop index variable
-    psS32 countInt = 0;           // # of data points being used
-    psF32 countFloat = 0.0;     // # of data points being used
-    psF32 mean = 0.0;           // The mean
-    psF32 diff = 0.0;           // Used in calculating stdev
-    psF32 sumSquares = 0.0;     // temporary variable
-    psF32 sumDiffs = 0.0;       // temporary variable
-
-    // This procedure requires the mean.  If it has not been already
-    // calculated, then call p_psVectorSampleMean()
-    if (0 != isnan(stats->sampleMean)) {
-        p_psVectorSampleMean(myVector, errors, maskVector, maskVal, stats);
-    }
-    // If the mean is NAN, then generate a warning and set the stdev to NAN.
-    if (0 != isnan(stats->sampleMean)) {
-        stats->sampleStdev = NAN;
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): p_psVectorSampleMean() reported a NAN mean.\n");
-        psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-        return;
-    }
-
-    mean = stats->sampleMean;
-    if (stats->options & PS_STAT_USE_RANGE) {
-        if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i]) &&
-                        (stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    diff = myVector->data.F32[i] - mean;
-                    sumSquares += (diff * diff);
-                    sumDiffs += diff;
-                    countInt++;
-                }
-            }
-        } else {
-            for (i = 0; i < myVector->n; i++) {
-                if ((stats->min <= myVector->data.F32[i]) &&
-                        (myVector->data.F32[i] <= stats->max)) {
-                    diff = myVector->data.F32[i] - mean;
-                    sumSquares += (diff * diff);
-                    sumDiffs += diff;
-                    countInt++;
-                }
-            }
-            countInt = myVector->n;
-        }
-    } else {
-        if (maskVector != NULL) {
-            for (i = 0; i < myVector->n; i++) {
-                if (!(maskVal & maskVector->data.U8[i])) {
-                    diff = myVector->data.F32[i] - mean;
-                    sumSquares += (diff * diff);
-                    sumDiffs += diff;
-                    countInt++;
-                }
-            }
-        } else {
-            for (i = 0; i < myVector->n; i++) {
-                diff = myVector->data.F32[i] - mean;
-                sumSquares += (diff * diff);
-                sumDiffs += diff;
-                countInt++;
-            }
-            countInt = myVector->n;
-        }
-    }
-    if (countInt == 0) {
-        stats->sampleStdev = NAN;
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): no valid psVector elements (%d).  Setting stats->sampleStdev = NAN.\n", countInt);
-    } else if (countInt == 1) {
-        stats->sampleStdev = 0.0;
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): only one valid psVector elements (%d).  Setting stats->sampleStdev = 0.0.\n", countInt);
-    } else {
-        countFloat = (psF32)countInt;
-        stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
-    }
-    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-}
 /******************************************************************************
 p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
@@ -1094,9 +989,5 @@
     -2: warning
  
-XXX: Do we really need to calculate median and mean?
- 
 XXX: Use static vectors for tmpMask.
- 
-XXX: This has not been tested.
  *****************************************************************************/
 psS32 p_psVectorClippedStats(const psVector* myVector,
@@ -1223,7 +1114,9 @@
         // Otherwise, use the new results and continue.
         if (isnan(statsTmp->sampleMean) || isnan(statsTmp->sampleStdev)) {
-            // Exit loop.  XXX: throw an error/warning here?
+            // Exit loop.
             iter = stats->clipIter;
-            rc = -1;
+            psError(PS_ERR_UNKNOWN, true, "p_psVectorSampleMean() or p_psVectorSampleStdev() returned a NAN.\n");
+            psFree(tmpMask);
+            return(-1);
         } else {
             clippedMean = statsTmp->sampleMean;
@@ -1246,116 +1139,6 @@
 }
 
-/*****************************************************************************
-These macros and functions define the following functions:
- 
-    p_psNormalizeVectorRange(myData, low, high)
- 
-That assumes that the low/high arguments are PS_TYPE_F64; the vector myData
-can be of any type.  Arguments low/high will be converted to the appropriate
-type and one of the type-specific functions below will be called:
- 
-    p_psNormalizeVectorRangeU8(myData, low, high)
-    p_psNormalizeVectorRangeU16(myData, low, high)
-    p_psNormalizeVectorRangeU32(myData, low, high)
-    p_psNormalizeVectorRangeU64(myData, low, high)
-    p_psNormalizeVectorRangeS8(myData, low, high)
-    p_psNormalizeVectorRangeS16(myData, low, high)
-    p_psNormalizeVectorRangeS32(myData, low, high)
-    p_psNormalizeVectorRangeS64(myData, low, high)
-    p_psNormalizeVectorRangeF32(myData, low, high)
-    p_psNormalizeVectorRangeF64(myData, low, high)
- *****************************************************************************/
-#define PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(TYPE) \
-void p_psNormalizeVectorRange##TYPE(psVector* myData, \
-                                    ps##TYPE outLow, \
-                                    ps##TYPE outHigh) \
-{ \
-    ps##TYPE min = (ps##TYPE) PS_MAX_##TYPE; \
-    ps##TYPE max = (ps##TYPE) -PS_MAX_##TYPE; \
-    psS32 i = 0; \
-    \
-    for (i = 0; i < myData->n; i++) { \
-        if (myData->data.TYPE[i] < min) { \
-            min = myData->data.TYPE[i]; \
-        } \
-        if (myData->data.TYPE[i] > max) { \
-            max = myData->data.TYPE[i]; \
-        } \
-    } \
-    \
-    /* Ensure that max!=min before we divide by (max-min) */ \
-    if (max != min) { \
-        for (i = 0; i < myData->n; i++) { \
-            myData->data.TYPE[i] = (outLow + (myData->data.TYPE[i] - min) * \
-                                    (outHigh - outLow) / (max - min)); \
-        } \
-    } else { \
-        psLogMsg(__func__, PS_LOG_WARN, "WARNING: (max==min).  Setting all elements to min.\n"); \
-        for (i = 0; i < myData->n; i++) \
-        { \
-            \
-            myData->data.TYPE[i] = outLow; \
-            \
-        } \
-    } \
-} \
-
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U8)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U16)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U32)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U64)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S8)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S16)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S32)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S64)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F32)
-PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F64)
-
-void p_psNormalizeVectorRange(psVector* myData,
-                              psF64 outLow,
-                              psF64 outHigh)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    switch (myData->type.type) {
-    case PS_TYPE_U8:
-        p_psNormalizeVectorRangeU8(myData, (psU8) outLow, (psU8) outHigh);
-        break;
-    case PS_TYPE_U16:
-        p_psNormalizeVectorRangeU16(myData, (psU16) outLow, (psU16) outHigh);
-        break;
-    case PS_TYPE_U32:
-        p_psNormalizeVectorRangeU32(myData, (psU32) outLow, (psU32) outHigh);
-        break;
-    case PS_TYPE_U64:
-        p_psNormalizeVectorRangeU64(myData, (psU64) outLow, (psU64) outHigh);
-        break;
-    case PS_TYPE_S8:
-        p_psNormalizeVectorRangeS8(myData, (psS8) outLow, (psS8) outHigh);
-        break;
-    case PS_TYPE_S16:
-        p_psNormalizeVectorRangeS16(myData, (psS16) outLow, (psS16) outHigh);
-        break;
-    case PS_TYPE_S32:
-        p_psNormalizeVectorRangeS32(myData, (psS32) outLow, (psS32) outHigh);
-        break;
-    case PS_TYPE_S64:
-        p_psNormalizeVectorRangeS64(myData, (psS64) outLow, (psS64) outHigh);
-        break;
-    case PS_TYPE_F32:
-        p_psNormalizeVectorRangeF32(myData, (psF32) outLow, (psF32) outHigh);
-        break;
-    case PS_TYPE_F64:
-        p_psNormalizeVectorRangeF64(myData, (psF64) outLow, (psF64) outHigh);
-        break;
-    case PS_TYPE_C32:
-    case PS_TYPE_C64:
-    default:
-        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
-                "Unallowable operation: %s has incorrect type.",
-                myData);
-        break;
-    }
-    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
-}
+
+
 
 /******************************************************************************
@@ -1371,8 +1154,9 @@
 XXX: Create a 2nd-order polynomial version and solve for X analytically.
  *****************************************************************************/
-psF32 p_ps1DPolyMedian(psPolynomial1D* myPoly,
-                       psF32 rangeLow,
-                       psF32 rangeHigh,
-                       psF32 getThisValue)
+psF32 p_ps1DPolyMedian(
+    psPolynomial1D* myPoly,
+    psF32 rangeLow,
+    psF32 rangeHigh,
+    psF32 getThisValue)
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
@@ -1405,4 +1189,5 @@
 
         f = psPolynomial1DEval(myPoly, midpoint);
+        printf("p_ps1DPolyMedian(): f(%f) is %f\n", midpoint, f);
         if (fabs(f - getThisValue) <= FLT_EPSILON) {
             return (midpoint);
@@ -1419,4 +1204,65 @@
     return (midpoint);
 }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#define PS_SQRT(x) sqrt(x)
+static psF32 QuadraticInverse(
+    psF32 a,
+    psF32 b,
+    psF32 c,
+    psF32 y)
+{
+    psF64 tmp = sqrt((y - c)/a + (b*b)/(4.0 * a * a));
+
+    psF64 x1 = -b/(2.0*a) + tmp;
+    psF64 x2 = -b/(2.0*a) - tmp;
+
+    psF64 y1 = (a * x1 * x1) + (b * x1) + c;
+    psF64 y2 = (a * x2 * x2) + (b * x2) + c;
+
+    printf("QuadraticInverse: %fx^2 + %fx + %f\n", a, b, c);
+    printf("QuadraticInverse: y is %f\n", y);
+    printf("QuadraticInverse: (x1, x2) is (%f %f)\n", x1, x2);
+    printf("QuadraticInverse: (y1, y2) is (%f %f)\n", y1, y2);
+
+    return(x1);
+}
+
 
 /******************************************************************************
@@ -1438,8 +1284,9 @@
 arbitrary vectors.  We should probably test that condition.
 *****************************************************************************/
-psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec,
-                                        psVector *yVec,
-                                        psS32 binNum,
-                                        psF32 yVal)
+psF32 fitQuadraticSearchForYThenReturnX(
+    psVector *xVec,
+    psVector *yVec,
+    psS32 binNum,
+    psF32 yVal)
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
@@ -1461,5 +1308,4 @@
     psVector *y = psVectorAlloc(3, PS_TYPE_F64);
     psVector *yErr = psVectorAlloc(3, PS_TYPE_F64);
-    psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
 
     psF32 tmpFloat = 0.0f;
@@ -1489,5 +1335,4 @@
             if (!(y->data.F64[1] <= y->data.F64[2])) {
                 psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n");
-                psFree(myPoly);
                 psFree(x);
                 psFree(y);
@@ -1505,5 +1350,4 @@
             if (!(y->data.F64[1] >= y->data.F64[2])) {
                 psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n");
-                psFree(myPoly);
                 psFree(x);
                 psFree(y);
@@ -1525,5 +1369,5 @@
 
         // Determine the coefficients of the polynomial.
-        //        myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
+        psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
         myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, yErr, x);
         if (myPoly == NULL) {
@@ -1531,5 +1375,4 @@
                     false,
                     PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLYNOMIAL_1D_FIT);
-            psFree(myPoly);
             psFree(x);
             psFree(y);
@@ -1541,5 +1384,6 @@
         psTrace(__func__, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
         psTrace(__func__, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
-        psTrace(__func__, 6, "Fitted y vec is (%f %f %f)\n", (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
+        psTrace(__func__, 6, "Fitted y vec is (%f %f %f)\n",
+                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
                 (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
                 (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
@@ -1548,5 +1392,11 @@
         // polynomial, looking for the value x such that f(x) = yVal
         psTrace(__func__, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
+        printf("Calling p_ps1DPolyMedian() for the y-value (%.2f) that has an x in the range (%.2f - %.2f)\n", yVal, x->data.F64[0], x->data.F64[2]);
         tmpFloat = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], yVal);
+        printf("Cool, tmpFloat is %.2f\n", tmpFloat);
+
+        QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal);
+        psFree(myPoly);
+
 
         if (isnan(tmpFloat)) {
@@ -1554,5 +1404,4 @@
                     false,
                     PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLY_MEDIAN);
-            psFree(myPoly);
             psFree(x);
             psFree(y);
@@ -1581,5 +1430,4 @@
 
     psTrace(__func__, 6, "FIT: return %f\n", tmpFloat);
-    psFree(myPoly);
     psFree(x);
     psFree(y);
@@ -1591,9 +1439,10 @@
 
 /******************************************************************************
-XXX: We assume unnormalized gaussians.
+NOTE: We assume unnormalized gaussians.
  *****************************************************************************/
-psF32 psMinimizeLMChi2Gauss1D(psVector *deriv,
-                              const psVector *params,
-                              const psVector *coords)
+psF32 psMinimizeLMChi2Gauss1D(
+    psVector *deriv,
+    const psVector *params,
+    const psVector *coords)
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
@@ -1627,4 +1476,5 @@
 }
 
+//XXX: Use the psLib function?
 psF32 LinInterpolate(psF32 x0,
                      psF32 x1,
@@ -1643,76 +1493,4 @@
     return(x0 + y * (x1 - x0) / (y1 - y0));
 }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
 
 
@@ -1754,5 +1532,5 @@
 XXX: Review and ensure that all memory is free'ed at premature exits.
 *****************************************************************************/
-#define INITIAL_NUM_BINS 1000.0
+#define INITIAL_NUM_BINS 500.0
 psS32 p_psVectorRobustStats(const psVector* myVector,
                             const psVector* errors,
@@ -1770,6 +1548,6 @@
     psHistogram *cumulativeRobustHistogram = NULL;
     psS32 numBins = 0;
-    psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32);
-    tmpScalar->type.type = PS_TYPE_F32;
+    psScalar tmpScalar;
+    tmpScalar.type.type = PS_TYPE_F32;
     psS32 totalDataPoints = 0;
     psS32 rc = 0;
@@ -1806,5 +1584,4 @@
                 psFree(tmpStatsMinMax);
                 psFree(tmpMaskVec);
-                psFree(tmpScalar);
                 psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
                 return(1);
@@ -1833,5 +1610,4 @@
             psFree(tmpStatsMinMax);
             psFree(tmpMaskVec);
-            psFree(tmpScalar);
 
             psTrace(__func__, 4, "---- %s(0) end  ----\n", __func__);
@@ -1885,6 +1661,6 @@
             binMedian = 0;
         } else {
-            tmpScalar->data.F32 = totalDataPoints/2.0;
-            binMedian = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+            tmpScalar.data.F32 = totalDataPoints/2.0;
+            binMedian = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
             if (binMedian < 0) {
                 psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 50 precent data point (%d).\n", binMedian);
@@ -1892,5 +1668,4 @@
                 psFree(robustHistogram);
                 psFree(cumulativeRobustHistogram);
-                psFree(tmpScalar);
                 psFree(tmpMaskVec);
                 psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
@@ -1905,4 +1680,5 @@
         // XXX: Check return codes.
         //
+        printf("ADD: step 3: Interpolate to the exact 50-percent position: this is the robust histogram median.\n");
         stats->robustMedian = fitQuadraticSearchForYThenReturnX(
                                   *(psVector* *)&cumulativeRobustHistogram->bounds,
@@ -1918,20 +1694,20 @@
         psS32 binLo = 0;
         psS32 binHi = 0;
-        tmpScalar->data.F32 = totalDataPoints * 0.158655f;
-        if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
+        tmpScalar.data.F32 = totalDataPoints * 0.158655f;
+        if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
             binLo = 0;
         } else {
             // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
-            binLo = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+            binLo = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
             if (binLo > cumulativeRobustHistogram->nums->n-1) {
                 binLo = cumulativeRobustHistogram->nums->n-1;
             }
         }
-        tmpScalar->data.F32 = totalDataPoints * 0.841345f;
-        if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
+        tmpScalar.data.F32 = totalDataPoints * 0.841345f;
+        if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
             binHi = 0;
         } else {
             // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
-            binHi = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+            binHi = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
             if (binHi > cumulativeRobustHistogram->nums->n-1) {
                 binHi = cumulativeRobustHistogram->nums->n-1;
@@ -1947,5 +1723,4 @@
             psFree(robustHistogram);
             psFree(cumulativeRobustHistogram);
-            psFree(tmpScalar);
             psFree(tmpMaskVec);
             psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
@@ -2014,17 +1789,6 @@
         if (sigma < (2 * binSize)) {
             psTrace(__func__, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize);
-            psF32 medianLo;
-            psF32 medianHi;
-            if ((binMedian - 25) > 0) {
-                medianLo = robustHistogram->bounds->data.F32[binMedian - 25];
-            } else {
-                medianLo = robustHistogram->bounds->data.F32[0];
-            }
-            if ((binMedian + 25) < robustHistogram->bounds->n) {
-                medianHi = robustHistogram->bounds->data.F32[binMedian + 25];
-            } else {
-                medianHi = robustHistogram->bounds->data.F32[robustHistogram->bounds->n - 1];
-            }
-
+            psF32 medianLo = robustHistogram->bounds->data.F32[PS_MAX(0, (binMedian - 25))];
+            psF32 medianHi = robustHistogram->bounds->data.F32[PS_MIN(robustHistogram->bounds->n - 1, (binMedian + 25))];
             psTrace(__func__, 6, "Masking data more than 25 bins from the median (%.2f)\n", stats->robustMedian);
             psTrace(__func__, 6, "The median is at bin number %d.  We mask bins outside the bin range (%d:%d)\n", binMedian, binMedian - 25, binMedian + 25);
@@ -2052,17 +1816,17 @@
     psS32 binHi25 = 0;
 
-    tmpScalar->data.F32 = totalDataPoints * 0.25f;
-    if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
+    tmpScalar.data.F32 = totalDataPoints * 0.25f;
+    if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
         // XXX: Special case where its in first bin.  Must code last bin.
         binLo25 = 0;
     } else {
-        binLo25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
-    }
-    tmpScalar->data.F32 = totalDataPoints * 0.75f;
-    if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
+        binLo25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
+    }
+    tmpScalar.data.F32 = totalDataPoints * 0.75f;
+    if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
         // XXX: Special case where its in first bin.  Must code last bin.
         binHi25 = 0;
     } else {
-        binHi25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
+        binHi25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
     }
     if ((binLo25 < 0) || (binHi25 < 0)) {
@@ -2071,5 +1835,4 @@
         psFree(robustHistogram);
         psFree(cumulativeRobustHistogram);
-        psFree(tmpScalar);
         psFree(tmpMaskVec);
         psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
@@ -2084,4 +1847,5 @@
     // XXX: Check for errors.
     //
+    printf("ADD: step 8: Interpolate for the lower quartile positions.\n");
     psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(
                            *(psVector* *)&cumulativeRobustHistogram->bounds,
@@ -2089,4 +1853,5 @@
                            binLo25,
                            totalDataPoints * 0.25f);
+    printf("ADD: step 8: Interpolate for the upper quartile positions.\n");
     psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(
                            *(psVector* *)&cumulativeRobustHistogram->bounds,
@@ -2133,5 +1898,4 @@
             psFree(robustHistogram);
             psFree(cumulativeRobustHistogram);
-            psFree(tmpScalar);
             psFree(tmpMaskVec);
             psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
@@ -2170,16 +1934,16 @@
         psS32 binMin = 0;
         psS32 binMax = 0;
-        tmpScalar->data.F32 = stats->robustMedian - (2.0 * sigma);
-        if (tmpScalar->data.F32 <= newHistogram->bounds->data.F32[0]) {
+        tmpScalar.data.F32 = stats->robustMedian - (2.0 * sigma);
+        if (tmpScalar.data.F32 <= newHistogram->bounds->data.F32[0]) {
             binMin = 0;
         } else {
-            binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
-        }
-
-        tmpScalar->data.F32 = stats->robustMedian + (2.0 + sigma);
-        if (tmpScalar->data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) {
+            binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar);
+        }
+
+        tmpScalar.data.F32 = stats->robustMedian + (2.0 + sigma);
+        if (tmpScalar.data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) {
             binMax = newHistogram->bounds->n-1;
         } else {
-            binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
+            binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar);
         }
         if ((binMin < 0) || (binMax < 0)) {
@@ -2207,16 +1971,16 @@
         // histogram median.
         //
-        tmpScalar->data.F32 = stats->robustMedian - (20.0 * sigma);
-        if (tmpScalar->data.F32 < tmpStatsMinMax->min) {
+        tmpScalar.data.F32 = stats->robustMedian - (20.0 * sigma);
+        if (tmpScalar.data.F32 < tmpStatsMinMax->min) {
             binMin = 0;
         } else {
-            binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
+            binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar);
             // XXX: check for errors here.
         }
-        tmpScalar->data.F32 = stats->robustMedian + (20.0 * sigma);
-        if (tmpScalar->data.F32 > tmpStatsMinMax->max) {
+        tmpScalar.data.F32 = stats->robustMedian + (20.0 * sigma);
+        if (tmpScalar.data.F32 > tmpStatsMinMax->max) {
             binMax = newHistogramSmoothed->n - 1;
         } else {
-            binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
+            binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar);
             // XXX: check for errors here.
         }
@@ -2275,5 +2039,4 @@
             psFree(robustHistogram);
             psFree(cumulativeRobustHistogram);
-            psFree(tmpScalar);
             psFree(newHistogram);
             psFree(x);
@@ -2314,5 +2077,4 @@
     psFree(tmpStatsMinMax);
     psFree(cumulativeRobustHistogram);
-    psFree(tmpScalar);
     psFree(tmpMaskVec);
     psFree(robustHistogram);
@@ -2321,12 +2083,4 @@
     return(0);
 }
-
-
-
-
-
-
-
-
 
 /*****************************************************************************/
@@ -2362,5 +2116,5 @@
     newStruct->robustUQ = NAN;
     newStruct->robustLQ = NAN;
-    newStruct->robustN50 = -1;            // XXX: This is never used
+    newStruct->robustN50 = -1;
     newStruct->fittedMean = NAN;
     newStruct->fittedStdev = NAN;
@@ -2373,5 +2127,5 @@
     newStruct->min = NAN;
     newStruct->max = NAN;
-    newStruct->binsize = NAN;          // XXX: This is never used
+    newStruct->binsize = NAN;
     newStruct->nSubsample = 100000;
     newStruct->options = options;
Index: /trunk/psLib/src/math/psStats.h
===================================================================
--- /trunk/psLib/src/math/psStats.h	(revision 6304)
+++ /trunk/psLib/src/math/psStats.h	(revision 6305)
@@ -14,6 +14,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.49 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-26 23:49:11 $
+ *  @version $Revision: 1.50 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-02 21:09:08 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -202,77 +202,4 @@
 );
 
-// XXX: Create a single, generic, version of the vector normalize function.
-// XXX: Ask IfA for a public psLib function.
-
-/** Normalize the range of a vector containing U8 values  */
-void p_psNormalizeVectorRangeU8(
-    psVector* myData,                  ///< Vector containing the U8 values
-    psU8 low,                          ///< Minimum value
-    psU8 high                          ///< Maximum value
-);
-
-/** Normalize the range of a vector containing U16 values  */
-void p_psNormalizeVectorRangeU16(
-    psVector* myData,                  ///< Vector containing the U16 values
-    psU16 low,                         ///< Minimum value
-    psU16 high                         ///< Maximum value
-);
-
-/** Normalize the range of a vector containing U32 values  */
-void p_psNormalizeVectorRangeU32(
-    psVector* myData,                  ///< Vector containing the U32 values
-    psU32 low,                         ///< Minimum value
-    psU32 high                         ///< Maximum value
-);
-
-/** Normalize the range of a vector containing U64 values  */
-void p_psNormalizeVectorRangeU64(
-    psVector* myData,                  ///< Vector containing the U64 values
-    psU64 low,                         ///< Minimum value
-    psU64 high                         ///< Maximum value
-);
-
-/** Normalize the range of a vector containing S8 values  */
-void p_psNormalizeVectorRangeS8(
-    psVector* myData,                  ///< Vector containing the S8 values
-    psS8 low,                          ///< Minimum value
-    psS8 high                          ///< Maximum value
-);
-
-/** Normalize the range of a vector containing S16 values  */
-void p_psNormalizeVectorRangeS16(
-    psVector* myData,                  ///< Vector containing the S16 values
-    psS16 low,                         ///< Minimum value
-    psS16 high                         ///< Maximum value
-);
-
-/** Normalize the range of a vector containing S32 values  */
-void p_psNormalizeVectorRangeS32(
-    psVector* myData,                  ///< Vector containing the S32 values
-    psS32 low,                         ///< Minimum value
-    psS32 high                         ///< Maximum value
-);
-
-/** Normalize the range of a vector containing S64 values  */
-void p_psNormalizeVectorRangeS64(
-    psVector* myData,                  ///< Vector containing the S64 values
-    psS64 low,                         ///< Minimum value
-    psS64 high                         ///< Maximum value
-);
-
-/** Normalize the range of a vector containing F32 values  */
-void p_psNormalizeVectorRangeF32(
-    psVector* myData,                  ///< Vector containing the F32 values
-    psF32 low,                         ///< Minimum value
-    psF32 high                         ///< Maximum value
-);
-
-/** Normalize the range of a vector containing F64 values  */
-void p_psNormalizeVectorRangeF64(
-    psVector* myData,                  ///< Vector containing the F64 values
-    psF64 low,                         ///< Minimum value
-    psF64 high                         ///< Maximum value
-);
-
 /// @}