Changeset 1859

Timestamp:

Sep 22, 2004, 6:56:41 PM (22 years ago)

Author:

gusciora

Message:

Implemented the generation spline functions.

Location:

trunk/psLib/src

Files:

• dataManip/psFunctions.c (modified) (13 diffs)
• dataManip/psFunctions.h (modified) (3 diffs)
• dataManip/psMinimize.c (modified) (22 diffs)
• math/psMinimize.c (modified) (22 diffs)
• math/psPolynomial.c (modified) (13 diffs)
• math/psPolynomial.h (modified) (3 diffs)
• math/psSpline.c (modified) (13 diffs)
• math/psSpline.h (modified) (3 diffs)

trunk/psLib/src/dataManip/psFunctions.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-22 01:30:21 $
+ *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-23 04:56:40 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -467 +467 @@
     for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
         if (myPoly->mask[loop_x] == 0) {
+            psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 10,
+                    "polysum+= sum*coeff [%f+= (%f * %f)\n", polySum, xSum, myPoly->coeff[loop_x]);
             polySum += xSum * myPoly->coeff[loop_x];
             xSum *= x;
@@ -1687 +1689 @@
 
 /*****************************************************************************
-VectorBinDisectF32(): This is a private function which takes as input a
+p_psVectorBinDisectF32(): This is a private function which takes as input a
 vector of floating point data as well as a single floating point values.
 The input vector values are assumed to be non-decreasing (v[i-1] <= v[j] for
@@ -1696 +1698 @@
 XXX: Macro this for a few different types.
  *****************************************************************************/
-int VectorBinDisectF32(float *bins,
-                       int numBins,
-                       float x)
+int p_psVectorBinDisectF32(float *bins,
+                           int numBins,
+                           float x)
 {
     int min;
@@ -1704 +1706 @@
     int mid;
 
-    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
-            "---- Calling VectorBinDisectF32(%f)\n", x);
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4,
+            "---- Calling p_psVectorBinDisectF32(%f)\n", x);
 
     if ((x < bins[0]) ||
             (x > bins[numBins-1])) {
         psLogMsg(__func__, PS_LOG_WARN,
-                 "VectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).",
+                 "p_psVectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).",
                  x, bins[0], bins[numBins-1]);
         return(-1);
@@ -1720 +1722 @@
 
     while (min != max) {
-        psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
+        psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4,
                 "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
                 min, mid, max, x, bins[mid]);
 
         if (x == bins[mid]) {
-            psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
-                    "---- Exiting VectorBinDisectF32(): bin %d\n", mid);
+            psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4,
+                    "---- Exiting p_psVectorBinDisectF32(): bin %d\n", mid);
             return(mid);
         } else if (x < bins[mid]) {
@@ -1736 +1738 @@
     }
 
-    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
-            "---- Exiting VectorBinDisectF32(): bin %d\n", min);
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4,
+            "---- Exiting p_psVectorBinDisectF32(): bin %d\n", min);
     return(min);
 }
 
 /*****************************************************************************
-VectorBinDisectS32(): integer version of above.
+p_psVectorBinDisectS32(): integer version of above.
  *****************************************************************************/
-int VectorBinDisectS32(int *bins,
-                       int numBins,
-                       int x)
+int p_psVectorBinDisectS32(int *bins,
+                           int numBins,
+                           int x)
 {
     int min;
@@ -1752 +1754 @@
     int mid;
 
-    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
-            "---- Calling VectorBinDisectS32(%f)\n", x);
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4,
+            "---- Calling p_psVectorBinDisectS32(%f)\n", x);
 
     if ((x < bins[0]) ||
             (x > bins[numBins-1])) {
         psLogMsg(__func__, PS_LOG_WARN,
-                 "VectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).",
+                 "p_psVectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).",
                  x, bins[0], bins[numBins-1]);
         return(-1);
@@ -1768 +1770 @@
 
     while (min != max) {
-        psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
+        psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4,
                 "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
                 min, mid, max, x, bins[mid]);
 
         if (x == bins[mid]) {
-            psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
-                    "---- Exiting VectorBinDisectS32(): bin %d\n", min);
+            psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4,
+                    "---- Exiting p_psVectorBinDisectS32(): bin %d\n", min);
             return(min);
         } else if (x < bins[mid]) {
@@ -1784 +1786 @@
     }
 
-    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
-            "---- Exiting VectorBinDisectS32(): bin %d\n", min);
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4,
+            "---- Exiting p_psVectorBinDisectS32(): bin %d\n", min);
     return(min);
 }
 
 /*****************************************************************************
-p_psVectorBinDisect(): A wrapper to the above VectorBinDisect().
+p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect().
  *****************************************************************************/
 int p_psVectorBinDisect(psVector *bins,
@@ -1801 +1803 @@
 
     if (x->type.type == PS_TYPE_S32) {
-        return(VectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));
+        return(p_psVectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));
     } else if (x->type.type == PS_TYPE_F32) {
-        return(VectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));
+        return(p_psVectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));
     } else {
         // XXX: Generate error message.
@@ -1879 +1881 @@
             "---- p_ps1DInterpolateF32() begin ----\n");
 
-    binNum = VectorBinDisectF32(domain, n, x);
+    binNum = p_psVectorBinDisectF32(domain, n, x);
 
     if (0 == numIntPoints%2) {
@@ -1953 +1955 @@
 
     n = spline->n;
-    binNum = VectorBinDisectF32(spline->domains, (spline->n)+1, x);
+    binNum = p_psVectorBinDisectF32(spline->domains, (spline->n)+1, x);
     if (binNum == -1) {
         psLogMsg(__func__, PS_LOG_WARN,

trunk/psLib/src/dataManip/psFunctions.h

@@ -12 +12 @@
 *  @author George Gusciora, MHPCC
 *
-*  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2004-09-22 01:30:21 $
+*  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2004-09-23 04:56:41 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -407 +407 @@
                           psVector *data);
 
+int p_psVectorBinDisectF32(float *bins,
+                           int numBins,
+                           float x);
+
+int p_psVectorBinDisectS32(int *bins,
+                           int numBins,
+                           int x);
+
 int p_psVectorBinDisect(psVector *bins,
                         psScalar *x);
@@ -415 +423 @@
                                 psScalar *x);
 
+float p_psNRSpline1DEval(psSpline1D *spline,
+                         const psVector* restrict x,
+                         const psVector* restrict y,
+                         float X);
+
 /* \} */// End of MathGroup Functions
 

trunk/psLib/src/dataManip/psMinimize.c

@@ -9 +9 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-22 02:42:49 $
+ *  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-23 04:56:41 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -43 +43 @@
 #include "psMinimize.h"
 #include "psMatrix.h"
-#include "psTrace.h"
+#include "psConstants.h"
 
 /*****************************************************************************/
 /* DEFINE STATEMENTS                                                         */
 /*****************************************************************************/
-
-#define MAX_LMM_ITERATIONS 100
-#define MAX_MINIMIZE_ITERATIONS 100
 
 /** Preprocessor macro to generate error on a NULL 1DPolynomial */
@@ -96 +93 @@
 /* GLOBAL VARIABLES                                                          */
 /*****************************************************************************/
-
 static psMinimizeChi2PowellFunc Chi2PowellFunc = NULL;
 //static psMinimizePowellFunc PowellFunc = NULL;
@@ -116 +112 @@
 input parameter "x" between 0 and input parameter polyOrder.  The result is
 returned as a psVector sums.
+ 
+XXX: change name
  *****************************************************************************/
-void p_psBuildSums1D(double x, int polyOrder, psVector* sums)
+void p_psBuildSums1D(double x,
+                     int polyOrder,
+                     psVector* sums)
 {
     int i = 0;
@@ -148 +148 @@
 
 /*****************************************************************************
- 
 CalculateSecondDerivs(): Given a set of x/y vectors corresponding to a
-tabulated function at n points, this routine calcualtes the second
+tabulated function at n points, this routine calculates the second
 derivatives of the interpolating cubic splines at those n points.
  
 The first and second derivatives at the endpoints, undefined in the SDR, are
-here defined to be 0.0.
- 
-XXX: This algorithm is very similar to that in Numerical Recipes.
+here defined to be 0.0.  They can be modified via ypo and yp1.
+ 
+XXX: This algorithm is derived from the Numerical Recipes.
  *****************************************************************************/
 float *CalculateSecondDerivs(const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
@@ -174 +173 @@
     float *X = (float *) & (x->data.F32[0]);
     float *Y = (float *) & (y->data.F32[0]);
+    float qn;
 
     if (x == NULL) {
@@ -183 +183 @@
     }
 
-    // XXX: The first and second derivatives at the endpoints, undefined in
-    // the SDR, are here defined to be 0.0.
-    u[0]= 0.0;
-    u[n-1]= 0.0;
-    derivs2[0] = 0.0;
-    derivs2[n-1] = 0.0;
+    // XXX: The second derivatives at the endpoints, undefined in the SDR,
+    // are set in psConstants.h: LEFT_SPLINE_DERIV, RIGHT_SPLINE_DERIV.
+    derivs2[0] = -0.5;
+    u[0]= (3.0/(X[1]-X[0])) * ((Y[1]-Y[0])/(X[1]-X[0]) - LEFT_SPLINE_DERIV);
 
     for (i=1;i<=(n-2);i++) {
@@ -194 +192 @@
         p = sig * derivs2[i-1] + 2.0;
         derivs2[i] = (sig - 1.0) / p;
-        u[i] = (Y[i+1] - Y[i])/(X[i+1]-X[i]) - (Y[i]-Y[i-1])/(X[i]-X[i-1]);
-        u[i] = (6.0 * u[i] / (X[i+1] - X[i-1]) - sig * u[i-1]) / p;
+        u[i] = ((Y[i+1] - Y[i])/(X[i+1]-X[i])) - ((Y[i]-Y[i-1])/(X[i]-X[i-1]));
+        u[i] = ((6.0 * u[i] / (X[i+1] - X[i-1])) - (sig * u[i-1])) / p;
 
         psTrace(".psLib.dataManip.CalculateSecondDerivs", 6,
@@ -205 +203 @@
     }
 
+    qn = 0.5;
+    u[n-1] = (3.0/(X[n-1]-X[n-2])) * (RIGHT_SPLINE_DERIV - (Y[n-1]-Y[n-2])/(X[n-1]-X[n-2]));
+    derivs2[n-1] = (u[n-1] - (qn * u[n-2])) / ((qn * derivs2[n-2]) + 1.0);
+
     for (k=(n-2);k>=0;k--) {
         derivs2[k] = derivs2[k] * derivs2[k+1] + u[k];
+
         psTrace(".psLib.dataManip.CalculateSecondDerivs", 6,
                 "derivs2[%d] is %f\n", k, derivs2[k]);
     }
+
     if (mustFreeX == true) {
         psFree(X);
@@ -224 +228 @@
 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
 /*****************************************************************************/
-
-
 
 /*****************************************************************************
@@ -238 +240 @@
  (3)         ((((A*A*A)-A) * mySpline->p_psDeriv2[0]) * H^2)/6.0 +
  (4)         ((((B*B*B)-B) * mySpline->p_psDeriv2[1]) * H^2)/6.0
-Where
- H = x[1]-x[0], A = (x[1]-x)/H, B = (x-x[0])/H
+Where:
+ H = x[1]-x[0]
+ A = (x[1]-x)/H
+ B = (x-x[0])/H
 The bulk of the code in this routine is the expansion of the above equation
-into a polynomial in terms of x.  This gets pretty complicated.
+into a polynomial in terms of x, and then saving the coefficients of the
+powers of x in the spline polynomials.  This gets pretty complicated.
  
 XXX: usage of yErr is not specified in IfA documentation.
@@ -262 +267 @@
     }
     mySpline->p_psDeriv2 = CalculateSecondDerivs(x, y);
+    for (i=0;i<y->n;i++)
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
 
     for (i=0;i<numSplines;i++) {
@@ -271 +279 @@
         //
         // From (1)
-        (mySpline->spline[i])->coeff[0] = y->data.F32[i] * (x->data.F32[i+1]/H);
+        (mySpline->spline[i])->coeff[0] = (y->data.F32[i] * x->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])-= ((y->data.F32[i+1] * x->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);
@@ -289 +297 @@
         //
         // From (1)
-        (mySpline->spline[i])->coeff[1] = - ((mySpline->p_psDeriv2)[i]) / H;
+        (mySpline->spline[i])->coeff[1] = -(y->data.F32[i]) / H;
         // From (2)
-        (mySpline->spline[i])->coeff[1]-= ((mySpline->p_psDeriv2)[i+1]) / H;
+        ((mySpline->spline[i])->coeff[1])+= (y->data.F32[i+1] / H);
         // From (3)
-        tmp = - (x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H);
+        tmp = -3.0 * (x->data.F32[i+1] * x->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;
+        ((mySpline->spline[i])->coeff[1])+= tmp;
         // From (4)
-        tmp = (x->data.F32[i] + x->data.F32[i]) / (H * H * H);
+        tmp = 3.0 * (x->data.F32[i] * x->data.F32[i]) / (H * H * H);
         tmp-= (1.0 / H);
         tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
-        (mySpline->spline[i])->coeff[1]+= tmp;
+        ((mySpline->spline[i])->coeff[1])+= tmp;
 
         //
@@ -307 +315 @@
         //
         // From (3)
-        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x->data.F32[i+1] / 6.0;
+        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x->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;
+        ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x->data.F32[i] / (6.0 * H));
 
         //
@@ -317 +325 @@
         (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
         // From (4)
-        (mySpline->spline[i])->coeff[3]+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
+        ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
 
         psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
@@ -334 +342 @@
     return(mySpline);
 }
+
+/******************************************************************************
+p_psNRSpline1DEval(): This routine does NR-style evaluation of cubic splines.
+It takes advantage of the 2nd derivatives of the cubic splines, which are
+stored in the psSPline1D data structure, and computes the interpolated value
+directly, without computing (or using) the interpolating cubic spline
+polynomial.
+ 
+This routine is here mostly for a sanity check on the psLib function
+evalSpline() which computes the interpolated value based on the cubic spline
+polynomials which are stored in psSpline1D.
+ *****************************************************************************/
+float p_psNRSpline1DEval(psSpline1D *spline,
+                         const psVector* restrict x,
+                         const psVector* restrict y,
+                         float X)
+{
+    int n;
+    int klo;
+    int khi;
+    float H;
+    float A;
+    float B;
+    float C;
+    float D;
+    float Y;
+
+    n = spline->n;
+    klo = p_psVectorBinDisectF32(spline->domains, (spline->n)+1, X);
+    khi = klo + 1;
+    H = (spline->domains)[khi] - (spline->domains)[klo];
+    A = ((spline->domains)[khi] - X) / H;
+    B = (X - (spline->domains)[klo]) / H;
+    C = ((A*A*A)-A) * (H*H/6.0);
+    D = ((B*B*B)-B) * (H*H/6.0);
+
+    Y = (A * y->data.F32[klo]) +
+        (B * y->data.F32[khi]) +
+        (C * (spline->p_psDeriv2)[klo]) +
+        (D * (spline->p_psDeriv2)[khi]);
+
+    return(Y);
+}
+
 
 /******************************************************************************
@@ -687 +739 @@
 }
 
-
-#define STEP_SIZE 0.10
 /******************************************************************************
-    This routine takes as input an arbitrary function, and the parameter to
-    vary, and the line along which it must vary.  This function produces as
-    output a bracket [a, b, c] such that
+p_psDetermineBracket():  This routine takes as input an arbitrary function,
+and the parameter to vary, and the line along which it must vary.  This
+function produces as output a bracket [a, b, c] such that
+ 
         f(param + b * line) is less than f(param + a * line) and
                                          f(param + c * line).
@@ -700 +751 @@
     smaller/larger than b.  Repeat this process until a local minimum is
     found.
+ 
+ 
  *****************************************************************************/
 psVector *p_psDetermineBracket(psVector *params,
@@ -719 +772 @@
     float new_cDir = 0.0;
     psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
-    float stepSize = STEP_SIZE;
+    float stepSize = DETERMINE_BRACKET_STEP_SIZE;
     psVector *tmp = NULL;
     int i = 0;
@@ -869 +922 @@
     return(NULL);
 }
-
-
-
-
-
-
-
-
-
-
-
 
 /******************************************************************************
@@ -1061 +1103 @@
 
 
-
-
-
-
-
-
-
-
-
-
-
-
-
 /******************************************************************************
 This routine must minimize a possibly multi-dimensional function.  The

trunk/psLib/src/math/psMinimize.c

@@ -9 +9 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-22 02:42:49 $
+ *  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-23 04:56:41 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -43 +43 @@
 #include "psMinimize.h"
 #include "psMatrix.h"
-#include "psTrace.h"
+#include "psConstants.h"
 
 /*****************************************************************************/
 /* DEFINE STATEMENTS                                                         */
 /*****************************************************************************/
-
-#define MAX_LMM_ITERATIONS 100
-#define MAX_MINIMIZE_ITERATIONS 100
 
 /** Preprocessor macro to generate error on a NULL 1DPolynomial */
@@ -96 +93 @@
 /* GLOBAL VARIABLES                                                          */
 /*****************************************************************************/
-
 static psMinimizeChi2PowellFunc Chi2PowellFunc = NULL;
 //static psMinimizePowellFunc PowellFunc = NULL;
@@ -116 +112 @@
 input parameter "x" between 0 and input parameter polyOrder.  The result is
 returned as a psVector sums.
+ 
+XXX: change name
  *****************************************************************************/
-void p_psBuildSums1D(double x, int polyOrder, psVector* sums)
+void p_psBuildSums1D(double x,
+                     int polyOrder,
+                     psVector* sums)
 {
     int i = 0;
@@ -148 +148 @@
 
 /*****************************************************************************
- 
 CalculateSecondDerivs(): Given a set of x/y vectors corresponding to a
-tabulated function at n points, this routine calcualtes the second
+tabulated function at n points, this routine calculates the second
 derivatives of the interpolating cubic splines at those n points.
  
 The first and second derivatives at the endpoints, undefined in the SDR, are
-here defined to be 0.0.
- 
-XXX: This algorithm is very similar to that in Numerical Recipes.
+here defined to be 0.0.  They can be modified via ypo and yp1.
+ 
+XXX: This algorithm is derived from the Numerical Recipes.
  *****************************************************************************/
 float *CalculateSecondDerivs(const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
@@ -174 +173 @@
     float *X = (float *) & (x->data.F32[0]);
     float *Y = (float *) & (y->data.F32[0]);
+    float qn;
 
     if (x == NULL) {
@@ -183 +183 @@
     }
 
-    // XXX: The first and second derivatives at the endpoints, undefined in
-    // the SDR, are here defined to be 0.0.
-    u[0]= 0.0;
-    u[n-1]= 0.0;
-    derivs2[0] = 0.0;
-    derivs2[n-1] = 0.0;
+    // XXX: The second derivatives at the endpoints, undefined in the SDR,
+    // are set in psConstants.h: LEFT_SPLINE_DERIV, RIGHT_SPLINE_DERIV.
+    derivs2[0] = -0.5;
+    u[0]= (3.0/(X[1]-X[0])) * ((Y[1]-Y[0])/(X[1]-X[0]) - LEFT_SPLINE_DERIV);
 
     for (i=1;i<=(n-2);i++) {
@@ -194 +192 @@
         p = sig * derivs2[i-1] + 2.0;
         derivs2[i] = (sig - 1.0) / p;
-        u[i] = (Y[i+1] - Y[i])/(X[i+1]-X[i]) - (Y[i]-Y[i-1])/(X[i]-X[i-1]);
-        u[i] = (6.0 * u[i] / (X[i+1] - X[i-1]) - sig * u[i-1]) / p;
+        u[i] = ((Y[i+1] - Y[i])/(X[i+1]-X[i])) - ((Y[i]-Y[i-1])/(X[i]-X[i-1]));
+        u[i] = ((6.0 * u[i] / (X[i+1] - X[i-1])) - (sig * u[i-1])) / p;
 
         psTrace(".psLib.dataManip.CalculateSecondDerivs", 6,
@@ -205 +203 @@
     }
 
+    qn = 0.5;
+    u[n-1] = (3.0/(X[n-1]-X[n-2])) * (RIGHT_SPLINE_DERIV - (Y[n-1]-Y[n-2])/(X[n-1]-X[n-2]));
+    derivs2[n-1] = (u[n-1] - (qn * u[n-2])) / ((qn * derivs2[n-2]) + 1.0);
+
     for (k=(n-2);k>=0;k--) {
         derivs2[k] = derivs2[k] * derivs2[k+1] + u[k];
+
         psTrace(".psLib.dataManip.CalculateSecondDerivs", 6,
                 "derivs2[%d] is %f\n", k, derivs2[k]);
     }
+
     if (mustFreeX == true) {
         psFree(X);
@@ -224 +228 @@
 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
 /*****************************************************************************/
-
-
 
 /*****************************************************************************
@@ -238 +240 @@
  (3)         ((((A*A*A)-A) * mySpline->p_psDeriv2[0]) * H^2)/6.0 +
  (4)         ((((B*B*B)-B) * mySpline->p_psDeriv2[1]) * H^2)/6.0
-Where
- H = x[1]-x[0], A = (x[1]-x)/H, B = (x-x[0])/H
+Where:
+ H = x[1]-x[0]
+ A = (x[1]-x)/H
+ B = (x-x[0])/H
 The bulk of the code in this routine is the expansion of the above equation
-into a polynomial in terms of x.  This gets pretty complicated.
+into a polynomial in terms of x, and then saving the coefficients of the
+powers of x in the spline polynomials.  This gets pretty complicated.
  
 XXX: usage of yErr is not specified in IfA documentation.
@@ -262 +267 @@
     }
     mySpline->p_psDeriv2 = CalculateSecondDerivs(x, y);
+    for (i=0;i<y->n;i++)
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
 
     for (i=0;i<numSplines;i++) {
@@ -271 +279 @@
         //
         // From (1)
-        (mySpline->spline[i])->coeff[0] = y->data.F32[i] * (x->data.F32[i+1]/H);
+        (mySpline->spline[i])->coeff[0] = (y->data.F32[i] * x->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])-= ((y->data.F32[i+1] * x->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);
@@ -289 +297 @@
         //
         // From (1)
-        (mySpline->spline[i])->coeff[1] = - ((mySpline->p_psDeriv2)[i]) / H;
+        (mySpline->spline[i])->coeff[1] = -(y->data.F32[i]) / H;
         // From (2)
-        (mySpline->spline[i])->coeff[1]-= ((mySpline->p_psDeriv2)[i+1]) / H;
+        ((mySpline->spline[i])->coeff[1])+= (y->data.F32[i+1] / H);
         // From (3)
-        tmp = - (x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H);
+        tmp = -3.0 * (x->data.F32[i+1] * x->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;
+        ((mySpline->spline[i])->coeff[1])+= tmp;
         // From (4)
-        tmp = (x->data.F32[i] + x->data.F32[i]) / (H * H * H);
+        tmp = 3.0 * (x->data.F32[i] * x->data.F32[i]) / (H * H * H);
         tmp-= (1.0 / H);
         tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
-        (mySpline->spline[i])->coeff[1]+= tmp;
+        ((mySpline->spline[i])->coeff[1])+= tmp;
 
         //
@@ -307 +315 @@
         //
         // From (3)
-        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x->data.F32[i+1] / 6.0;
+        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x->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;
+        ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x->data.F32[i] / (6.0 * H));
 
         //
@@ -317 +325 @@
         (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
         // From (4)
-        (mySpline->spline[i])->coeff[3]+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
+        ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
 
         psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
@@ -334 +342 @@
     return(mySpline);
 }
+
+/******************************************************************************
+p_psNRSpline1DEval(): This routine does NR-style evaluation of cubic splines.
+It takes advantage of the 2nd derivatives of the cubic splines, which are
+stored in the psSPline1D data structure, and computes the interpolated value
+directly, without computing (or using) the interpolating cubic spline
+polynomial.
+ 
+This routine is here mostly for a sanity check on the psLib function
+evalSpline() which computes the interpolated value based on the cubic spline
+polynomials which are stored in psSpline1D.
+ *****************************************************************************/
+float p_psNRSpline1DEval(psSpline1D *spline,
+                         const psVector* restrict x,
+                         const psVector* restrict y,
+                         float X)
+{
+    int n;
+    int klo;
+    int khi;
+    float H;
+    float A;
+    float B;
+    float C;
+    float D;
+    float Y;
+
+    n = spline->n;
+    klo = p_psVectorBinDisectF32(spline->domains, (spline->n)+1, X);
+    khi = klo + 1;
+    H = (spline->domains)[khi] - (spline->domains)[klo];
+    A = ((spline->domains)[khi] - X) / H;
+    B = (X - (spline->domains)[klo]) / H;
+    C = ((A*A*A)-A) * (H*H/6.0);
+    D = ((B*B*B)-B) * (H*H/6.0);
+
+    Y = (A * y->data.F32[klo]) +
+        (B * y->data.F32[khi]) +
+        (C * (spline->p_psDeriv2)[klo]) +
+        (D * (spline->p_psDeriv2)[khi]);
+
+    return(Y);
+}
+
 
 /******************************************************************************
@@ -687 +739 @@
 }
 
-
-#define STEP_SIZE 0.10
 /******************************************************************************
-    This routine takes as input an arbitrary function, and the parameter to
-    vary, and the line along which it must vary.  This function produces as
-    output a bracket [a, b, c] such that
+p_psDetermineBracket():  This routine takes as input an arbitrary function,
+and the parameter to vary, and the line along which it must vary.  This
+function produces as output a bracket [a, b, c] such that
+ 
         f(param + b * line) is less than f(param + a * line) and
                                          f(param + c * line).
@@ -700 +751 @@
     smaller/larger than b.  Repeat this process until a local minimum is
     found.
+ 
+ 
  *****************************************************************************/
 psVector *p_psDetermineBracket(psVector *params,
@@ -719 +772 @@
     float new_cDir = 0.0;
     psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
-    float stepSize = STEP_SIZE;
+    float stepSize = DETERMINE_BRACKET_STEP_SIZE;
     psVector *tmp = NULL;
     int i = 0;
@@ -869 +922 @@
     return(NULL);
 }
-
-
-
-
-
-
-
-
-
-
-
 
 /******************************************************************************
@@ -1061 +1103 @@
 
 
-
-
-
-
-
-
-
-
-
-
-
-
-
 /******************************************************************************
 This routine must minimize a possibly multi-dimensional function.  The

trunk/psLib/src/math/psPolynomial.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-22 01:30:21 $
+ *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-23 04:56:40 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -467 +467 @@
     for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
         if (myPoly->mask[loop_x] == 0) {
+            psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 10,
+                    "polysum+= sum*coeff [%f+= (%f * %f)\n", polySum, xSum, myPoly->coeff[loop_x]);
             polySum += xSum * myPoly->coeff[loop_x];
             xSum *= x;
@@ -1687 +1689 @@
 
 /*****************************************************************************
-VectorBinDisectF32(): This is a private function which takes as input a
+p_psVectorBinDisectF32(): This is a private function which takes as input a
 vector of floating point data as well as a single floating point values.
 The input vector values are assumed to be non-decreasing (v[i-1] <= v[j] for
@@ -1696 +1698 @@
 XXX: Macro this for a few different types.
  *****************************************************************************/
-int VectorBinDisectF32(float *bins,
-                       int numBins,
-                       float x)
+int p_psVectorBinDisectF32(float *bins,
+                           int numBins,
+                           float x)
 {
     int min;
@@ -1704 +1706 @@
     int mid;
 
-    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
-            "---- Calling VectorBinDisectF32(%f)\n", x);
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4,
+            "---- Calling p_psVectorBinDisectF32(%f)\n", x);
 
     if ((x < bins[0]) ||
             (x > bins[numBins-1])) {
         psLogMsg(__func__, PS_LOG_WARN,
-                 "VectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).",
+                 "p_psVectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).",
                  x, bins[0], bins[numBins-1]);
         return(-1);
@@ -1720 +1722 @@
 
     while (min != max) {
-        psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
+        psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4,
                 "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
                 min, mid, max, x, bins[mid]);
 
         if (x == bins[mid]) {
-            psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
-                    "---- Exiting VectorBinDisectF32(): bin %d\n", mid);
+            psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4,
+                    "---- Exiting p_psVectorBinDisectF32(): bin %d\n", mid);
             return(mid);
         } else if (x < bins[mid]) {
@@ -1736 +1738 @@
     }
 
-    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
-            "---- Exiting VectorBinDisectF32(): bin %d\n", min);
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4,
+            "---- Exiting p_psVectorBinDisectF32(): bin %d\n", min);
     return(min);
 }
 
 /*****************************************************************************
-VectorBinDisectS32(): integer version of above.
+p_psVectorBinDisectS32(): integer version of above.
  *****************************************************************************/
-int VectorBinDisectS32(int *bins,
-                       int numBins,
-                       int x)
+int p_psVectorBinDisectS32(int *bins,
+                           int numBins,
+                           int x)
 {
     int min;
@@ -1752 +1754 @@
     int mid;
 
-    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
-            "---- Calling VectorBinDisectS32(%f)\n", x);
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4,
+            "---- Calling p_psVectorBinDisectS32(%f)\n", x);
 
     if ((x < bins[0]) ||
             (x > bins[numBins-1])) {
         psLogMsg(__func__, PS_LOG_WARN,
-                 "VectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).",
+                 "p_psVectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).",
                  x, bins[0], bins[numBins-1]);
         return(-1);
@@ -1768 +1770 @@
 
     while (min != max) {
-        psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
+        psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4,
                 "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
                 min, mid, max, x, bins[mid]);
 
         if (x == bins[mid]) {
-            psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
-                    "---- Exiting VectorBinDisectS32(): bin %d\n", min);
+            psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4,
+                    "---- Exiting p_psVectorBinDisectS32(): bin %d\n", min);
             return(min);
         } else if (x < bins[mid]) {
@@ -1784 +1786 @@
     }
 
-    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
-            "---- Exiting VectorBinDisectS32(): bin %d\n", min);
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4,
+            "---- Exiting p_psVectorBinDisectS32(): bin %d\n", min);
     return(min);
 }
 
 /*****************************************************************************
-p_psVectorBinDisect(): A wrapper to the above VectorBinDisect().
+p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect().
  *****************************************************************************/
 int p_psVectorBinDisect(psVector *bins,
@@ -1801 +1803 @@
 
     if (x->type.type == PS_TYPE_S32) {
-        return(VectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));
+        return(p_psVectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));
     } else if (x->type.type == PS_TYPE_F32) {
-        return(VectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));
+        return(p_psVectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));
     } else {
         // XXX: Generate error message.
@@ -1879 +1881 @@
             "---- p_ps1DInterpolateF32() begin ----\n");
 
-    binNum = VectorBinDisectF32(domain, n, x);
+    binNum = p_psVectorBinDisectF32(domain, n, x);
 
     if (0 == numIntPoints%2) {
@@ -1953 +1955 @@
 
     n = spline->n;
-    binNum = VectorBinDisectF32(spline->domains, (spline->n)+1, x);
+    binNum = p_psVectorBinDisectF32(spline->domains, (spline->n)+1, x);
     if (binNum == -1) {
         psLogMsg(__func__, PS_LOG_WARN,

trunk/psLib/src/math/psPolynomial.h

@@ -12 +12 @@
 *  @author George Gusciora, MHPCC
 *
-*  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2004-09-22 01:30:21 $
+*  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2004-09-23 04:56:41 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -407 +407 @@
                           psVector *data);
 
+int p_psVectorBinDisectF32(float *bins,
+                           int numBins,
+                           float x);
+
+int p_psVectorBinDisectS32(int *bins,
+                           int numBins,
+                           int x);
+
 int p_psVectorBinDisect(psVector *bins,
                         psScalar *x);
@@ -415 +423 @@
                                 psScalar *x);
 
+float p_psNRSpline1DEval(psSpline1D *spline,
+                         const psVector* restrict x,
+                         const psVector* restrict y,
+                         float X);
+
 /* \} */// End of MathGroup Functions
 

trunk/psLib/src/math/psSpline.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-22 01:30:21 $
+ *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-23 04:56:40 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -467 +467 @@
     for (loop_x = 0; loop_x < myPoly->n; loop_x++) {
         if (myPoly->mask[loop_x] == 0) {
+            psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 10,
+                    "polysum+= sum*coeff [%f+= (%f * %f)\n", polySum, xSum, myPoly->coeff[loop_x]);
             polySum += xSum * myPoly->coeff[loop_x];
             xSum *= x;
@@ -1687 +1689 @@
 
 /*****************************************************************************
-VectorBinDisectF32(): This is a private function which takes as input a
+p_psVectorBinDisectF32(): This is a private function which takes as input a
 vector of floating point data as well as a single floating point values.
 The input vector values are assumed to be non-decreasing (v[i-1] <= v[j] for
@@ -1696 +1698 @@
 XXX: Macro this for a few different types.
  *****************************************************************************/
-int VectorBinDisectF32(float *bins,
-                       int numBins,
-                       float x)
+int p_psVectorBinDisectF32(float *bins,
+                           int numBins,
+                           float x)
 {
     int min;
@@ -1704 +1706 @@
     int mid;
 
-    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
-            "---- Calling VectorBinDisectF32(%f)\n", x);
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4,
+            "---- Calling p_psVectorBinDisectF32(%f)\n", x);
 
     if ((x < bins[0]) ||
             (x > bins[numBins-1])) {
         psLogMsg(__func__, PS_LOG_WARN,
-                 "VectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).",
+                 "p_psVectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).",
                  x, bins[0], bins[numBins-1]);
         return(-1);
@@ -1720 +1722 @@
 
     while (min != max) {
-        psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
+        psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4,
                 "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
                 min, mid, max, x, bins[mid]);
 
         if (x == bins[mid]) {
-            psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
-                    "---- Exiting VectorBinDisectF32(): bin %d\n", mid);
+            psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4,
+                    "---- Exiting p_psVectorBinDisectF32(): bin %d\n", mid);
             return(mid);
         } else if (x < bins[mid]) {
@@ -1736 +1738 @@
     }
 
-    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4,
-            "---- Exiting VectorBinDisectF32(): bin %d\n", min);
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4,
+            "---- Exiting p_psVectorBinDisectF32(): bin %d\n", min);
     return(min);
 }
 
 /*****************************************************************************
-VectorBinDisectS32(): integer version of above.
+p_psVectorBinDisectS32(): integer version of above.
  *****************************************************************************/
-int VectorBinDisectS32(int *bins,
-                       int numBins,
-                       int x)
+int p_psVectorBinDisectS32(int *bins,
+                           int numBins,
+                           int x)
 {
     int min;
@@ -1752 +1754 @@
     int mid;
 
-    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
-            "---- Calling VectorBinDisectS32(%f)\n", x);
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4,
+            "---- Calling p_psVectorBinDisectS32(%f)\n", x);
 
     if ((x < bins[0]) ||
             (x > bins[numBins-1])) {
         psLogMsg(__func__, PS_LOG_WARN,
-                 "VectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).",
+                 "p_psVectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).",
                  x, bins[0], bins[numBins-1]);
         return(-1);
@@ -1768 +1770 @@
 
     while (min != max) {
-        psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
+        psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4,
                 "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n",
                 min, mid, max, x, bins[mid]);
 
         if (x == bins[mid]) {
-            psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
-                    "---- Exiting VectorBinDisectS32(): bin %d\n", min);
+            psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4,
+                    "---- Exiting p_psVectorBinDisectS32(): bin %d\n", min);
             return(min);
         } else if (x < bins[mid]) {
@@ -1784 +1786 @@
     }
 
-    psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4,
-            "---- Exiting VectorBinDisectS32(): bin %d\n", min);
+    psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4,
+            "---- Exiting p_psVectorBinDisectS32(): bin %d\n", min);
     return(min);
 }
 
 /*****************************************************************************
-p_psVectorBinDisect(): A wrapper to the above VectorBinDisect().
+p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect().
  *****************************************************************************/
 int p_psVectorBinDisect(psVector *bins,
@@ -1801 +1803 @@
 
     if (x->type.type == PS_TYPE_S32) {
-        return(VectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));
+        return(p_psVectorBinDisectS32(bins->data.S32, bins->n, x->data.S32));
     } else if (x->type.type == PS_TYPE_F32) {
-        return(VectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));
+        return(p_psVectorBinDisectF32(bins->data.F32, bins->n, x->data.F32));
     } else {
         // XXX: Generate error message.
@@ -1879 +1881 @@
             "---- p_ps1DInterpolateF32() begin ----\n");
 
-    binNum = VectorBinDisectF32(domain, n, x);
+    binNum = p_psVectorBinDisectF32(domain, n, x);
 
     if (0 == numIntPoints%2) {
@@ -1953 +1955 @@
 
     n = spline->n;
-    binNum = VectorBinDisectF32(spline->domains, (spline->n)+1, x);
+    binNum = p_psVectorBinDisectF32(spline->domains, (spline->n)+1, x);
     if (binNum == -1) {
         psLogMsg(__func__, PS_LOG_WARN,

trunk/psLib/src/math/psSpline.h

@@ -12 +12 @@
 *  @author George Gusciora, MHPCC
 *
-*  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2004-09-22 01:30:21 $
+*  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2004-09-23 04:56:41 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -407 +407 @@
                           psVector *data);
 
+int p_psVectorBinDisectF32(float *bins,
+                           int numBins,
+                           float x);
+
+int p_psVectorBinDisectS32(int *bins,
+                           int numBins,
+                           int x);
+
 int p_psVectorBinDisect(psVector *bins,
                         psScalar *x);
@@ -415 +423 @@
                                 psScalar *x);
 
+float p_psNRSpline1DEval(psSpline1D *spline,
+                         const psVector* restrict x,
+                         const psVector* restrict y,
+                         float X);
+
 /* \} */// End of MathGroup Functions