Index: trunk/psLib/src/math/psStats.c
===================================================================
--- trunk/psLib/src/math/psStats.c	(revision 2341)
+++ trunk/psLib/src/math/psStats.c	(revision 2342)
@@ -9,6 +9,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.92 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-11-12 19:29:39 $
+ *  @version $Revision: 1.93 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-11-12 19:36:07 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -925,4 +925,6 @@
  
 XXX: Should we make the arguments psScalars?
+ 
+XXX: Should we make the arguments PS_TYPE_F64?
  *****************************************************************************/
 #define PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(TYPE) \
@@ -983,5 +985,5 @@
 XXX: Terminate when f(x)-getThisValue is within some error tolerance.
  
-XXX: Solve for X analytically.
+XXX: Create a 2nd-order polynomial version and solve for X analytically.
  *****************************************************************************/
 float p_ps1DPolyMedian(psPolynomial1D* myPoly,
@@ -1020,9 +1022,11 @@
 
 /******************************************************************************
-fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): This routine
-takes psVectors of x/y pairs as input, and fits a quadratic to the 3 points
-surrounding element binNum in the vectors (the midpoint between element i
-and i+1 is used for x[i]).  It then determines for what value x does that
-quadratic f(x) = yVal (the input parameter).
+fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): A general
+routine which fits a quadratic to three points and returns the x-value
+corresponding to the input y-value.  This routine takes psVectors of x/y pairs
+as input, and fits a quadratic to the 3 points surrounding element binNum in
+the vectors (the midpoint between element i and i+1 is used for x[i]).  It
+then determines for what value x does that quadratic f(x) = yVal (the input
+parameter).
  
 XXX: After you fit the polynomial, solve for X analytically.
@@ -1035,16 +1039,19 @@
                                         float yVal)
 {
+    PS_VECTOR_CHECK_NULL(xVec, NAN);
+    PS_VECTOR_CHECK_NULL(yVec, NAN);
+    PS_VECTOR_CHECK_TYPE(xVec, PS_TYPE_F64, NAN);
+    PS_VECTOR_CHECK_TYPE(yVec, PS_TYPE_F64, NAN);
+    PS_VECTOR_CHECK_SIZE_EQUAL(xVec, yVec, NAN);
+    PS_INT_CHECK_RANGE(binNum, 0, (xVec->n - 1), NAN);
+
     PS_VECTOR_DECLARE_ALLOC_STATIC(x, 3, PS_TYPE_F64);
     PS_VECTOR_DECLARE_ALLOC_STATIC(y, 3, PS_TYPE_F64);
     PS_VECTOR_DECLARE_ALLOC_STATIC(yErr, 3, PS_TYPE_F64);
-    static psPolynomial1D* myPoly = NULL;
-
-    if (myPoly == NULL) {
-        myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);
-        p_psMemSetPersistent(myPoly, true);
-    }
+    PS_POLY_1D_DECLARE_ALLOC_STATIC(myPoly, 2, PS_POLYNOMIAL_ORD);
     float tmpFloat;
 
     if ((binNum > 0) && (binNum < (yVec->n - 2))) {
+        // The general case.  We have all three points.
         x->data.F64[0] = (double) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));
         x->data.F64[1] = (double) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1]));
@@ -1070,8 +1077,12 @@
         tmpFloat = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], yVal);
     } else {
+        // The special case where we have two points only at the beginning of
+        // the vectors x and y.
         if (binNum == 0) {
             tmpFloat = 0.5 * (xVec->data.F32[binNum] +
                               xVec->data.F32[binNum + 1]);
         } else if (binNum == (xVec->n - 1)) {
+            // The special case where we have two points only at the end of
+            // the vectors x and y.
             // XXX: Is this right?
             tmpFloat = xVec->data.F32[binNum];