Changeset 1858

Timestamp:

Sep 22, 2004, 6:55:48 PM (22 years ago)

Author:

gusciora

Message:

...

File:

• trunk/psLib/test/dataManip/tst_psFunc07.c (modified) (3 diffs)

trunk/psLib/test/dataManip/tst_psFunc07.c

@@ -15 +15 @@
 #include "psFunctions.h"
 
-#define N 8
+#define NUM_SPLINES 4
 #define A 3.0
 #define B -4.0
 #define C -5.0
 #define D 6.0
+#define MIN 1.0
+#define MAX 30
 
 float myFunc(float x)
 {
     return(A + x * (B + x * (C + x * (D))));
+    //    return(A + x * (B + x * (C)));
+    //    return(A + x * (B));
+    //    return(A);
 }
+
+
+/*****************************************************************************
+XXX: This is Numerical Recipes in C code.  Remove from distribution.
+ *****************************************************************************/
+void spline(float *x, float *y, int n, float yp1, float ypn, float *y2)
+{
+    int i,k;
+    float p,qn,sig,un,u[NUM_SPLINES+1];
+
+    if (yp1 > 0.99e30) {
+        y2[1]=u[1]=0.0;
+    } else {
+        y2[1] = -0.5;
+        u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
+    }
+
+    for (i=2;i<=n-1;i++) {
+        sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
+        p=sig*y2[i-1]+2.0;
+        y2[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;
+    }
+
+    if (ypn > 0.99e30) {
+        qn=un=0.0;
+    } else {
+        qn=0.5;
+        un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
+    }
+    y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
+
+    for (k=n-1;k>=1;k--) {
+        y2[k]=y2[k]*y2[k+1]+u[k];
+    }
+}
+
+/*****************************************************************************
+XXX: This is Numerical Recipes in C code.  Remove from distribution.
+ *****************************************************************************/
+void splint(float xa[], float ya[], float y2a[], int n, float x, float *y)
+{
+    int klo,khi,k;
+    float h,b,a;
+
+    klo=1;
+    khi=n;
+    while (khi-klo > 1) {
+        k=(khi+klo) >> 1;
+        if (xa[k] > x)
+            khi=k;
+        else
+            klo=k;
+    }
+    h=xa[khi]-xa[klo];
+    a=(xa[khi]-x)/h;
+    b=(x-xa[klo])/h;
+    *y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
+}
+
 
 int main()
@@ -33 +99 @@
     //    int  currentId = psMemGetId();
     psSpline1D *tmpSpline;
-    psVector *data;
     psVector *x;
+    psVector *newX;
     psVector *y;
+    psVector *newY;
+    float NRX[NUM_SPLINES+2];
+    float NRY[NUM_SPLINES+2];
+    float tmpY;
+    float derivs[NUM_SPLINES+2];
 
     psTraceSetLevel(".", 0);
@@ -44 +115 @@
                             "psFunction functions",
                             "psSpline1DAlloc()");
-    data = psVectorAlloc(N+1, PS_TYPE_F32);
-    x = psVectorAlloc(N+1, PS_TYPE_F32);
-    y = psVectorAlloc(N+1, PS_TYPE_F32);
-    tmpSpline = psSpline1DAlloc(N, 3, 1.0, 10.0);
+    x = psVectorAlloc(NUM_SPLINES+1, PS_TYPE_F32);
+    newX = psVectorAlloc(NUM_SPLINES+1, PS_TYPE_F32);
+    y = psVectorAlloc(NUM_SPLINES+1, PS_TYPE_F32);
+    tmpSpline = psSpline1DAlloc(NUM_SPLINES, 3, MIN, MAX);
 
-    for (i=0;i<N+1;i++) {
-        data->data.F32[i] = tmpSpline->domains[i];
-        x->data.F32[i] = tmpSpline->domains[i];
-        y->data.F32[i] = tmpSpline->domains[i];
-        y->data.F32[i] = myFunc(x->data.F32[i]);
-        printf("HMMM: (x, y) is (%f, %f)\n", x->data.F32[i], y->data.F32[i]);
+    for (i=0;i<NUM_SPLINES+1;i++) {
+        NRX[i+1] = x->data.F32[i] = tmpSpline->domains[i];
+        NRY[i+1] = y->data.F32[i] = myFunc(x->data.F32[i]);
+        printf("DATA: (x, y) is (%f, %f)\n", x->data.F32[i], y->data.F32[i]);
+    }
+    for (i=0;i<NUM_SPLINES;i++) {
+        newX->data.F32[i]= ((x->data.F32[i] + x->data.F32[i+1])/2.0);
+    }
+    newX->data.F32[NUM_SPLINES] = x->data.F32[NUM_SPLINES];
+    /****************************************************************************/
+    /*   NR Code       */
+    /****************************************************************************/
+    printf("========================== NR Code =========================\n");
+    spline(&NRX[0],
+           &NRY[0],
+           NUM_SPLINES+1,
+           0.0,
+           0.0,
+           &derivs[0]);
+    for (i=1;i<=NUM_SPLINES+1;i++) {
+        printf("NR: derivs[%d] is %f\n", i, derivs[i]);
     }
 
-    //    psSpline1DGen(tmpSpline, data);
-    psVectorFitSpline1D(tmpSpline, x, y, NULL);
+    for (i=1;i<=NUM_SPLINES+1;i++) {
+        printf("===================================================\n");
+        splint(&NRX[0],
+               &NRY[0],
+               &derivs[0],
+               NUM_SPLINES+1,
+               newX->data.F32[i-1],
+               &tmpY);
 
-    for (i=0;i<N+1;i++) {
-        //        x->data.F32[i] = 0.5 + (float) i+1;
+        printf("NR: p_psNRSpline1DEval(%f) is %f.  Should be %f\n",
+               newX->data.F32[i-1], tmpY, myFunc(newX->data.F32[i-1]));
     }
 
-    y = psSpline1DEvalVector(x, tmpSpline);
 
-    for (i=0;i<N+1;i++) {
-        if (fabs(x->data.F32[i]-y->data.F32[i]) > FLT_EPSILON) {
-            printf("ERROR: f(%f) is %f.  Should be %f\n", x->data.F32[i], y->data.F32[i], myFunc(x->data.F32[i]));
-            testStatus = true;
-        } else {
-            printf("COOL: f(%f) is %f\n", x->data.F32[i], y->data.F32[i]);
+    /****************************************************************************/
+    /*   psLib Code      */
+    /****************************************************************************/
+    psVectorFitSpline1D(tmpSpline, x, y, NULL);
+    for (i=0;i<NUM_SPLINES+1;i++) {
+        printf("psLib: derivs[%d] is %f\n", i, tmpSpline->p_psDeriv2[i]);
+    }
+    newY = psSpline1DEvalVector(newX, tmpSpline);
 
-        }
+    for (i=0;i<NUM_SPLINES+1;i++) {
+        printf("===================================================\n");
+        printf("EVAL: f(%f) is %f.  Should be %f\n",
+               newX->data.F32[i], newY->data.F32[i], myFunc(newX->data.F32[i]));
+
+        tmpY = p_psNRSpline1DEval(tmpSpline, x, y, newX->data.F32[i]);
+        printf("HMMM: p_psNRSpline1DEval(%f) is %f.  Should be %f\n",
+               newX->data.F32[i], tmpY, myFunc(newX->data.F32[i]));
     }