Index: trunk/psLib/src/math/psSpline.c
===================================================================
--- trunk/psLib/src/math/psSpline.c	(revision 1823)
+++ trunk/psLib/src/math/psSpline.c	(revision 1831)
@@ -7,6 +7,6 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-17 02:14:52 $
+ *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-18 01:50:45 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -1809,10 +1809,10 @@
 p_psInterpolate1D(): 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-order LaGrange interpolated
+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?
  *****************************************************************************/
-float p_psFullInterpolate1DF32(float *domain,
+float p_ps1DFullInterpolateF32(float *domain,
                                float *range,
                                int n,
@@ -1822,21 +1822,39 @@
     int m;
     static psVector *p = NULL;
-    psVectorRecycle(p, n, PS_TYPE_F32);
+    p = psVectorRecycle(p, n, PS_TYPE_F32);
     p_psMemSetPersistent(p, true);
+    p_psMemSetPersistent(p->data.F32, true);
+
+    psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
+            "---- p_ps1DFullInterpolateF32() begin (%d-order at x=%f) (%d data points)----\n", n-1, x, n);
+
+    for (i=0;i<n;i++) {
+        psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
+                "domain/range is (%f %f)\n", domain[i], range[i]);
+    }
+
+    for (i=0;i<n;i++) {
+        p->data.F32[i] = range[i];
+        psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
+                "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
+
+    }
 
     // From NR, during each iteration of the m loop, we are computing the
     // p_{i ... i+m} terms.
-    for (m=0;m<n;m++) {
+    for (m=1;m<n;m++) {
         for (i=0;i<n-m;i++) {
-            if (m == 0) {
-                p->data.F32[i] = range[i];
-            } else {
-                // From NR: we are computing P_{i ... i+m}
-                p->data.F32[i] = (((x-range[i+m]) * p->data.F32[i]) +
-                                  ((range[i]-x) * p->data.F32[i+1])) /
-                                 (domain[i] - domain[i+m]);
-            }
-        }
-    }
+            // From NR: we are computing P_{i ... i+m}
+            p->data.F32[i] = (((x-domain[i+m]) * p->data.F32[i]) +
+                              ((domain[i]-x) * p->data.F32[i+1])) /
+                             (domain[i] - domain[i+m]);
+            //printf("((%f-%f * %f) + (%f-%f * %f)) / (%f - %f)\n", x, domain[i+m], p->data.F32[i], domain[i], x, p->data.F32[i+1], domain[i], domain[i+m]);
+            psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 6,
+                    "p->data.F32[%d] is %f\n", i, p->data.F32[i]);
+        }
+    }
+    psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
+            "---- p_ps1DFullInterpolateF32() end ----\n");
+
     return(p->data.F32[0]);
 }
@@ -1844,5 +1862,5 @@
 
 // This is a
-float p_psInterpolate1DF32(float *domain,
+float p_ps1DInterpolateF32(float *domain,
                            float *range,
                            int n,
@@ -1854,5 +1872,9 @@
     int origin;
 
+    psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4,
+            "---- p_ps1DInterpolateF32() begin ----\n");
+
     binNum = VectorBinDisectF32(domain, n, x);
+
     if (0 == numIntPoints%2) {
         origin = binNum - ((numIntPoints/2) - 1);
@@ -1871,5 +1893,7 @@
     }
 
-    return(p_psFullInterpolate1DF32(&domain[origin], &range[origin], n, x));
+    psTrace(".psLib.dataManip.psFunctions.p_ps1DInterpolateF32", 4,
+            "---- p_ps1DInterpolateF32() end ----\n");
+    return(p_ps1DFullInterpolateF32(&domain[origin], &range[origin], order+1, x));
 }
 
@@ -1882,9 +1906,12 @@
 XXX: This stuff does not work with a mask.
  *****************************************************************************/
-float p_psInterpolate1D(psVector *domain,
-                        psVector *range,
-                        int order,
-                        psScalar *x)
-{
+psScalar *p_psVectorInterpolate(psVector *domain,
+                                psVector *range,
+                                int order,
+                                psScalar *x)
+{
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
+            "---- p_psVectorInterpolate() begin ----\n");
+
     if (domain->type.type != range->type.type != x->type.type) {
         // XXX psError
@@ -1898,11 +1925,18 @@
 
     if (x->type.type == PS_TYPE_F32) {
-        return(p_psInterpolate1DF32(domain->data.F32, range->data.F32,
-                                    domain->n, order, x->data.F32));
+        psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
+                "---- p_psVectorInterpolate() end ----\n");
+        return(psScalarAlloc(p_ps1DInterpolateF32(domain->data.F32,
+                             range->data.F32,
+                             domain->n,
+                             order,
+                             x->data.F32), PS_TYPE_F32));
     } else {
         // XXX psError: type not supported
     }
 
-    return(-1.0);
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4,
+            "---- p_psVectorInterpolate() end ----\n");
+    return(NULL);
 }