Changeset 1846 for trunk/psLib/src/math/psMinimize.c

Timestamp:

Sep 21, 2004, 3:30:21 PM (22 years ago)

Author:

gusciora

Message:

Added the spline generation.

File:

• trunk/psLib/src/math/psMinimize.c (modified) (8 diffs)

trunk/psLib/src/math/psMinimize.c

@@ -1 @@
-/** @file  psMinimize.c
+M/** @file  psMinimize.c
  *  \brief basic minimization functions
  *  @ingroup Math
@@ -9 +9 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.41 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-20 23:16:10 $
+ *  @version $Revision: 1.42 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-22 01:30:21 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -148 +148 @@
 }
 
+/*****************************************************************************
+ 
+CalculateSecondDerivs(): Given a set of x/y vectors corresponding to a
+tabulated function at n points, this routine calcualtes 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.
+ *****************************************************************************/
+float *CalculateSecondDerivs(const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
+                             const psVector* restrict y)        ///< Coordinates
+{
+    int i;
+    int k;
+    bool mustFreeX = false;
+    float sig;
+    float p;
+    int n = y->n;
+    float *u = (float *) psAlloc(n * sizeof(float));
+    float *derivs2 = (float *) psAlloc(n * sizeof(float));
+    float *X = (float *) & (x->data.F32[0]);
+    float *Y = (float *) & (y->data.F32[0]);
+
+    if (x == NULL) {
+        X = (float *) psAlloc(n * sizeof(float));
+        for(i=0;i<n;i++) {
+            X[i] = (float) i;
+        }
+        mustFreeX = true;
+    }
+
+    // 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;
+
+    for (i=1;i<=(n-2);i++) {
+        sig = (X[i] - X[i-1]) / (X[i+1] - X[i-1]);
+        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;
+    }
+
+    for (k=(n-2);k>=0;k--) {
+        derivs2[k] = derivs2[k] * derivs2[k+1] + u[k];
+    }
+
+
+
+
+
+
+
+
+
+    /*
+        if (mustFreeX == true) {
+            psFree(X);
+        }
+    */
+    psFree(u);
+
+    return(derivs2);
+}
+
+
 /*****************************************************************************/
 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
 /*****************************************************************************/
+
+
+
+/*****************************************************************************
+psVectorFitSpline1D(): given a psSpline1D data structure and a set of x/y
+vectors, this routine generates the cublic splines which satisfy those data
+points.
+ 
+The formula for calculating the spline polynomials is derived from Numerical
+Recipes in C.  The basic idea is that the polynomial is
+ (1)     y = (A * y[0]) +
+ (2)         (B * y[1]) +
+ (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
+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.
+ 
+XXX: usage of yErr is not specified in IfA documentation.
+ *****************************************************************************/
+
+psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline,              ///< The spline which will be generated.
+                                const psVector* restrict x,        ///< Ordinates (or NULL to just use the indices)
+                                const psVector* restrict y,        ///< Coordinates
+                                const psVector* restrict yErr)     ///< Errors in coordinates, or NULL
+{
+    int numSplines = (y->n)-1;
+    float tmp;
+    float H;
+    int i;
+
+    if (mySpline == NULL) {
+        //XXX psErrorMsg()
+    }
+    mySpline->p_psDeriv2 = CalculateSecondDerivs(x, y);
+
+    for (i=0;i<numSplines;i++) {
+        H = x->data.F32[i+1] - x->data.F32[i];
+        //
+        // ******** Calulate 0-order term ********
+        //
+        // From (1)
+        (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);
+        // From (3)
+        tmp = (x->data.F32[i+1] * x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H);
+        tmp-= (x->data.F32[i+1] / H);
+        tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0;
+        ((mySpline->spline[i])->coeff[0])+= tmp;
+        // From (4)
+        tmp = -(x->data.F32[i] * x->data.F32[i] * x->data.F32[i]) / (H * H * H);
+        tmp+= (x->data.F32[i] / H);
+        tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0;
+        ((mySpline->spline[i])->coeff[0])+= tmp;
+
+        //
+        // ******** Calulate 1-order term ********
+        //
+        // From (1)
+        (mySpline->spline[i])->coeff[1] = - ((mySpline->p_psDeriv2)[i]) / H;
+        // From (2)
+        (mySpline->spline[i])->coeff[1]-= ((mySpline->p_psDeriv2)[i+1]) / H;
+        // From (3)
+        tmp = - (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;
+        // From (4)
+        tmp = (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;
+
+        //
+        // ******** Calulate 2-order term ********
+        //
+        // From (3)
+        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x->data.F32[i+1] / 6.0;
+        // From (4)
+        (mySpline->spline[i])->coeff[2]-= ((mySpline->p_psDeriv2)[i+1]) * 3.0 * x->data.F32[i] / 6.0;
+
+        //
+        // ******** Calulate 3-order term ********
+        //
+        // From (3)
+        (mySpline->spline[i])->coeff[2] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
+        // From (4)
+        (mySpline->spline[i])->coeff[2]+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
+    }
+
+
+
+
+
+
+
+
+
+    return(mySpline);
+}
 
 /******************************************************************************
@@ -388 +561 @@
 
 
+
+
+
+
+
+
+
+
+
+
+
 /******************************************************************************
 psVectorFitPolynomial1D():  This routine must fit a polynomial of degree
@@ -504 +688 @@
     return(min);
 }
+
 
 #define STEP_SIZE 0.10
@@ -551 +736 @@
         }
     }
+
     if (null == 0) {
         psTrace(".psLib.dataManip.p_psDetermineBracket", 2,
@@ -685 +871 @@
     return(NULL);
 }
+
+
+
+
+
+
+
+
+
+
 
 
@@ -865 +1061 @@
     return(0.0);
 }
+
+
+
+
+
+
+
+
+
+
+
+
+