Changeset 5180

Timestamp:

Sep 29, 2005, 10:06:58 AM (21 years ago)

Author:

gusciora

Message:

Misc mods to psSPline.c

Location:

trunk/psLib/src/math

Files:

• psPolynomial.c (modified) (2 diffs)
• psSpline.c (modified) (12 diffs)

trunk/psLib/src/math/psPolynomial.c

@@ -7 +7 @@
 *  polynomials.  It also contains a Gaussian functions.
 *
-*  @version $Revision: 1.128 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-09-23 23:01:30 $
+*  @version $Revision: 1.129 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-09-29 20:06:58 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -784 +784 @@
 }
 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 psF64 psPolynomial1DEval(const psPolynomial1D* poly,
                          psF64 x)

trunk/psLib/src/math/psSpline.c

@@ -7 +7 @@
 *  splines.
 *
-*  @version $Revision: 1.130 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2005-09-29 19:40:59 $
+*  @version $Revision: 1.131 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2005-09-29 20:06:58 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -79 +79 @@
 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.  They can be modified via ypo and yp1.
+The first and second derivatives at the endpoints were previously undefined in
+the SDR.  From bugzilla #???, they are required to be 0.0, impementing natural
+splines.
+ 
+Endpoints are defined by
+    PS_LEFT_SPLINE_DERIV
+    PS_RIGHT_SPLINE_DERIV
  
 This routine assumes that vectors x and y are of the appropriate types/sizes
@@ -89 +94 @@
 XXX: do an F64 version?
  *****************************************************************************/
-static psF32 *calculateSecondDerivs(const psVector* x,        ///< Ordinates
-                                    const psVector* y)        ///< Coordinates
+#define PS_LEFT_SPLINE_DERIV 0.0
+#define PS_RIGHT_SPLINE_DERIV 0.0
+static psF32 *calculateSecondDerivs(
+    const psVector* x,                  ///< Ordinates
+    const psVector* y)                  ///< Coordinates
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     if (psTraceGetLevel(__func__) >= 6) {
-        PS_VECTOR_PRINT_F32(x);
-        PS_VECTOR_PRINT_F32(y);
+        p_psVectorPrint(1, (psVector *) x, "x");
+        p_psVectorPrint(1, (psVector *) y, "y");
     }
     psS32 n = y->n;
@@ -339 +347 @@
 
     psS32 binNum;
-    unsigned int numIntPoints = order+1;
+    psS32 numIntPoints = order+1;
     psS32 origin;
 
@@ -460 +468 @@
     return(tmpSpline);
 }
-
 
 /*****************************************************************************
@@ -481 +488 @@
  
 XXX: What types must be supported?
- 
-XXX: Should we allow x==NULL?
  *****************************************************************************/
 psSpline1D *psVectorFitSpline1D(
-    const psVector* x,        ///< Ordinates.
-    const psVector* y)         ///< Coordinates.
+    const psVector* x,                  ///< Ordinates.
+    const psVector* y)                  ///< Coordinates.
 {
     psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
@@ -498 +503 @@
     psVector *xPtr = NULL;
     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);
+        // Convert x to F32 if necessary.
         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];
-            }
+            xPtr = psVectorCopy(NULL, x, PS_TYPE_F32);
         } else if (PS_TYPE_F32 == x->type.type) {
             xPtr = (psVector *) x;
@@ -513 +515 @@
         // 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;
     }
 
@@ -520 +520 @@
     // Convert y to F32 if necessary.
     if (PS_TYPE_F64 == y->type.type) {
-        yPtr = psVectorAlloc(y->n, PS_TYPE_F32);
-        for (psS32 i = 0 ; i < y->n ; i++) {
-            yPtr->data.F32[i] = (psF32) y->data.F64[i];
-        }
+        yPtr = psVectorCopy(NULL, y, PS_TYPE_F32);
     } else {
         yPtr = (psVector *) y;
@@ -570 +567 @@
         tmp+= xPtr->data.F32[i] / H;
         tmp*= spline->p_psDeriv2[i+1] * H * H / 6.0;
-        (spline->spline[i]->coeff[0])+= tmp;
+        spline->spline[i]->coeff[0]+= tmp;
 
         //
@@ -594 +591 @@
         //
         // 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);
 
         //
@@ -602 +599 @@
         //
         // 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);
+        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]);