IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 19, 2005, 3:21:13 PM (21 years ago)
Author:
drobbin
Message:

changed psPolynomial fxns to use F64 and removed psDPoly fxns

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/math/psFunctions.c

    r4580 r4581  
    77*  polynomials.  It also contains a Gaussian functions.
    88*
    9 *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2005-07-19 02:55:54 $
     9*  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2005-07-20 01:21:13 $
    1111*
    1212*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    5050static void polynomial3DFree(psPolynomial3D* poly);
    5151static void polynomial4DFree(psPolynomial4D* poly);
    52 static void dPolynomial1DFree(psDPolynomial1D* poly);
    53 static void dPolynomial2DFree(psDPolynomial2D* poly);
    54 static void dPolynomial3DFree(psDPolynomial3D* poly);
    55 static void dPolynomial4DFree(psDPolynomial4D* poly);
    5652static void spline1DFree(psSpline1D *tmpSpline);
    5753static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x);
    … …  
    166162}
    167163
    168 static void dPolynomial1DFree(psDPolynomial1D* poly)
    169 {
    170     psFree(poly->coeff);
    171     psFree(poly->coeffErr);
    172     psFree(poly->mask);
    173 }
    174 
    175 static void dPolynomial2DFree(psDPolynomial2D* poly)
    176 {
    177     for (unsigned int x = 0; x < poly->nX; x++) {
    178         psFree(poly->coeff[x]);
    179         psFree(poly->coeffErr[x]);
    180         psFree(poly->mask[x]);
    181     }
    182     psFree(poly->coeff);
    183     psFree(poly->coeffErr);
    184     psFree(poly->mask);
    185 }
    186 
    187 static void dPolynomial3DFree(psDPolynomial3D* poly)
    188 {
    189     unsigned int x = 0;
    190     unsigned int y = 0;
    191 
    192     for (x = 0; x < poly->nX; x++) {
    193         for (y = 0; y < poly->nY; y++) {
    194             psFree(poly->coeff[x][y]);
    195             psFree(poly->coeffErr[x][y]);
    196             psFree(poly->mask[x][y]);
    197         }
    198         psFree(poly->coeff[x]);
    199         psFree(poly->coeffErr[x]);
    200         psFree(poly->mask[x]);
    201     }
    202 
    203     psFree(poly->coeff);
    204     psFree(poly->coeffErr);
    205     psFree(poly->mask);
    206 }
    207 
    208 static void dPolynomial4DFree(psDPolynomial4D* poly)
    209 {
    210     unsigned int x = 0;
    211     unsigned int y = 0;
    212     unsigned int z = 0;
    213 
    214     for (x = 0; x < poly->nX; x++) {
    215         for (y = 0; y < poly->nY; y++) {
    216             for (z = 0; z < poly->nZ; z++) {
    217                 psFree(poly->coeff[x][y][z]);
    218                 psFree(poly->coeffErr[x][y][z]);
    219                 psFree(poly->mask[x][y][z]);
    220             }
    221             psFree(poly->coeff[x][y]);
    222             psFree(poly->coeffErr[x][y]);
    223             psFree(poly->mask[x][y]);
    224         }
    225         psFree(poly->coeff[x]);
    226         psFree(poly->coeffErr[x]);
    227         psFree(poly->mask[x]);
    228     }
    229 
    230     psFree(poly->coeff);
    231     psFree(poly->coeffErr);
    232     psFree(poly->mask);
    233 }
    234 
    235164/*****************************************************************************
    236165createChebyshevPolys(n): this routine takes as input the required order n,
    … …  
    283212{
    284213    psS32 loop_x = 0;
    285     psF32 polySum = 0.0;
    286     psF32 xSum = 1.0;
     214    psF64 polySum = 0.0;
     215    psF64 xSum = 1.0;
    287216
    288217    psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4,
    289             "---- Calling ordPolynomial1DEval(%f)\n", x);
     218            "---- Calling ordPolynomial1DEval(%lf)\n", x);
    290219    psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4,
    291220            "Polynomial order is %d\n", poly->n);
    292221    for (loop_x = 0; loop_x < poly->n; loop_x++) {
    293222        psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4,
    294                 "Polynomial coeff[%d] is %f\n", loop_x, poly->coeff[loop_x]);
     223                "Polynomial coeff[%d] is %lf\n", loop_x, poly->coeff[loop_x]);
    295224    }
    296225
    … …  
    298227        if (poly->mask[loop_x] == 0) {
    299228            psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 10,
    300                     "polysum+= sum*coeff [%f+= (%f * %f)\n", polySum, xSum, poly->coeff[loop_x]);
     229                    "polysum+= sum*coeff [%lf+= (%lf * %lf)\n", polySum, xSum, poly->coeff[loop_x]);
    301230            polySum += xSum * poly->coeff[loop_x];
    302231        }
    … …  
    312241static psF64 chebPolynomial1DEval(psF64 x, const psPolynomial1D* poly)
    313242{
    314     PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
     243    PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
    315244    // XXX: Create a macro for this in psConstants.h
    316245    if (poly->n < 1) {
    … …  
    321250    psS32 n = poly->n;
    322251    psS32 i;
    323     psF32 tmp = 0.0;
     252    psF64 tmp = 0.0;
    324253
    325254    // Special case where the Chebyshev poly is constant.
    … …  
    343272
    344273    // General case where the Chebyshev poly has 2 or more terms.
    345     d = psVectorAlloc(n, PS_TYPE_F32);
     274    d = psVectorAlloc(n, PS_TYPE_F64);
    346275    if(poly->mask[n-1] == 0) {
    347         d->data.F32[n-1] = poly->coeff[n-1];
     276        d->data.F64[n-1] = poly->coeff[n-1];
    348277    } else {
    349         d->data.F32[n-1] = 0.0;
    350     }
    351 
    352     d->data.F32[n-2] = (2.0 * x * d->data.F32[n-1]);
     278        d->data.F64[n-1] = 0.0;
     279    }
     280
     281    d->data.F64[n-2] = (2.0 * x * d->data.F64[n-1]);
    353282    if(poly->mask[n-2] == 0) {
    354         d->data.F32[n-2] += poly->coeff[n-2];
     283        d->data.F64[n-2] += poly->coeff[n-2];
    355284    }
    356285
    357286    for (i=n-3;i>=1;i--) {
    358         d->data.F32[i] = (2.0 * x * d->data.F32[i+1]) -
    359                          (d->data.F32[i+2]);
     287        d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) -
     288                         (d->data.F64[i+2]);
    360289        if(poly->mask[i] == 0) {
    361             d->data.F32[i] += poly->coeff[i];
    362         }
    363     }
    364 
    365     tmp = (x * d->data.F32[1]) -
    366           (d->data.F32[2]);
     290            d->data.F64[i] += poly->coeff[i];
     291        }
     292    }
     293
     294    tmp = (x * d->data.F64[1]) -
     295          (d->data.F64[2]);
    367296    if(poly->mask[0] == 0) {
    368297        tmp += (0.5 * poly->coeff[0]);
    … …  
    400329    psS32 loop_x = 0;
    401330    psS32 loop_y = 0;
    402     psF32 polySum = 0.0;
    403     psF32 xSum = 1.0;
    404     psF32 ySum = 1.0;
     331    psF64 polySum = 0.0;
     332    psF64 xSum = 1.0;
     333    psF64 ySum = 1.0;
    405334
    406335    for (loop_x = 0; loop_x < poly->nX; loop_x++) {
    … …  
    420349static psF64 chebPolynomial2DEval(psF64 x, psF64 y, const psPolynomial2D* poly)
    421350{
    422     PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
    423     PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0);
     351    PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
     352    PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0);
    424353    PS_ASSERT_POLY_NON_NULL(poly, NAN);
    425354
    … …  
    427356    psS32 loop_y = 0;
    428357    psS32 i = 0;
    429     psF32 polySum = 0.0;
     358    psF64 polySum = 0.0;
    430359    psPolynomial1D* *chebPolys = NULL;
    431360    psS32 maxChebyPoly = 0;
    … …  
    460389    psS32 loop_y = 0;
    461390    psS32 loop_z = 0;
    462     psF32 polySum = 0.0;
    463     psF32 xSum = 1.0;
    464     psF32 ySum = 1.0;
    465     psF32 zSum = 1.0;
     391    psF64 polySum = 0.0;
     392    psF64 xSum = 1.0;
     393    psF64 ySum = 1.0;
     394    psF64 zSum = 1.0;
    466395
    467396    for (loop_x = 0; loop_x < poly->nX; loop_x++) {
    … …  
    485414static psF64 chebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psPolynomial3D* poly)
    486415{
    487     PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
    488     PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0);
    489     PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0);
     416    PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
     417    PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0);
     418    PS_ASSERT_DOUBLE_WITHIN_RANGE(z, -1.0, 1.0, 0.0);
    490419    psS32 loop_x = 0;
    491420    psS32 loop_y = 0;
    492421    psS32 loop_z = 0;
    493422    psS32 i = 0;
    494     psF32 polySum = 0.0;
     423    psF64 polySum = 0.0;
    495424    psPolynomial1D* *chebPolys = NULL;
    496425    psS32 maxChebyPoly = 0;
    … …  
    533462    psS32 loop_z = 0;
    534463    psS32 loop_t = 0;
    535     psF32 polySum = 0.0;
    536     psF32 xSum = 1.0;
    537     psF32 ySum = 1.0;
    538     psF32 zSum = 1.0;
    539     psF32 tSum = 1.0;
     464    psF64 polySum = 0.0;
     465    psF64 xSum = 1.0;
     466    psF64 ySum = 1.0;
     467    psF64 zSum = 1.0;
     468    psF64 tSum = 1.0;
    540469
    541470    for (loop_x = 0; loop_x < poly->nX; loop_x++) {
    … …  
    563492static psF64 chebPolynomial4DEval(psF64 x, psF64 y, psF64 z, psF64 t, const psPolynomial4D* poly)
    564493{
    565     PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
    566     PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0);
    567     PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0);
    568     PS_ASSERT_FLOAT_WITHIN_RANGE(t, -1.0, 1.0, 0.0);
     494    PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
     495    PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0);
     496    PS_ASSERT_DOUBLE_WITHIN_RANGE(z, -1.0, 1.0, 0.0);
     497    PS_ASSERT_DOUBLE_WITHIN_RANGE(t, -1.0, 1.0, 0.0);
    569498    psS32 loop_x = 0;
    570499    psS32 loop_y = 0;
    … …  
    572501    psS32 loop_t = 0;
    573502    psS32 i = 0;
    574     psF32 polySum = 0.0;
     503    psF64 polySum = 0.0;
    575504    psPolynomial1D* *chebPolys = NULL;
    576505    psS32 maxChebyPoly = 0;
    … …  
    612541    return(polySum);
    613542}
    614 
    615 /*****************************************************************************
    616     Polynomial coefficients will be accessed in [w][x][y][z] fashion.
    617  *****************************************************************************/
    618 static psF64 dOrdPolynomial1DEval(psF64 x, const psDPolynomial1D* poly)
    619 {
    620     psS32 loop_x = 0;
    621     psF64 polySum = 0.0;
    622     psF64 xSum = 1.0;
    623 
    624     for (loop_x = 0; loop_x < poly->n; loop_x++) {
    625         if (poly->mask[loop_x] == 0) {
    626             polySum += xSum * poly->coeff[loop_x];
    627         }
    628         xSum *= x;
    629     }
    630 
    631     return(polySum);
    632 }
    633 
    634 // XXX: You can do this without having to psAlloc() vector d.
    635 // XXX: How does the mask vector effect Crenshaw's formula?
    636 static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* poly)
    637 {
    638     PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
    639     psVector *d;
    640     psS32 n;
    641     psS32 i;
    642     psF64 tmp;
    643 
    644     n = poly->n;
    645     d = psVectorAlloc(n, PS_TYPE_F64);
    646     if(poly->mask[n-1] == 0) {
    647         d->data.F64[n-1] = poly->coeff[n-1];
    648     } else {
    649         d->data.F64[n-1] = 0.0;
    650     }
    651     d->data.F64[n-2] = (2.0 * x * d->data.F64[n-1]);
    652     if(poly->mask[n-2] == 0) {
    653         d->data.F64[n-2] += poly->coeff[n-2];
    654     }
    655     for (i=n-3;i>=1;i--) {
    656         d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) -
    657                          (d->data.F64[i+2]);
    658         if(poly->mask[i] == 0) {
    659             d->data.F64[i] += poly->coeff[i];
    660         }
    661     }
    662 
    663     tmp = (x * d->data.F64[1]) -
    664           (d->data.F64[2]);
    665     if(poly->mask[0] == 0) {
    666         tmp += (0.5 * poly->coeff[0]);
    667     }
    668 
    669     psFree(d);
    670     return(tmp);
    671 }
    672 
    673 static psF64 dOrdPolynomial2DEval(psF64 x,
    674                                   psF64 y,
    675                                   const psDPolynomial2D* poly)
    676 {
    677     psS32 loop_x = 0;
    678     psS32 loop_y = 0;
    679     psF64 polySum = 0.0;
    680     psF64 xSum = 1.0;
    681     psF64 ySum = 1.0;
    682 
    683     for (loop_x = 0; loop_x < poly->nX; loop_x++) {
    684         ySum = xSum;
    685         for (loop_y = 0; loop_y < poly->nY; loop_y++) {
    686             if (poly->mask[loop_x][loop_y] == 0) {
    687                 polySum += ySum * poly->coeff[loop_x][loop_y];
    688             }
    689             ySum *= y;
    690         }
    691         xSum *= x;
    692     }
    693 
    694     return(polySum);
    695 }
    696 
    697 static psF64 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* poly)
    698 {
    699     PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
    700     PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0);
    701     psS32 loop_x = 0;
    702     psS32 loop_y = 0;
    703     psS32 i = 0;
    704     psF64 polySum = 0.0;
    705     psPolynomial1D* *chebPolys = NULL;
    706     psS32 maxChebyPoly = 0;
    707 
    708     // Determine how many Chebyshev polynomials
    709     // are needed, then create them.
    710     maxChebyPoly = poly->nX;
    711     if (poly->nY > maxChebyPoly) {
    712         maxChebyPoly = poly->nY;
    713     }
    714     chebPolys = createChebyshevPolys(maxChebyPoly);
    715 
    716     for (loop_x = 0; loop_x < poly->nX; loop_x++) {
    717         for (loop_y = 0; loop_y < poly->nY; loop_y++) {
    718             if (poly->mask[loop_x][loop_y] == 0) {
    719                 polySum += poly->coeff[loop_x][loop_y] *
    720                            psPolynomial1DEval(chebPolys[loop_x], x) *
    721                            psPolynomial1DEval(chebPolys[loop_y], y);
    722             }
    723         }
    724     }
    725 
    726     for (i=0;i<maxChebyPoly;i++) {
    727         psFree(chebPolys[i]);
    728     }
    729     psFree(chebPolys);
    730     return(polySum);
    731 }
    732 
    733 static psF64 dOrdPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* poly)
    734 {
    735     psS32 loop_x = 0;
    736     psS32 loop_y = 0;
    737     psS32 loop_z = 0;
    738     psF64 polySum = 0.0;
    739     psF64 xSum = 1.0;
    740     psF64 ySum = 1.0;
    741     psF64 zSum = 1.0;
    742 
    743     for (loop_x = 0; loop_x < poly->nX; loop_x++) {
    744         ySum = xSum;
    745         for (loop_y = 0; loop_y < poly->nY; loop_y++) {
    746             zSum = ySum;
    747             for (loop_z = 0; loop_z < poly->nZ; loop_z++) {
    748                 if (poly->mask[loop_x][loop_y][loop_z] == 0) {
    749                     polySum += zSum * poly->coeff[loop_x][loop_y][loop_z];
    750                 }
    751                 zSum *= z;
    752             }
    753             ySum *= y;
    754         }
    755         xSum *= x;
    756     }
    757 
    758     return(polySum);
    759 }
    760 
    761 static psF64 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* poly)
    762 {
    763     PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
    764     PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0);
    765     PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0);
    766     psS32 loop_x = 0;
    767     psS32 loop_y = 0;
    768     psS32 loop_z = 0;
    769     psS32 i = 0;
    770     psF64 polySum = 0.0;
    771     psPolynomial1D* *chebPolys = NULL;
    772     psS32 maxChebyPoly = 0;
    773 
    774     // Determine how many Chebyshev polynomials
    775     // are needed, then create them.
    776     maxChebyPoly = poly->nX;
    777     if (poly->nY > maxChebyPoly) {
    778         maxChebyPoly = poly->nY;
    779     }
    780     if (poly->nZ > maxChebyPoly) {
    781         maxChebyPoly = poly->nZ;
    782     }
    783     chebPolys = createChebyshevPolys(maxChebyPoly);
    784 
    785     for (loop_x = 0; loop_x < poly->nX; loop_x++) {
    786         for (loop_y = 0; loop_y < poly->nY; loop_y++) {
    787             for (loop_z = 0; loop_z < poly->nZ; loop_z++) {
    788                 if (poly->mask[loop_x][loop_y][loop_z] == 0) {
    789                     polySum += poly->coeff[loop_x][loop_y][loop_z] *
    790                                psPolynomial1DEval(chebPolys[loop_x], x) *
    791                                psPolynomial1DEval(chebPolys[loop_y], y) *
    792                                psPolynomial1DEval(chebPolys[loop_z], z);
    793                 }
    794             }
    795         }
    796     }
    797 
    798     for (i=0;i<maxChebyPoly;i++) {
    799         psFree(chebPolys[i]);
    800     }
    801     psFree(chebPolys);
    802     return(polySum);
    803 }
    804 
    805 static psF64 dOrdPolynomial4DEval(psF64 x, psF64 y, psF64 z, psF64 t, const psDPolynomial4D* poly)
    806 {
    807     psS32 loop_x = 0;
    808     psS32 loop_y = 0;
    809     psS32 loop_z = 0;
    810     psS32 loop_t = 0;
    811     psF64 polySum = 0.0;
    812     psF64 xSum = 1.0;
    813     psF64 ySum = 1.0;
    814     psF64 zSum = 1.0;
    815     psF64 tSum = 1.0;
    816 
    817     for (loop_x = 0; loop_x < poly->nX; loop_x++) {
    818         ySum = xSum;
    819         for (loop_y = 0; loop_y < poly->nY; loop_y++) {
    820             zSum = ySum;
    821             for (loop_z = 0; loop_z < poly->nZ; loop_z++) {
    822                 tSum = zSum;
    823                 for (loop_t = 0; loop_t < poly->nT; loop_t++) {
    824                     if (poly->mask[loop_x][loop_y][loop_z][loop_t] == 0) {
    825                         polySum += tSum * poly->coeff[loop_x][loop_y][loop_z][loop_t];
    826                     }
    827                     tSum *= t;
    828                 }
    829                 zSum *= z;
    830             }
    831             ySum *= y;
    832         }
    833         xSum *= x;
    834     }
    835 
    836     return(polySum);
    837 }
    838 
    839 static psF64 dChebPolynomial4DEval(psF64 x, psF64 y, psF64 z, psF64 t, const psDPolynomial4D* poly)
    840 {
    841     PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
    842     PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0);
    843     PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0);
    844     PS_ASSERT_FLOAT_WITHIN_RANGE(t, -1.0, 1.0, 0.0);
    845     psS32 loop_x = 0;
    846     psS32 loop_y = 0;
    847     psS32 loop_z = 0;
    848     psS32 loop_t = 0;
    849     psS32 i = 0;
    850     psF64 polySum = 0.0;
    851     psPolynomial1D* *chebPolys = NULL;
    852     psS32 maxChebyPoly = 0;
    853 
    854     // Determine how many Chebyshev polynomials
    855     // are needed, then create them.
    856     maxChebyPoly = poly->nX;
    857     if (poly->nY > maxChebyPoly) {
    858         maxChebyPoly = poly->nY;
    859     }
    860     if (poly->nZ > maxChebyPoly) {
    861         maxChebyPoly = poly->nZ;
    862     }
    863     if (poly->nT > maxChebyPoly) {
    864         maxChebyPoly = poly->nT;
    865     }
    866     chebPolys = createChebyshevPolys(maxChebyPoly);
    867 
    868     for (loop_x = 0; loop_x < poly->nX; loop_x++) {
    869         for (loop_y = 0; loop_y < poly->nY; loop_y++) {
    870             for (loop_z = 0; loop_z < poly->nZ; loop_z++) {
    871                 for (loop_t = 0; loop_t < poly->nT; loop_t++) {
    872                     if (poly->mask[loop_x][loop_y][loop_z][loop_t] == 0) {
    873                         polySum += poly->coeff[loop_x][loop_y][loop_z][loop_t] *
    874                                    psPolynomial1DEval(chebPolys[loop_x], x) *
    875                                    psPolynomial1DEval(chebPolys[loop_y], y) *
    876                                    psPolynomial1DEval(chebPolys[loop_z], z) *
    877                                    psPolynomial1DEval(chebPolys[loop_t], t);
    878                     }
    879                 }
    880             }
    881         }
    882     }
    883 
    884     for (i=0;i<maxChebyPoly;i++) {
    885         psFree(chebPolys[i]);
    886     }
    887     psFree(chebPolys);
    888     return(polySum);
    889 }
    890 
    891543
    892544/*****************************************************************************
    … …  
    1082734    newPoly->type = type;
    1083735    newPoly->n = n;
    1084     newPoly->coeff = (psF32 *)psAlloc(n * sizeof(psF32));
    1085     newPoly->coeffErr = (psF32 *)psAlloc(n * sizeof(psF32));
     736    newPoly->coeff = psAlloc(n * sizeof(psF64));
     737    newPoly->coeffErr = psAlloc(n * sizeof(psF64));
    1086738    newPoly->mask = (char *)psAlloc(n * sizeof(char));
    1087739    for (i = 0; i < n; i++) {
    … …  
    1111763    newPoly->nY = nY;
    1112764
    1113     newPoly->coeff = (psF32 **)psAlloc(nX * sizeof(psF32 *));
    1114     newPoly->coeffErr = (psF32 **)psAlloc(nX * sizeof(psF32 *));
     765    newPoly->coeff = psAlloc(nX * sizeof(psF64 *));
     766    newPoly->coeffErr = psAlloc(nX * sizeof(psF64 *));
    1115767    newPoly->mask = (char **)psAlloc(nX * sizeof(char *));
    1116768    for (x = 0; x < nX; x++) {
    1117         newPoly->coeff[x] = (psF32 *)psAlloc(nY * sizeof(psF32));
    1118         newPoly->coeffErr[x] = (psF32 *)psAlloc(nY * sizeof(psF32));
     769        newPoly->coeff[x] = psAlloc(nY * sizeof(psF64));
     770        newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64));
    1119771        newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char));
    1120772    }
    … …  
    1150802    newPoly->nZ = nZ;
    1151803
    1152     newPoly->coeff = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
    1153     newPoly->coeffErr = (psF32 ***)psAlloc(nX * sizeof(psF32 **));
     804    newPoly->coeff = psAlloc(nX * sizeof(psF64 **));
     805    newPoly->coeffErr = psAlloc(nX * sizeof(psF64 **));
    1154806    newPoly->mask = (char ***)psAlloc(nX * sizeof(char **));
    1155807    for (x = 0; x < nX; x++) {
    1156         newPoly->coeff[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
    1157         newPoly->coeffErr[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *));
     808        newPoly->coeff[x] = psAlloc(nY * sizeof(psF64 *));
     809        newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64 *));
    1158810        newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *));
    1159811        for (y = 0; y < nY; y++) {
    1160             newPoly->coeff[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
    1161             newPoly->coeffErr[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32));
     812            newPoly->coeff[x][y] = psAlloc(nZ * sizeof(psF64));
     813            newPoly->coeffErr[x][y] = psAlloc(nZ * sizeof(psF64));
    1162814            newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char));
    1163815        }
    … …  
    1199851    newPoly->nT = nT;
    1200852
    1201     newPoly->coeff = (psF32 ****)psAlloc(nX * sizeof(psF32 ***));
    1202     newPoly->coeffErr = (psF32 ****)psAlloc(nX * sizeof(psF32 ***));
     853    newPoly->coeff = psAlloc(nX * sizeof(psF64 ***));
     854    newPoly->coeffErr = psAlloc(nX * sizeof(psF64 ***));
    1203855    newPoly->mask = (char ****)psAlloc(nX * sizeof(char ***));
    1204856    for (x = 0; x < nX; x++) {
    1205         newPoly->coeff[x] = (psF32 ***)psAlloc(nY * sizeof(psF32 **));
    1206         newPoly->coeffErr[x] = (psF32 ***)psAlloc(nY * sizeof(psF32 **));
     857        newPoly->coeff[x] = psAlloc(nY * sizeof(psF64 **));
     858        newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64 **));
    1207859        newPoly->mask[x] = (char ***)psAlloc(nY * sizeof(char **));
    1208860        for (y = 0; y < nY; y++) {
    1209             newPoly->coeff[x][y] = (psF32 **)psAlloc(nZ * sizeof(psF32 *));
    1210             newPoly->coeffErr[x][y] = (psF32 **)psAlloc(nZ * sizeof(psF32 *));
     861            newPoly->coeff[x][y] = psAlloc(nZ * sizeof(psF64 *));
     862            newPoly->coeffErr[x][y] = psAlloc(nZ * sizeof(psF64 *));
    1211863            newPoly->mask[x][y] = (char **)psAlloc(nZ * sizeof(char *));
    1212864            for (z = 0; z < nZ; z++) {
    1213                 newPoly->coeff[x][y][z] = (psF32 *)psAlloc(nT * sizeof(psF32));
    1214                 newPoly->coeffErr[x][y][z] = (psF32 *)psAlloc(nT * sizeof(psF32));
     865                newPoly->coeff[x][y][z] = psAlloc(nT * sizeof(psF64));
     866                newPoly->coeffErr[x][y][z] = psAlloc(nT * sizeof(psF64));
    1215867                newPoly->mask[x][y][z] = (char *)psAlloc(nT * sizeof(char));
    1216868            }
    … …  
    1253905    PS_ASSERT_POLY_NON_NULL(poly, NULL);
    1254906    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
    1255     PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
     907    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
    1256908
    1257909    psVector *tmp;
    1258910
    1259     tmp = psVectorAlloc(x->n, PS_TYPE_F32);
     911    tmp = psVectorAlloc(x->n, PS_TYPE_F64);
    1260912    for (psS32 i=0;i<x->n;i++) {
    1261         tmp->data.F32[i] = psPolynomial1DEval(poly, x->data.F32[i]);
     913        tmp->data.F64[i] = psPolynomial1DEval(poly, x->data.F64[i]);
    1262914    }
    1263915
    … …  
    1288940    PS_ASSERT_POLY_NON_NULL(poly, NULL);
    1289941    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
    1290     PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
     942    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
    1291943    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
    1292     PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
     944    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL);
    1293945
    1294946    psVector *tmp;
    … …  
    1301953
    1302954    // Create output vector to return
    1303     tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
     955    tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
    1304956
    1305957    // Evaluate the polynomial at the specified points
    1306958    for (psS32 i=0; i<vecLen; i++) {
    1307         tmp->data.F32[i] = psPolynomial2DEval(poly,x->data.F32[i],y->data.F32[i]);
     959        tmp->data.F64[i] = psPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]);
    1308960    }
    1309961
    … …  
    1336988    PS_ASSERT_POLY_NON_NULL(poly, NULL);
    1337989    PS_ASSERT_VECTOR_NON_NULL(x, NULL);
    1338     PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
     990    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
    1339991    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
    1340     PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
     992    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL);
    1341993    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
    1342     PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL);
     994    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL);
    1343995
    1344996    psVector *tmp;
    … …  
    13541006
    13551007    // Allocate output vector
    1356     tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
     1008    tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
    13571009
    13581010    // Evaluate polynomial
    13591011    for (psS32 i = 0; i < vecLen; i++) {
    1360         tmp->data.F32[i] = psPolynomial3DEval(poly,
    1361                                               x->data.F32[i],
    1362                                               y->data.F32[i],
    1363                                               z->data.F32[i]);
     1012        tmp->data.F64[i] = psPolynomial3DEval(poly,
     1013                                              x->data.F64[i],
     1014                                              y->data.F64[i],
     1015                                              z->data.F64[i]);
    13641016    }
    13651017
    … …  
    13891041                                   const psVector *z,
    13901042                                   const psVector *t)
    1391 {
    1392     PS_ASSERT_POLY_NON_NULL(poly, NULL);
    1393     PS_ASSERT_VECTOR_NON_NULL(x, NULL);
    1394     PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
    1395     PS_ASSERT_VECTOR_NON_NULL(y, NULL);
    1396     PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
    1397     PS_ASSERT_VECTOR_NON_NULL(z, NULL);
    1398     PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL);
    1399     PS_ASSERT_VECTOR_NON_NULL(t, NULL);
    1400     PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F32, NULL);
    1401 
    1402     psVector *tmp;
    1403     psS32 vecLen=x->n;
    1404 
    1405     // Determine output vector size from min of input vectors
    1406     if (z->n < vecLen) {
    1407         vecLen = z->n;
    1408     }
    1409     if (y->n < vecLen) {
    1410         vecLen = y->n;
    1411     }
    1412     if (t->n < vecLen) {
    1413         vecLen = t->n;
    1414     }
    1415 
    1416     // Allocate output vector
    1417     tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
    1418 
    1419     // Evaluate polynomial
    1420     for (psS32 i = 0; i < vecLen; i++) {
    1421         tmp->data.F32[i] = psPolynomial4DEval(poly,
    1422                                               x->data.F32[i],
    1423                                               y->data.F32[i],
    1424                                               z->data.F32[i],
    1425                                               t->data.F32[i]);
    1426     }
    1427 
    1428     // Return output vector
    1429     return(tmp);
    1430 }
    1431 
    1432 
    1433 psDPolynomial1D* psDPolynomial1DAlloc( int n,
    1434                                        psPolynomialType type)
    1435 {
    1436     PS_ASSERT_INT_POSITIVE(n, NULL);
    1437 
    1438     unsigned int i = 0;
    1439     psDPolynomial1D* newPoly = NULL;
    1440 
    1441     newPoly = (psDPolynomial1D* ) psAlloc(sizeof(psDPolynomial1D));
    1442     psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial1DFree);
    1443 
    1444     newPoly->type = type;
    1445     newPoly->n = n;
    1446     newPoly->coeff = (psF64 *)psAlloc(n * sizeof(psF64));
    1447     newPoly->coeffErr = (psF64 *)psAlloc(n * sizeof(psF64));
    1448     newPoly->mask = (char *)psAlloc(n * sizeof(char));
    1449     for (i = 0; i < n; i++) {
    1450         newPoly->coeff[i] = 0.0;
    1451         newPoly->coeffErr[i] = 0.0;
    1452         newPoly->mask[i] = 0;
    1453     }
    1454 
    1455     return(newPoly);
    1456 }
    1457 
    1458 psDPolynomial2D* psDPolynomial2DAlloc( int nX,  int nY,
    1459                                        psPolynomialType type)
    1460 {
    1461     PS_ASSERT_INT_POSITIVE(nX, NULL);
    1462     PS_ASSERT_INT_POSITIVE(nY, NULL);
    1463 
    1464     unsigned int x = 0;
    1465     unsigned int y = 0;
    1466     psDPolynomial2D* newPoly = NULL;
    1467 
    1468     newPoly = (psDPolynomial2D* ) psAlloc(sizeof(psDPolynomial2D));
    1469     psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial2DFree);
    1470 
    1471     newPoly->type = type;
    1472     newPoly->nX = nX;
    1473     newPoly->nY = nY;
    1474 
    1475     newPoly->coeff = (psF64 **)psAlloc(nX * sizeof(psF64 *));
    1476     newPoly->coeffErr = (psF64 **)psAlloc(nX * sizeof(psF64 *));
    1477     newPoly->mask = (char **)psAlloc(nX * sizeof(char *));
    1478     for (x = 0; x < nX; x++) {
    1479         newPoly->coeff[x] = (psF64 *)psAlloc(nY * sizeof(psF64));
    1480         newPoly->coeffErr[x] = (psF64 *)psAlloc(nY * sizeof(psF64));
    1481         newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char));
    1482     }
    1483     for (x = 0; x < nX; x++) {
    1484         for (y = 0; y < nY; y++) {
    1485             newPoly->coeff[x][y] = 0.0;
    1486             newPoly->coeffErr[x][y] = 0.0;
    1487             newPoly->mask[x][y] = 0;
    1488         }
    1489     }
    1490 
    1491     return(newPoly);
    1492 }
    1493 
    1494 psDPolynomial3D* psDPolynomial3DAlloc( int nX,  int nY,  int nZ,
    1495                                        psPolynomialType type)
    1496 {
    1497     PS_ASSERT_INT_POSITIVE(nX, NULL);
    1498     PS_ASSERT_INT_POSITIVE(nY, NULL);
    1499     PS_ASSERT_INT_POSITIVE(nZ, NULL);
    1500 
    1501     unsigned int x = 0;
    1502     unsigned int y = 0;
    1503     unsigned int z = 0;
    1504     psDPolynomial3D* newPoly = NULL;
    1505 
    1506     newPoly = (psDPolynomial3D* ) psAlloc(sizeof(psDPolynomial3D));
    1507     psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial3DFree);
    1508 
    1509     newPoly->type = type;
    1510     newPoly->nX = nX;
    1511     newPoly->nY = nY;
    1512     newPoly->nZ = nZ;
    1513 
    1514     newPoly->coeff = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
    1515     newPoly->coeffErr = (psF64 ***)psAlloc(nX * sizeof(psF64 **));
    1516     newPoly->mask = (char ***)psAlloc(nX * sizeof(char **));
    1517     for (x = 0; x < nX; x++) {
    1518         newPoly->coeff[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
    1519         newPoly->coeffErr[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *));
    1520         newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *));
    1521         for (y = 0; y < nY; y++) {
    1522             newPoly->coeff[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
    1523             newPoly->coeffErr[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64));
    1524             newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char));
    1525         }
    1526     }
    1527     for (x = 0; x < nX; x++) {
    1528         for (y = 0; y < nY; y++) {
    1529             for (z = 0; z < nZ; z++) {
    1530                 newPoly->coeff[x][y][z] = 0.0;
    1531                 newPoly->coeffErr[x][y][z] = 0.0;
    1532                 newPoly->mask[x][y][z] = 0;
    1533             }
    1534         }
    1535     }
    1536 
    1537     return(newPoly);
    1538 }
    1539 
    1540 psDPolynomial4D* psDPolynomial4DAlloc( int nX,  int nY,  int nZ,  int nT,
    1541                                        psPolynomialType type)
    1542 {
    1543     PS_ASSERT_INT_POSITIVE(nX, NULL);
    1544     PS_ASSERT_INT_POSITIVE(nY, NULL);
    1545     PS_ASSERT_INT_POSITIVE(nZ, NULL);
    1546     PS_ASSERT_INT_POSITIVE(nT, NULL);
    1547 
    1548     unsigned int x = 0;
    1549     unsigned int y = 0;
    1550     unsigned int z = 0;
    1551     unsigned int t = 0;
    1552     psDPolynomial4D* newPoly = NULL;
    1553 
    1554     newPoly = (psDPolynomial4D* ) psAlloc(sizeof(psDPolynomial4D));
    1555     psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial4DFree);
    1556 
    1557     newPoly->type = type;
    1558     newPoly->nX = nX;
    1559     newPoly->nY = nY;
    1560     newPoly->nZ = nZ;
    1561     newPoly->nT = nT;
    1562 
    1563     newPoly->coeff = (psF64 ****)psAlloc(nX * sizeof(psF64 ***));
    1564     newPoly->coeffErr = (psF64 ****)psAlloc(nX * sizeof(psF64 ***));
    1565     newPoly->mask = (char ****)psAlloc(nX * sizeof(char ***));
    1566     for (x = 0; x < nX; x++) {
    1567         newPoly->coeff[x] = (psF64 ***)psAlloc(nY * sizeof(psF64 **));
    1568         newPoly->coeffErr[x] = (psF64 ***)psAlloc(nY * sizeof(psF64 **));
    1569         newPoly->mask[x] = (char ***)psAlloc(nY * sizeof(char **));
    1570         for (y = 0; y < nY; y++) {
    1571             newPoly->coeff[x][y] = (psF64 **)psAlloc(nZ * sizeof(psF64 *));
    1572             newPoly->coeffErr[x][y] = (psF64 **)psAlloc(nZ * sizeof(psF64 *));
    1573             newPoly->mask[x][y] = (char **)psAlloc(nZ * sizeof(char *));
    1574             for (z = 0; z < nZ; z++) {
    1575                 newPoly->coeff[x][y][z] = (psF64 *)psAlloc(nT * sizeof(psF64));
    1576                 newPoly->coeffErr[x][y][z] = (psF64 *)psAlloc(nT * sizeof(psF64));
    1577                 newPoly->mask[x][y][z] = (char *)psAlloc(nT * sizeof(char));
    1578             }
    1579         }
    1580     }
    1581     for (x = 0; x < nX; x++) {
    1582         for (y = 0; y < nY; y++) {
    1583             for (z = 0; z < nZ; z++) {
    1584                 for (t = 0; t < nT; t++) {
    1585                     newPoly->coeff[x][y][z][t] = 0.0;
    1586                     newPoly->coeffErr[x][y][z][t] = 0.0;
    1587                     newPoly->mask[x][y][z][t] = 0;
    1588                 }
    1589             }
    1590         }
    1591     }
    1592 
    1593     return(newPoly);
    1594 }
    1595 
    1596 
    1597 psF64 psDPolynomial1DEval(const psDPolynomial1D* poly, psF64 x)
    1598 {
    1599     PS_ASSERT_POLY_NON_NULL(poly, NAN);
    1600 
    1601     if (poly->type == PS_POLYNOMIAL_ORD) {
    1602         return(dOrdPolynomial1DEval(x, poly));
    1603     } else if (poly->type == PS_POLYNOMIAL_CHEB) {
    1604         return(dChebPolynomial1DEval(x, poly));
    1605     } else {
    1606         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    1607                 PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE,
    1608                 poly->type);
    1609     }
    1610     return(NAN);
    1611 }
    1612 
    1613 psVector *psDPolynomial1DEvalVector(const psDPolynomial1D *poly,
    1614                                     const psVector *x)
    1615 
    1616 {
    1617     PS_ASSERT_POLY_NON_NULL(poly, NULL);
    1618     PS_ASSERT_VECTOR_NON_NULL(x, NULL);
    1619     PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
    1620 
    1621     psVector *tmp;
    1622 
    1623     tmp = psVectorAlloc(x->n, PS_TYPE_F64);
    1624     for (psS32 i=0;i<x->n;i++) {
    1625         tmp->data.F64[i] = psDPolynomial1DEval(poly,
    1626                                                x->data.F64[i]);
    1627     }
    1628 
    1629     return(tmp);
    1630 }
    1631 
    1632 
    1633 psF64 psDPolynomial2DEval(const psDPolynomial2D* poly,
    1634                           psF64 x,
    1635                           psF64 y)
    1636 {
    1637     PS_ASSERT_POLY_NON_NULL(poly, NAN);
    1638     if (poly->type == PS_POLYNOMIAL_ORD) {
    1639         return(dOrdPolynomial2DEval(x, y, poly));
    1640     } else if (poly->type == PS_POLYNOMIAL_CHEB) {
    1641         return(dChebPolynomial2DEval(x, y, poly));
    1642     } else {
    1643         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    1644                 PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE,
    1645                 poly->type);
    1646     }
    1647     return(NAN);
    1648 }
    1649 
    1650 psVector *psDPolynomial2DEvalVector(const psDPolynomial2D *poly,
    1651                                     const psVector *x,
    1652                                     const psVector *y)
    1653 {
    1654     PS_ASSERT_POLY_NON_NULL(poly, NULL);
    1655     PS_ASSERT_VECTOR_NON_NULL(x, NULL);
    1656     PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
    1657     PS_ASSERT_VECTOR_NON_NULL(y, NULL);
    1658     PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL);
    1659 
    1660     psVector *tmp;
    1661     psS32 vecLen=x->n;
    1662 
    1663     // Determine the output vector length from minimum length of input vectors
    1664     if (y->n < vecLen) {
    1665         vecLen = y->n;
    1666     }
    1667 
    1668     // Allocate output vector
    1669     tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
    1670 
    1671     // Evaluate the polynomial
    1672     for (psS32 i = 0; i < vecLen; i++) {
    1673         tmp->data.F64[i] = psDPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]);
    1674     }
    1675 
    1676     // Return output vector
    1677     return(tmp);
    1678 }
    1679 
    1680 
    1681 psF64 psDPolynomial3DEval(const psDPolynomial3D* poly,
    1682                           psF64 x,
    1683                           psF64 y,
    1684                           psF64 z)
    1685 {
    1686     PS_ASSERT_POLY_NON_NULL(poly, NAN);
    1687 
    1688     if (poly->type == PS_POLYNOMIAL_ORD) {
    1689         return(dOrdPolynomial3DEval(x, y, z, poly));
    1690     } else if (poly->type == PS_POLYNOMIAL_CHEB) {
    1691         return(dChebPolynomial3DEval(x, y, z, poly));
    1692     } else {
    1693         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    1694                 PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE,
    1695                 poly->type);
    1696     }
    1697     return(NAN);
    1698 }
    1699 
    1700 psVector *psDPolynomial3DEvalVector(const psDPolynomial3D *poly,
    1701                                     const psVector *x,
    1702                                     const psVector *y,
    1703                                     const psVector *z)
    1704 
    1705 {
    1706     PS_ASSERT_POLY_NON_NULL(poly, NULL);
    1707     PS_ASSERT_VECTOR_NON_NULL(x, NULL);
    1708     PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
    1709     PS_ASSERT_VECTOR_NON_NULL(y, NULL);
    1710     PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL);
    1711     PS_ASSERT_VECTOR_NON_NULL(z, NULL);
    1712     PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL);
    1713 
    1714     psVector *tmp;
    1715     psS32 vecLen=x->n;
    1716 
    1717     // Determine the size of output vector from min of input vectors
    1718     if (y->n < vecLen) {
    1719         vecLen = y->n;
    1720     }
    1721     if (z->n < vecLen) {
    1722         vecLen = z->n;
    1723     }
    1724 
    1725     // Allocate output vector
    1726     tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
    1727 
    1728     // Evaluate polynomial
    1729     for (psS32 i = 0; i < vecLen; i++) {
    1730         tmp->data.F64[i] = psDPolynomial3DEval(poly,
    1731                                                x->data.F64[i],
    1732                                                y->data.F64[i],
    1733                                                z->data.F64[i]);
    1734     }
    1735 
    1736     // Return output vector
    1737     return(tmp);
    1738 }
    1739 
    1740 psF64 psDPolynomial4DEval(const psDPolynomial4D* poly,
    1741                           psF64 x,
    1742                           psF64 y,
    1743                           psF64 z,
    1744                           psF64 t)
    1745 {
    1746     PS_ASSERT_POLY_NON_NULL(poly, NAN);
    1747 
    1748     if (poly->type == PS_POLYNOMIAL_ORD) {
    1749         return(dOrdPolynomial4DEval(x,y,z,t, poly));
    1750     } else if (poly->type == PS_POLYNOMIAL_CHEB) {
    1751         return(dChebPolynomial4DEval(x,y,z,t, poly));
    1752     } else {
    1753         psError(PS_ERR_BAD_PARAMETER_TYPE, true,
    1754                 PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE,
    1755                 poly->type);
    1756     }
    1757     return(NAN);
    1758 }
    1759 
    1760 psVector *psDPolynomial4DEvalVector(const psDPolynomial4D *poly,
    1761                                     const psVector *x,
    1762                                     const psVector *y,
    1763                                     const psVector *z,
    1764                                     const psVector *t)
    17651043{
    17661044    PS_ASSERT_POLY_NON_NULL(poly, NULL);
    … …  
    17771055    psS32 vecLen=x->n;
    17781056
    1779     // Determine the output vector size from min of input vectors
     1057    // Determine output vector size from min of input vectors
    17801058    if (z->n < vecLen) {
    17811059        vecLen = z->n;
    … …  
    17911069    tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
    17921070
    1793     // Evaluate the polynomial
     1071    // Evaluate polynomial
    17941072    for (psS32 i = 0; i < vecLen; i++) {
    1795         tmp->data.F64[i] = psDPolynomial4DEval(poly,
    1796                                                x->data.F64[i],
    1797                                                y->data.F64[i],
    1798                                                z->data.F64[i],
    1799                                                t->data.F64[i]);
     1073        tmp->data.F64[i] = psPolynomial4DEval(poly,
     1074                                              x->data.F64[i],
     1075                                              y->data.F64[i],
     1076                                              z->data.F64[i],
     1077                                              t->data.F64[i]);
    18001078    }
    18011079
    … …  
    18031081    return(tmp);
    18041082}
    1805 
    1806 
    1807 
    18081083
    18091084//typedef struct {
Note: See TracChangeset for help on using the changeset viewer.