IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 11, 2005, 12:18:40 PM (21 years ago)
Author:
gusciora
Message:

....

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/psLib/src/math/psSpline.c

    r4971 r4991  
    77*  splines.
    88*
    9 *  @version $Revision: 1.122 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2005-09-08 00:14:32 $
     9*  @version $Revision: 1.123 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2005-09-11 22:18:40 $
    1111*
    1212*
     
    207207}
    208208
     209
     210/*****************************************************************************
     211CalculateSecondDerivs(): Given a set of x/y vectors corresponding to a
     212tabulated function at n points, this routine calculates the second
     213derivatives of the interpolating cubic splines at those n points.
     214 
     215The first and second derivatives at the endpoints, undefined in the SDR, are
     216here defined to be 0.0.  They can be modified via ypo and yp1.
     217 
     218This routine assumes that vectors x and y are of the appropriate types/sizes
     219(F32).
     220 
     221XXX: This algorithm is derived from the Numerical Recipes.
     222XXX: use recycled vectors for internal data.
     223XXX: do an F64 version?
     224 *****************************************************************************/
     225static psF32 *calculateSecondDerivs(const psVector* x,        ///< Ordinates (or NULL to just use the indices)
     226                                    const psVector* y)        ///< Coordinates
     227{
     228    psTrace(".psLib.dataManip.calculateSecondDerivs", 4,
     229            "---- calculateSecondDerivs() begin ----\n");
     230
     231    psS32 i;
     232    psS32 k;
     233    psF32 sig;
     234    psF32 p;
     235    psS32 n = y->n;
     236    psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32));
     237    psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32));
     238    psF32 *X = (psF32 *) & (x->data.F32[0]);
     239    psF32 *Y = (psF32 *) & (y->data.F32[0]);
     240    psF32 qn;
     241
     242    // XXX: The second derivatives at the endpoints, undefined in the SDR,
     243    // are set in psConstants.h: PS_LEFT_SPLINE_DERIV, PS_RIGHT_SPLINE_DERIV.
     244    derivs2[0] = -0.5;
     245    u[0]= (3.0/(X[1]-X[0])) * ((Y[1]-Y[0])/(X[1]-X[0]) - PS_LEFT_SPLINE_DERIV);
     246
     247    for (i=1;i<=(n-2);i++) {
     248        sig = (X[i] - X[i-1]) / (X[i+1] - X[i-1]);
     249        p = sig * derivs2[i-1] + 2.0;
     250        derivs2[i] = (sig - 1.0) / p;
     251        u[i] = ((Y[i+1] - Y[i])/(X[i+1]-X[i])) - ((Y[i]-Y[i-1])/(X[i]-X[i-1]));
     252        u[i] = ((6.0 * u[i] / (X[i+1] - X[i-1])) - (sig * u[i-1])) / p;
     253
     254        psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
     255                "X[%d] is %f\n", i, X[i]);
     256        psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
     257                "Y[%d] is %f\n", i, Y[i]);
     258        psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
     259                "u[%d] is %f\n", i, u[i]);
     260    }
     261
     262    qn = 0.5;
     263    u[n-1] = (3.0/(X[n-1]-X[n-2])) * (PS_RIGHT_SPLINE_DERIV - (Y[n-1]-Y[n-2])/(X[n-1]-X[n-2]));
     264    derivs2[n-1] = (u[n-1] - (qn * u[n-2])) / ((qn * derivs2[n-2]) + 1.0);
     265
     266    for (k=(n-2);k>=0;k--) {
     267        derivs2[k] = derivs2[k] * derivs2[k+1] + u[k];
     268
     269        psTrace(".psLib.dataManip.calculateSecondDerivs", 6,
     270                "derivs2[%d] is %f\n", k, derivs2[k]);
     271    }
     272
     273    psFree(u);
     274    psTrace(".psLib.dataManip.calculateSecondDerivs", 4,
     275            "---- calculateSecondDerivs() end ----\n");
     276    return(derivs2);
     277}
     278
     279/*****************************************************************************/
     280/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
     281/*****************************************************************************/
     282
     283/*****************************************************************************
     284psVectorFitSpline1D(): given a psSpline1D data structure and a set of x/y
     285vectors, this routine generates the linear or cublic splines which satisfy
     286those data points.
     287 
     288The formula for calculating the spline polynomials is derived from Numerical
     289Recipes in C.  The basic idea is that the polynomial is
     290 (1)     y = (A * y[0]) +
     291 (2)         (B * y[1]) +
     292 (3)         ((((A*A*A)-A) * mySpline->p_psDeriv2[0]) * H^2)/6.0 +
     293 (4)         ((((B*B*B)-B) * mySpline->p_psDeriv2[1]) * H^2)/6.0
     294Where:
     295 H = x[1]-x[0]
     296 A = (x[1]-x)/H
     297 B = (x-x[0])/H
     298The bulk of the code in this routine is the expansion of the above equation
     299into a polynomial in terms of x, and then saving the coefficients of the
     300powers of x in the spline polynomials.  This gets pretty complicated.
     301 
     302XXX: usage of yErr is not specified in IfA documentation.
     303 
     304XXX: Is the x argument redundant?  What do we do if the x argument is
     305supplied, but does not equal the knots specified in mySpline?
     306 
     307XXX: can psSpline be NULL?
     308 
     309XXX: reimplement this assuming that mySpline is NULL?
     310 
     311XXX: What happens if X is NULL, then an index vector is generated for X, but
     312that index vector lies outside the range vectors in mySpline?
     313 
     314XXX: Assumes mySpline->knots is psF32.  Must add psU32 and psF64.
     315 *****************************************************************************/
     316psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline,     ///< The spline which will be generated.
     317                                const psVector* x,        ///< Ordinates (or NULL to just use the indices)
     318                                const psVector* y,        ///< Coordinates
     319                                const psVector* yErr)     ///< Errors in coordinates, or NULL
     320{
     321    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
     322    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
     323    if (mySpline != NULL) {
     324        PS_ASSERT_VECTOR_TYPE(mySpline->knots, PS_TYPE_F32, NULL);
     325    }
     326    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
     327            "---- psVectorFitSpline1D() begin ----\n");
     328    psS32 numSplines = (y->n)-1;
     329    psF32 tmp;
     330    psF32 H;
     331    psS32 i;
     332    psF32 slope;
     333    psVector *x32 = NULL;
     334    psVector *y32 = NULL;
     335    psVector *yErr32 = NULL;
     336    static psVector *x32Static = NULL;
     337    static psVector *y32Static = NULL;
     338    static psVector *yErr32Static = NULL;
     339
     340    PS_VECTOR_CONVERT_F64_TO_F32_STATIC(y, y32, y32Static);
     341
     342    // If yErr==NULL, set all errors equal.
     343    if (yErr == NULL) {
     344        PS_VECTOR_GEN_YERR_STATIC_F32(yErr32Static, y->n);
     345        yErr32 = yErr32Static;
     346    } else {
     347        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(yErr, NULL);
     348        PS_VECTOR_CONVERT_F64_TO_F32_STATIC(yErr, yErr32, yErr32Static);
     349    }
     350
     351    // If x==NULL, create an x32 vector with x values set to (0:n).
     352    if (x == NULL) {
     353        PS_VECTOR_GEN_X_INDEX_STATIC_F32(x32Static, y->n);
     354        x32 = x32Static;
     355    } else {
     356        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
     357        PS_VECTOR_CONVERT_F64_TO_F32_STATIC(x, x32, x32Static);
     358    }
     359    PS_ASSERT_VECTORS_SIZE_EQUAL(x32, y32, NULL);
     360    PS_ASSERT_VECTORS_SIZE_EQUAL(yErr32, y32, NULL);
     361
     362    /*
     363        XXX:
     364        This can not be implemented until SDR states what order spline should be
     365        created.
     366        Should we error if mySpline is not NULL?
     367        Should we error if mySPline is not NULL?
     368    */
     369    if (mySpline == NULL) {
     370        mySpline = psSpline1DAllocGeneric(x32, 3);
     371    }
     372    PS_ASSERT_PTR_NON_NULL(mySpline, NULL);
     373    PS_ASSERT_INT_NONNEGATIVE(mySpline->n, NULL);
     374
     375    if (y32->n != (1 + mySpline->n)) {
     376        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
     377                "data size / spline size mismatch (%d %d)\n",
     378                y32->n, mySpline->n);
     379        return(NULL);
     380    }
     381
     382    // If these are linear splines, which means their polynomials will have
     383    // two coefficients, then we do the simple calculation.
     384    if (2 == (mySpline->spline[0])->n) {
     385        for (i=0;i<mySpline->n;i++) {
     386            slope = (y32->data.F32[i+1] - y32->data.F32[i]) /
     387                    (mySpline->knots->data.F32[i+1] - mySpline->knots->data.F32[i]);
     388            (mySpline->spline[i])->coeff[0] = y32->data.F32[i] -
     389                                              (slope * mySpline->knots->data.F32[i]);
     390
     391            (mySpline->spline[i])->coeff[1] = slope;
     392            psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
     393                    "---- mySpline %d coeffs are (%f, %f)\n", i,
     394                    (mySpline->spline[i])->coeff[0],
     395                    (mySpline->spline[i])->coeff[1]);
     396        }
     397        psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
     398                "---- Exiting psVectorFitSpline1D()()\n");
     399        return((psSpline1D *) mySpline);
     400    }
     401
     402    // Check if these are cubic splines (n==4).  If not, psError.
     403    if (4 != (mySpline->spline[0])->n) {
     404        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
     405                "Don't know how to generate %d-order splines.",
     406                (mySpline->spline[0])->n-1);
     407        return(NULL);
     408    }
     409
     410    // If we get here, then we know these are cubic splines.  We first
     411    // generate the second derivatives at each data point.
     412    mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32);
     413    for (i=0;i<y32->n;i++)
     414        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     415                "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
     416
     417    // We generate the coefficients of the spline polynomials.  I can't
     418    // concisely explain how this code works.  See above function comments
     419    // and Numerical Recipes in C.
     420    for (i=0;i<numSplines;i++) {
     421        H = x32->data.F32[i+1] - x32->data.F32[i];
     422        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
     423                "x data (%f - %f) (%f)\n",
     424                x32->data.F32[i],
     425                x32->data.F32[i+1], H);
     426        //
     427        // ******** Calculate 0-order term ********
     428        //
     429        // From (1)
     430        (mySpline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H);
     431        // From (2)
     432        ((mySpline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H);
     433        // From (3)
     434        tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
     435        tmp-= (x32->data.F32[i+1] / H);
     436        tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0;
     437        ((mySpline->spline[i])->coeff[0])+= tmp;
     438        // From (4)
     439        tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
     440        tmp+= (x32->data.F32[i] / H);
     441        tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0;
     442        ((mySpline->spline[i])->coeff[0])+= tmp;
     443
     444        //
     445        // ******** Calculate 1-order term ********
     446        //
     447        // From (1)
     448        (mySpline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H;
     449        // From (2)
     450        ((mySpline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H);
     451        // From (3)
     452        tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
     453        tmp+= (1.0 / H);
     454        tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0;
     455        ((mySpline->spline[i])->coeff[1])+= tmp;
     456        // From (4)
     457        tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
     458        tmp-= (1.0 / H);
     459        tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
     460        ((mySpline->spline[i])->coeff[1])+= tmp;
     461
     462        //
     463        // ******** Calculate 2-order term ********
     464        //
     465        // From (3)
     466        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H);
     467        // From (4)
     468        ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H));
     469
     470        //
     471        // ******** Calculate 3-order term ********
     472        //
     473        // From (3)
     474        (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
     475        // From (4)
     476        ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
     477
     478        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     479                "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
     480        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     481                "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
     482        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     483                "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
     484        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     485                "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
     486
     487    }
     488
     489    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
     490            "---- psVectorFitSpline1D() end ----\n");
     491    return(mySpline);
     492}
     493
     494
     495
     496
     497
     498
     499
     500
     501
     502/*****************************************************************************
     503psVectorFitSpline1DNEW(): given a psSpline1D data structure and a set of x/y
     504 
     505xF32 and yF32 are internal psVectors which are used to hold the psF32 versions
     506of the input data, if necessary.  xPtr and yPtr are pointers to either xF32 or
     507the x argument.  All computation is done on xPtr and yPtr.  xF32 and yF32 will
     508simply be psFree() at the end.
     509 
     510XXX: nKnots makes no sense.  This number is always equal to the size of the x
     511an y vectors.
     512 
     513XXX: How do we specify the spline order?  For now, order=3.
     514 *****************************************************************************/
     515#define PS_XXX_SPLINE_ORDER 3
     516psSpline1D *psVectorFitSpline1DNEW(const psVector* x,        ///< Ordinates (or NULL to just use the indices)
     517                                   const psVector* y,        ///< Coordinates
     518                                   int nKnots)
     519{
     520    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, "---- psVectorFitSpline1DNEW() begin ----\n");
     521    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
     522    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
     523    //    PS_ASSERT_INT_EQUAL(y->n, nKnots);
     524
     525    //
     526    // The following code ensures that xPtr points to a psF32 version of the
     527    // ordinate data.
     528    //
     529    psVector *xF32 = NULL;
     530    psVector *xPtr = NULL;
     531    if (x != NULL) {
     532        PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
     533        if (PS_TYPE_F64 == x->type.type) {
     534            xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
     535            for (psS32 i = 0 ; i < x->n ; i++) {
     536                xF32->data.F32[i] = (psF32) x->data.F64[i];
     537            }
     538            xPtr = xF32;
     539        } else if (PS_TYPE_F32 == x->type.type) {
     540            xPtr = (psVector *) x;
     541        } else {
     542            psError(PS_ERR_UNKNOWN, true, "psVector x is wrong type.\n");
     543            return(NULL);
     544        }
     545    } else {
     546        // If x==NULL, create an x32 vector with x values set to (0:n).
     547        xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
     548        for (psS32 i = 0 ; i < x->n ; i++) {
     549            xF32->data.F32[i] = (psF32) i;
     550        }
     551        xPtr = xF32;
     552    }
     553
     554    //
     555    // If y is of type psF64, then create a new vector yF32 and convert the
     556    // y elements.  Regardless of y's type, we create a yPtr which will be
     557    // used in the remainder of this function.
     558    //
     559    psVector *yF32 = NULL;
     560    psVector *yPtr = NULL;
     561    if (PS_TYPE_F64 == y->type.type) {
     562        yF32 = psVectorAlloc(y->n, PS_TYPE_F32);
     563        for (psS32 i = 0 ; i < y->n ; i++) {
     564            yF32->data.F32[i] = (psF32) y->data.F64[i];
     565        }
     566        yPtr = yF32;
     567    } else {
     568        yPtr = (psVector *) y;
     569    }
     570
     571    psSpline1D *mySpline = psSpline1DAllocGeneric(xPtr, PS_XXX_SPLINE_ORDER);
     572
     573    psS32 numSplines = nKnots - 1;
     574    psF32 tmp;
     575    psF32 H;
     576    psS32 i;
     577    psF32 slope;
     578    // XXX: get rid of x32 and y32 (this is from old code)
     579    psVector *x32 = xPtr;
     580    psVector *y32 = yPtr;
     581
     582    // If these are linear splines, which means their polynomials will have
     583    // two coefficients, then we do the simple calculation.
     584    if (1 == PS_XXX_SPLINE_ORDER) {
     585        for (i=0;i<mySpline->n;i++) {
     586            slope = (y32->data.F32[i+1] - y32->data.F32[i]) /
     587                    (mySpline->knots->data.F32[i+1] - mySpline->knots->data.F32[i]);
     588            (mySpline->spline[i])->coeff[0] = y32->data.F32[i] -
     589                                              (slope * mySpline->knots->data.F32[i]);
     590
     591            (mySpline->spline[i])->coeff[1] = slope;
     592            psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
     593                    "---- mySpline %d coeffs are (%f, %f)\n", i,
     594                    (mySpline->spline[i])->coeff[0],
     595                    (mySpline->spline[i])->coeff[1]);
     596        }
     597        psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
     598                "---- Exiting psVectorFitSpline1D()()\n");
     599        return((psSpline1D *) mySpline);
     600    }
     601
     602    //
     603    // Check if these are cubic splines (n==4).  If not, psError.
     604    //
     605    if (3 != PS_XXX_SPLINE_ORDER) {
     606        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
     607                "Don't know how to generate %d-order splines.", PS_XXX_SPLINE_ORDER);
     608        return(NULL);
     609    }
     610
     611    //
     612    // If we get here, then we know these are cubic splines.  We first
     613    // generate the second derivatives at each data point.
     614    //
     615    mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32);
     616    for (i=0;i<y32->n;i++)
     617        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     618                "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
     619
     620    //
     621    // We generate the coefficients of the spline polynomials.  I can't
     622    // concisely explain how this code works.  See above function comments
     623    // and Numerical Recipes in C.
     624    //
     625    for (i=0;i<numSplines;i++) {
     626        H = x32->data.F32[i+1] - x32->data.F32[i];
     627        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
     628                "x data (%f - %f) (%f)\n",
     629                x32->data.F32[i],
     630                x32->data.F32[i+1], H);
     631        //
     632        // ******** Calculate 0-order term ********
     633        //
     634        // From (1)
     635        (mySpline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H);
     636        // From (2)
     637        ((mySpline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H);
     638        // From (3)
     639        tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
     640        tmp-= (x32->data.F32[i+1] / H);
     641        tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0;
     642        ((mySpline->spline[i])->coeff[0])+= tmp;
     643        // From (4)
     644        tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
     645        tmp+= (x32->data.F32[i] / H);
     646        tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0;
     647        ((mySpline->spline[i])->coeff[0])+= tmp;
     648
     649        //
     650        // ******** Calculate 1-order term ********
     651        //
     652        // From (1)
     653        (mySpline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H;
     654        // From (2)
     655        ((mySpline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H);
     656        // From (3)
     657        tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
     658        tmp+= (1.0 / H);
     659        tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0;
     660        ((mySpline->spline[i])->coeff[1])+= tmp;
     661        // From (4)
     662        tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
     663        tmp-= (1.0 / H);
     664        tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
     665        ((mySpline->spline[i])->coeff[1])+= tmp;
     666
     667        //
     668        // ******** Calculate 2-order term ********
     669        //
     670        // From (3)
     671        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H);
     672        // From (4)
     673        ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H));
     674
     675        //
     676        // ******** Calculate 3-order term ********
     677        //
     678        // From (3)
     679        (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
     680        // From (4)
     681        ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
     682
     683        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     684                "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
     685        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     686                "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
     687        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     688                "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
     689        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
     690                "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
     691
     692    }
     693
     694    psFree(xF32);
     695    psFree(yF32);
     696
     697    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
     698            "---- psVectorFitSpline1D() end ----\n");
     699
     700    return(mySpline);
     701}
     702
     703
     704
     705
     706
    209707/*****************************************************************************/
    210708/*  FUNCTION IMPLEMENTATION - PUBLIC                                         */
     
    283781        }
    284782    } else {
    285         printf("XXX: Generate an error here.\n");
     783        psError(PS_ERR_UNKNOWN, true, "psVector in has an incorrect type.\n");
    286784        return(NULL);
    287785    }
Note: See TracChangeset for help on using the changeset viewer.