Index: trunk/psLib/src/math/psSpline.c
===================================================================
--- trunk/psLib/src/math/psSpline.c	(revision 5155)
+++ trunk/psLib/src/math/psSpline.c	(revision 5179)
@@ -7,6 +7,6 @@
 *  splines.
 *
-*  @version $Revision: 1.129 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-09-27 23:16:59 $
+*  @version $Revision: 1.130 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-09-29 19:40:59 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -491,33 +491,37 @@
     PS_ASSERT_VECTOR_NON_NULL(y, NULL);
     PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
-    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
-    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
     psS32 numSplines = (y->n)-1;
-
     psTrace(__func__, 5, "numSplines is %d\n", numSplines);
     //
     // The following code ensures that xPtr and yPtr points to a psF32 psVector.
     //
-    psVector *xF32 = NULL;
     psVector *xPtr = NULL;
-    if (PS_TYPE_F64 == x->type.type) {
-        xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
-        for (psS32 i = 0 ; i < x->n ; i++) {
-            xF32->data.F32[i] = (psF32) x->data.F64[i];
+    if (x != NULL) {
+        // Convert x to F32 if necessary.
+        PS_ASSERT_VECTOR_NON_NULL(x, NULL);
+        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
+        if (PS_TYPE_F64 == x->type.type) {
+            xPtr = psVectorAlloc(y->n, PS_TYPE_F32);
+            for (psS32 i = 0 ; i < x->n ; i++) {
+                xPtr->data.F32[i] = (psF32) x->data.F64[i];
+            }
+        } else if (PS_TYPE_F32 == x->type.type) {
+            xPtr = (psVector *) x;
         }
-        xPtr = xF32;
-    } else if (PS_TYPE_F32 == x->type.type) {
-        xPtr = (psVector *) x;
-    }
-
-    psVector *yF32 = NULL;
+    } else {
+        // Allocate an index vector for x
+        xPtr = psVectorCreate(NULL, 0.0, (psF64) y->n, 1.0, PS_TYPE_F32);
+        printf("x, y is (%d, %d)\n", (psS32) xPtr->n, (psS32) y->n);
+        xPtr = xPtr;
+    }
+
     psVector *yPtr = NULL;
+    // Convert y to F32 if necessary.
     if (PS_TYPE_F64 == y->type.type) {
-        yF32 = psVectorAlloc(y->n, PS_TYPE_F32);
+        yPtr = psVectorAlloc(y->n, PS_TYPE_F32);
         for (psS32 i = 0 ; i < y->n ; i++) {
-            yF32->data.F32[i] = (psF32) y->data.F64[i];
+            yPtr->data.F32[i] = (psF32) y->data.F64[i];
         }
-        yPtr = yF32;
     } else {
         yPtr = (psVector *) y;
@@ -531,14 +535,7 @@
     spline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
     for (psS32 i=0;i<numSplines;i++) {
-        (spline->spline)[i] = psPolynomial1DAlloc(3, PS_POLYNOMIAL_ORD);
-    }
-    //
-    // XXX: Ensure that the knots are distinct, and monotonic.
-    // XXX: Use a vector dup macro, or function here.
-    //
-    spline->knots = psVectorAlloc(x->n, PS_TYPE_F32);
-    for (psS32 i = 0 ; i < x->n ; i++) {
-        spline->knots->data.F32[i] = xPtr->data.F32[i];
-    }
+        spline->spline[i] = psPolynomial1DAlloc(3, PS_POLYNOMIAL_ORD);
+    }
+    spline->knots = psVectorCopy(NULL, xPtr, PS_TYPE_F32);
     //
     // Generate the second derivatives at each data point.
@@ -554,25 +551,24 @@
         psF32 H = xPtr->data.F32[i+1] - xPtr->data.F32[i];
         if (fabs(H) <= FLT_EPSILON) {
-            printf("XXX: Generate error: x data points are not distinct (%d %d).\n", i, i+1);
+            psError(PS_ERR_UNKNOWN, false, "x data points are not distinct (%d %d).\n", i, i+1);
         }
         psTrace(__func__, 6, "x data (%f - %f) (%f)\n", xPtr->data.F32[i], xPtr->data.F32[i+1], H);
-
         //
         // ******** Calculate 0-order term ********
         //
         // From (1)
-        (spline->spline[i])->coeff[0] = (yPtr->data.F32[i] * xPtr->data.F32[i+1]/H);
+        spline->spline[i]->coeff[0] = yPtr->data.F32[i] * xPtr->data.F32[i+1]/H;
         // From (2)
-        ((spline->spline[i])->coeff[0])-= ((yPtr->data.F32[i+1] * xPtr->data.F32[i])/H);
+        spline->spline[i]->coeff[0]-= (yPtr->data.F32[i+1] * xPtr->data.F32[i])/H;
         // From (3)
         psF32 tmp = (xPtr->data.F32[i+1] * xPtr->data.F32[i+1] * xPtr->data.F32[i+1]) / (H * H * H);
-        tmp-= (xPtr->data.F32[i+1] / H);
-        tmp*= (spline->p_psDeriv2)[i] * H * H / 6.0;
-        ((spline->spline[i])->coeff[0])+= tmp;
+        tmp-= xPtr->data.F32[i+1] / H;
+        tmp*= spline->p_psDeriv2[i] * H * H / 6.0;
+        spline->spline[i]->coeff[0]+= tmp;
         // From (4)
         tmp = -(xPtr->data.F32[i] * xPtr->data.F32[i] * xPtr->data.F32[i]) / (H * H * H);
-        tmp+= (xPtr->data.F32[i] / H);
-        tmp*= (spline->p_psDeriv2)[i+1] * H * H / 6.0;
-        ((spline->spline[i])->coeff[0])+= tmp;
+        tmp+= xPtr->data.F32[i] / H;
+        tmp*= spline->p_psDeriv2[i+1] * H * H / 6.0;
+        (spline->spline[i]->coeff[0])+= tmp;
 
         //
@@ -580,17 +576,17 @@
         //
         // From (1)
-        (spline->spline[i])->coeff[1] = -(yPtr->data.F32[i]) / H;
+        spline->spline[i]->coeff[1] = -(yPtr->data.F32[i]) / H;
         // From (2)
-        ((spline->spline[i])->coeff[1])+= (yPtr->data.F32[i+1] / H);
+        spline->spline[i]->coeff[1]+= yPtr->data.F32[i+1] / H;
         // From (3)
-        tmp = -3.0 * (xPtr->data.F32[i+1] * xPtr->data.F32[i+1]) / (H * H * H);
+        tmp = -3.0 * xPtr->data.F32[i+1] * xPtr->data.F32[i+1] / (H * H * H);
         tmp+= (1.0 / H);
-        tmp*= ((spline->p_psDeriv2)[i]) * H * H / 6.0;
-        ((spline->spline[i])->coeff[1])+= tmp;
+        tmp*= spline->p_psDeriv2[i] * H * H / 6.0;
+        spline->spline[i]->coeff[1]+= tmp;
         // From (4)
-        tmp = 3.0 * (xPtr->data.F32[i] * xPtr->data.F32[i]) / (H * H * H);
-        tmp-= (1.0 / H);
-        tmp*= ((spline->p_psDeriv2)[i+1]) * H * H / 6.0;
-        ((spline->spline[i])->coeff[1])+= tmp;
+        tmp = 3.0 * xPtr->data.F32[i] * xPtr->data.F32[i] / (H * H * H);
+        tmp-= 1.0 / H;
+        tmp*= spline->p_psDeriv2[i+1] * H * H / 6.0;
+        spline->spline[i]->coeff[1]+= tmp;
 
         //
@@ -598,7 +594,7 @@
         //
         // From (3)
-        (spline->spline[i])->coeff[2] = ((spline->p_psDeriv2)[i]) * 3.0 * xPtr->data.F32[i+1] / (6.0 * H);
+        spline->spline[i]->coeff[2] = ((spline->p_psDeriv2)[i]) * 3.0 * xPtr->data.F32[i+1] / (6.0 * H);
         // From (4)
-        ((spline->spline[i])->coeff[2])-= (((spline->p_psDeriv2)[i+1]) * 3.0 * xPtr->data.F32[i] / (6.0 * H));
+        (spline->spline[i]->coeff[2])-= (((spline->p_psDeriv2)[i+1]) * 3.0 * xPtr->data.F32[i] / (6.0 * H));
 
         //
@@ -606,19 +602,19 @@
         //
         // From (3)
-        (spline->spline[i])->coeff[3] = -((spline->p_psDeriv2)[i]) / (6.0 * H);
+        spline->spline[i]->coeff[3] = -((spline->p_psDeriv2)[i]) / (6.0 * H);
         // From (4)
-        ((spline->spline[i])->coeff[3])+=  ((spline->p_psDeriv2)[i+1]) / (6.0 * H);
-
-        psTrace(__func__, 6, "(spline->spline[%u])->coeff[0] is %f\n", i, (spline->spline[i])->coeff[0]);
-        psTrace(__func__, 6, "(spline->spline[%u])->coeff[1] is %f\n", i, (spline->spline[i])->coeff[1]);
-        psTrace(__func__, 6, "(spline->spline[%u])->coeff[2] is %f\n", i, (spline->spline[i])->coeff[2]);
-        psTrace(__func__, 6, "(spline->spline[%u])->coeff[3] is %f\n", i, (spline->spline[i])->coeff[3]);
-    }
-
-    if (PS_TYPE_F64 == x->type.type) {
-        psFree(xF32);
+        (spline->spline[i]->coeff[3])+=  ((spline->p_psDeriv2)[i+1]) / (6.0 * H);
+
+        psTrace(__func__, 6, "(spline->spline[%u])->coeff[0] is %f\n", i, spline->spline[i]->coeff[0]);
+        psTrace(__func__, 6, "(spline->spline[%u])->coeff[1] is %f\n", i, spline->spline[i]->coeff[1]);
+        psTrace(__func__, 6, "(spline->spline[%u])->coeff[2] is %f\n", i, spline->spline[i]->coeff[2]);
+        psTrace(__func__, 6, "(spline->spline[%u])->coeff[3] is %f\n", i, spline->spline[i]->coeff[3]);
+    }
+
+    if ((x == NULL) || (PS_TYPE_F64 == x->type.type)) {
+        psFree(xPtr);
     }
     if (PS_TYPE_F64 == y->type.type) {
-        psFree(yF32);
+        psFree(yPtr);
     }
     psTrace(__func__, 3, "---- %s() end ----\n", __func__);
@@ -637,5 +633,5 @@
      the spline fit functions require F32 and F64.
  
-XXX: This only works if spline0>knots if psF32.  Must we add support for psU32 and
+XXX: This only works if spline->knots if psF32.  Must we add support for psU32 and
 psF64?
  *****************************************************************************/