Index: trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- trunk/psLib/src/dataManip/psMinimize.c (revision 1836) +++ trunk/psLib/src/dataManip/psMinimize.c (revision 1846) @@ -1,3 +1,3 @@ -/** @file psMinimize.c +M/** @file psMinimize.c * \brief basic minimization functions * @ingroup Math @@ -9,6 +9,6 @@ * @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,7 +148,180 @@ } +/***************************************************************************** + +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=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;idata.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,4 +561,15 @@ + + + + + + + + + + + /****************************************************************************** psVectorFitPolynomial1D(): This routine must fit a polynomial of degree @@ -504,4 +688,5 @@ return(min); } + #define STEP_SIZE 0.10 @@ -551,4 +736,5 @@ } } + if (null == 0) { psTrace(".psLib.dataManip.p_psDetermineBracket", 2, @@ -685,4 +871,14 @@ return(NULL); } + + + + + + + + + + @@ -865,4 +1061,17 @@ return(0.0); } + + + + + + + + + + + + +