Index: trunk/psLib/src/math/psMinimize.c
===================================================================
--- trunk/psLib/src/math/psMinimize.c	(revision 1964)
+++ trunk/psLib/src/math/psMinimize.c	(revision 1982)
@@ -7,8 +7,8 @@
  *  fit a 1-D polynomial to a set of data points.
  *
- *  @author George Gusciora, MHPCC
+ *  @author GLF, MHPCC
  *
- *  @version $Revision: 1.54 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-05 23:14:15 $
+ *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-06 21:30:53 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -227,4 +227,9 @@
 supplied, but does not equal the domains specified in mySpline?
 XXX: can psSpline be NULL?
+ 
+XXX: reimplement this assuming that mySpline is NULL?
+ 
+XXX: What happens if X is NULL, then an index vector is generated for X, but
+that index vector lies outside the range vectors in mySpline?
  *****************************************************************************/
 psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline,              ///< The spline which will be generated.
@@ -269,10 +274,10 @@
         CONVERT_VECTOR_F64_TO_F32(x, x32, x32Static);
     }
-    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(x, y);
-    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(yErr, y);
-
-    if (y->n != (1 + mySpline->n)) {
+    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(x32, y32);
+    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(yErr32, y32);
+
+    if (y32->n != (1 + mySpline->n)) {
         psError(__func__, "data size / spline size mismatch (%d %d)\n",
-                y->n, mySpline->n);
+                y32->n, mySpline->n);
         return(NULL);
     }
@@ -282,7 +287,7 @@
     if (2 == (mySpline->spline[0])->n) {
         for (i=0;i<mySpline->n;i++) {
-            slope = (y->data.F32[i+1] - y->data.F32[i]) /
+            slope = (y32->data.F32[i+1] - y32->data.F32[i]) /
                     (mySpline->domains[i+1] - mySpline->domains[i]);
-            (mySpline->spline[i])->coeff[0] = y->data.F32[i] -
+            (mySpline->spline[i])->coeff[0] = y32->data.F32[i] -
                                               (slope * mySpline->domains[i]);
 
@@ -307,6 +312,6 @@
     // If we get here, then we know these are cubic splines.  We first
     // generate the second derivatives at each data point.
-    mySpline->p_psDeriv2 = CalculateSecondDerivs(x, y);
-    for (i=0;i<y->n;i++)
+    mySpline->p_psDeriv2 = CalculateSecondDerivs(x32, y32);
+    for (i=0;i<y32->n;i++)
         psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
                 "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
@@ -316,5 +321,5 @@
     // and Numerical Recipes in C.
     for (i=0;i<numSplines;i++) {
-        H = x->data.F32[i+1] - x->data.F32[i];
+        H = x32->data.F32[i+1] - x32->data.F32[i];
         psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
                 "x data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H);
@@ -323,15 +328,15 @@
         //
         // From (1)
-        (mySpline->spline[i])->coeff[0] = (y->data.F32[i] * x->data.F32[i+1]/H);
+        (mySpline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H);
         // From (2)
-        ((mySpline->spline[i])->coeff[0])-= ((y->data.F32[i+1] * x->data.F32[i])/H);
+        ((mySpline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H);
         // From (3)
-        tmp = (x->data.F32[i+1] * x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H);
-        tmp-= (x->data.F32[i+1] / H);
+        tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
+        tmp-= (x32->data.F32[i+1] / H);
         tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0;
         ((mySpline->spline[i])->coeff[0])+= tmp;
         // From (4)
-        tmp = -(x->data.F32[i] * x->data.F32[i] * x->data.F32[i]) / (H * H * H);
-        tmp+= (x->data.F32[i] / H);
+        tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
+        tmp+= (x32->data.F32[i] / H);
         tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0;
         ((mySpline->spline[i])->coeff[0])+= tmp;
@@ -341,14 +346,14 @@
         //
         // From (1)
-        (mySpline->spline[i])->coeff[1] = -(y->data.F32[i]) / H;
+        (mySpline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H;
         // From (2)
-        ((mySpline->spline[i])->coeff[1])+= (y->data.F32[i+1] / H);
+        ((mySpline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H);
         // From (3)
-        tmp = -3.0 * (x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H);
+        tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
         tmp+= (1.0 / H);
         tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0;
         ((mySpline->spline[i])->coeff[1])+= tmp;
         // From (4)
-        tmp = 3.0 * (x->data.F32[i] * x->data.F32[i]) / (H * H * H);
+        tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
         tmp-= (1.0 / H);
         tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
@@ -359,7 +364,7 @@
         //
         // From (3)
-        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x->data.F32[i+1] / (6.0 * H);
+        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H);
         // From (4)
-        ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x->data.F32[i] / (6.0 * H));
+        ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H));
 
         //