IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 5090 for trunk/psLib/src/math


Ignore:
Timestamp:
Sep 21, 2005, 4:47:16 PM (21 years ago)
Author:
gusciora
Message:

The nX terms for the polynomials now refer to the order, not the number
of terms.

Location:
trunk/psLib/src/math
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r5041 r5090  
    1010 *  @author EAM, IfA
    1111 *
    12  *  @version $Revision: 1.138 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2005-09-13 21:43:47 $
     12 *  @version $Revision: 1.139 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2005-09-22 02:47:16 $
    1414 *
    1515 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    14591459    // XXX: these ASSERTS are redundant.
    14601460    PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
    1461     PS_ASSERT_INT_NONNEGATIVE(myPoly->n, NULL);
     1461    PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_1D_n, NULL);
    14621462    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
    14631463    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL);
    … …  
    15231523    fac = 2.0/((psF32) n);
    15241524    // XXX: is this loop bound correct?
    1525     for (j=0;j<myPoly->n;j++) {
     1525    for (j=0;j<myPoly->COOL_1D_n+1;j++) {
    15261526        sum = 0.0;
    15271527        for (k=0;k<n;k++) {
    … …  
    15521552    // XXX: these ASSERTS are redundant.
    15531553    PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
    1554     PS_ASSERT_INT_NONNEGATIVE(myPoly->n, NULL);
     1554    PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_1D_n, NULL);
    15551555    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
    15561556    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL);
    … …  
    15971597    }
    15981598
    1599     nTerm = myPoly->n;
     1599    nTerm = 1 + myPoly->COOL_1D_n;
    16001600    nOrder = nTerm - 1;
    16011601
    … …  
    16971697
    16981698    PS_ASSERT_POLY_NON_NULL(poly, NULL);
    1699     PS_ASSERT_INT_NONNEGATIVE(poly->n, NULL);
     1699    PS_ASSERT_INT_NONNEGATIVE(poly->COOL_1D_n, NULL);
    17001700    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
    17011701    PS_ASSERT_VECTOR_NON_EMPTY(f, NULL);
    … …  
    18371837    // These ASSERTS are redundant.
    18381838    PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
    1839     PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
    1840     PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL);
     1839    PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_2D_nX, NULL);
     1840    PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_2D_nY, NULL);
    18411841
    18421842    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
    … …  
    18651865
    18661866    // XXX:Watch for changes to the psPolys: nTerm != nOrder.
    1867     psS32 nXterm = myPoly->nX;
    1868     psS32 nYterm = myPoly->nY;
     1867    psS32 nXterm = 1 + myPoly->COOL_2D_nX;
     1868    psS32 nYterm = 1 + myPoly->COOL_2D_nY;
    18691869    nTerm = nXterm * nYterm;
    18701870
    … …  
    22642264    // These ASSERTS are redundant.
    22652265    PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
    2266     PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
    2267     PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL);
    2268     PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, NULL);
     2266    PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_3D_nX, NULL);
     2267    PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_3D_nY, NULL);
     2268    PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_3D_nZ, NULL);
    22692269
    22702270    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
    … …  
    25272527    // These ASSERTS are redundant.
    25282528    PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
    2529     PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
    2530     PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL);
    2531     PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, NULL);
    2532     PS_ASSERT_INT_NONNEGATIVE(myPoly->nT, NULL);
     2529    PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_4D_nX, NULL);
     2530    PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_4D_nY, NULL);
     2531    PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_4D_nZ, NULL);
     2532    PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_4D_nT, NULL);
    25332533    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
    25342534    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL);

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

    r5074 r5090  
    77*  polynomials.  It also contains a Gaussian functions.
    88*
    9 *  @version $Revision: 1.122 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2005-09-20 20:48:35 $
     9*  @version $Revision: 1.123 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2005-09-22 02:47:16 $
    1111*
    1212*  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    9595    unsigned int x = 0;
    9696
    97     for (x = 0; x < poly->nX; x++) {
     97    for (x = 0; x < poly->COOL_2D_nX; x++) {
    9898        psFree(poly->coeff[x]);
    9999        psFree(poly->coeffErr[x]);
    … …  
    110110    unsigned int y = 0;
    111111
    112     for (x = 0; x < poly->nX; x++) {
    113         for (y = 0; y < poly->nY; y++) {
     112    for (x = 0; x < poly->COOL_3D_nX; x++) {
     113        for (y = 0; y < poly->COOL_3D_nY; y++) {
    114114            psFree(poly->coeff[x][y]);
    115115            psFree(poly->coeffErr[x][y]);
    … …  
    132132    unsigned int z = 0;
    133133
    134     for (x = 0; x < poly->nX; x++) {
    135         for (y = 0; y < poly->nY; y++) {
    136             for (z = 0; z < poly->nZ; z++) {
     134    for (x = 0; x < poly->COOL_4D_nX; x++) {
     135        for (y = 0; y < poly->COOL_4D_nY; y++) {
     136            for (z = 0; z < poly->COOL_4D_nZ; z++) {
    137137                psFree(poly->coeff[x][y][z]);
    138138                psFree(poly->coeffErr[x][y][z]);
    … …  
    170170    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *));
    171171    for (unsigned int i = 0; i < maxChebyPoly; i++) {
    172         chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD);
     172        chebPolys[i] = psPolynomial1DAlloc(i, PS_POLYNOMIAL_ORD);
    173173    }
    174174
    … …  
    182182
    183183        for (psS32 i = 2; i < maxChebyPoly; i++) {
    184             for (psS32 j = 0; j < chebPolys[i - 1]->n; j++) {
     184            for (psS32 j = 0; j < chebPolys[i - 1]->COOL_1D_n; j++) {
    185185                chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j];
    186186            }
    187             for (psS32 j = 0; j < chebPolys[i - 2]->n; j++) {
     187            for (psS32 j = 0; j < chebPolys[i - 2]->COOL_1D_n; j++) {
    188188                chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j];
    189189            }
    … …  
    210210            "---- Calling ordPolynomial1DEval(%lf)\n", x);
    211211    psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 4,
    212             "Polynomial order is %u\n", poly->n);
    213     for (loop_x = 0; loop_x < poly->n; loop_x++) {
     212            "Polynomial order is %u\n", poly->COOL_1D_n);
     213    for (loop_x = 0; loop_x < poly->COOL_1D_n+1; loop_x++) {
    214214        psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 4,
    215215                "Polynomial coeff[%u] is %lf\n", loop_x, poly->coeff[loop_x]);
    216216    }
    217217
    218     for (loop_x = 0; loop_x < poly->n; loop_x++) {
     218    for (loop_x = 0; loop_x < poly->COOL_1D_n+1; loop_x++) {
    219219        if (poly->mask[loop_x] == 0) {
    220220            // XXX: If you set the tracelevel to 10 here, and later set the tracelevel to
    … …  
    238238    PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
    239239    // XXX: Create a macro for this in psConstants.h
    240     if (poly->n < 1) {
    241         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error: Chebyshev polynomial as order %u.", poly->n);
     240    if (poly->COOL_1D_n < 1) {
     241        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error: Chebyshev polynomial is order %u.", poly->COOL_1D_n);
    242242        return(NAN);
    243243    }
    244244    psVector *d;
    245     unsigned int n = poly->n;
     245    // XXX: n should be nTerms here (for clarity).  Or get rid of the variable.
     246    unsigned int n = 1 + poly->COOL_1D_n;
    246247    unsigned int i;
    247248    psF64 tmp = 0.0;
    … …  
    302303    psPolynomial1D **chebPolys = NULL;
    303304
    304     n = poly->n;
     305    n = 1 + poly->COOL_1D_n;
    305306    chebPolys = createChebyshevPolys(n);
    306307
    307308    tmp = 0.0;
    308     for (i=0;i<poly->n;i++) {
     309    for (i=0;i<(1 + poly->COOL_1D_n);i++) {
    309310        tmp+= (poly->coeff[i] * psPolynomial1DEval(x, chebPolys[i]));
    310311    }
    … …  
    328329    psF64 ySum = 1.0;
    329330
    330     for (loop_x = 0; loop_x < poly->nX; loop_x++) {
     331    for (loop_x = 0; loop_x < (1 + poly->COOL_2D_nX); loop_x++) {
    331332        ySum = xSum;
    332         for (loop_y = 0; loop_y < poly->nY; loop_y++) {
     333        for (loop_y = 0; loop_y < (1 + poly->COOL_2D_nY); loop_y++) {
    333334            if (poly->mask[loop_x][loop_y] == 0) {
    334335                polySum += ySum * poly->coeff[loop_x][loop_y];
    … …  
    359360    // Determine how many Chebyshev polynomials
    360361    // are needed, then create them.
    361     maxChebyPoly = poly->nX;
    362     if (poly->nY > maxChebyPoly) {
    363         maxChebyPoly = poly->nY;
     362    maxChebyPoly = 1 + poly->COOL_2D_nX;
     363    if ((1 + poly->COOL_2D_nY) > maxChebyPoly) {
     364        maxChebyPoly = (1 + poly->COOL_2D_nY);
    364365    }
    365366    chebPolys = createChebyshevPolys(maxChebyPoly);
    366367
    367     for (loop_x = 0; loop_x < poly->nX; loop_x++) {
    368         for (loop_y = 0; loop_y < poly->nY; loop_y++) {
     368    for (loop_x = 0; loop_x < (1 + poly->COOL_2D_nX); loop_x++) {
     369        for (loop_y = 0; loop_y < (1 + poly->COOL_2D_nY); loop_y++) {
    369370            if (poly->mask[loop_x][loop_y] == 0) {
    370371                polySum += poly->coeff[loop_x][loop_y] *
    … …  
    394395    psF64 zSum = 1.0;
    395396
    396     for (loop_x = 0; loop_x < poly->nX; loop_x++) {
     397    for (loop_x = 0; loop_x < (1 + poly->COOL_3D_nX); loop_x++) {
    397398        ySum = xSum;
    398         for (loop_y = 0; loop_y < poly->nY; loop_y++) {
     399        for (loop_y = 0; loop_y < (1 + poly->COOL_3D_nY); loop_y++) {
    399400            zSum = ySum;
    400             for (loop_z = 0; loop_z < poly->nZ; loop_z++) {
     401            for (loop_z = 0; loop_z < (1 + poly->COOL_3D_nZ); loop_z++) {
    401402                if (poly->mask[loop_x][loop_y][loop_z] == 0) {
    402403                    polySum += zSum * poly->coeff[loop_x][loop_y][loop_z];
    … …  
    430431    // Determine how many Chebyshev polynomials
    431432    // are needed, then create them.
    432     maxChebyPoly = poly->nX;
    433     if (poly->nY > maxChebyPoly) {
    434         maxChebyPoly = poly->nY;
    435     }
    436     if (poly->nZ > maxChebyPoly) {
    437         maxChebyPoly = poly->nZ;
    438     }
    439     chebPolys = createChebyshevPolys(maxChebyPoly);
    440 
    441     for (loop_x = 0; loop_x < poly->nX; loop_x++) {
    442         for (loop_y = 0; loop_y < poly->nY; loop_y++) {
    443             for (loop_z = 0; loop_z < poly->nZ; loop_z++) {
     433    maxChebyPoly = poly->COOL_3D_nX;
     434    if (poly->COOL_3D_nY > maxChebyPoly) {
     435        maxChebyPoly = poly->COOL_3D_nY;
     436    }
     437    if (poly->COOL_3D_nZ > maxChebyPoly) {
     438        maxChebyPoly = poly->COOL_3D_nZ;
     439    }
     440    chebPolys = createChebyshevPolys(maxChebyPoly + 1);
     441
     442    for (loop_x = 0; loop_x < (1 + poly->COOL_3D_nX); loop_x++) {
     443        for (loop_y = 0; loop_y < (1 + poly->COOL_3D_nY); loop_y++) {
     444            for (loop_z = 0; loop_z < (1 + poly->COOL_3D_nZ); loop_z++) {
    444445                if (poly->mask[loop_x][loop_y][loop_z] == 0) {
    445446                    polySum += poly->coeff[loop_x][loop_y][loop_z] *
    … …  
    475476    psF64 tSum = 1.0;
    476477
    477     for (loop_x = 0; loop_x < poly->nX; loop_x++) {
     478    for (loop_x = 0; loop_x < (1 + poly->COOL_4D_nX); loop_x++) {
    478479        ySum = xSum;
    479         for (loop_y = 0; loop_y < poly->nY; loop_y++) {
     480        for (loop_y = 0; loop_y < (1 + poly->COOL_4D_nY); loop_y++) {
    480481            zSum = ySum;
    481             for (loop_z = 0; loop_z < poly->nZ; loop_z++) {
     482            for (loop_z = 0; loop_z < (1 + poly->COOL_4D_nZ); loop_z++) {
    482483                tSum = zSum;
    483                 for (loop_t = 0; loop_t < poly->nT; loop_t++) {
     484                for (loop_t = 0; loop_t < (1 + poly->COOL_4D_nT); loop_t++) {
    484485                    if (poly->mask[loop_x][loop_y][loop_z][loop_t] == 0) {
    485486                        polySum += tSum * poly->coeff[loop_x][loop_y][loop_z][loop_t];
    … …  
    518519    // Determine how many Chebyshev polynomials
    519520    // are needed, then create them.
    520     maxChebyPoly = poly->nX;
    521     if (poly->nY > maxChebyPoly) {
    522         maxChebyPoly = poly->nY;
    523     }
    524     if (poly->nZ > maxChebyPoly) {
    525         maxChebyPoly = poly->nZ;
    526     }
    527     if (poly->nT > maxChebyPoly) {
    528         maxChebyPoly = poly->nT;
     521    maxChebyPoly = poly->COOL_4D_nX;
     522    if (poly->COOL_4D_nY > maxChebyPoly) {
     523        maxChebyPoly = poly->COOL_4D_nY;
     524    }
     525    if (poly->COOL_4D_nZ > maxChebyPoly) {
     526        maxChebyPoly = poly->COOL_4D_nZ;
     527    }
     528    if (poly->COOL_4D_nT > maxChebyPoly) {
     529        maxChebyPoly = poly->COOL_4D_nT;
    529530    }
    530531    chebPolys = createChebyshevPolys(maxChebyPoly);
    531532
    532     for (loop_x = 0; loop_x < poly->nX; loop_x++) {
    533         for (loop_y = 0; loop_y < poly->nY; loop_y++) {
    534             for (loop_z = 0; loop_z < poly->nZ; loop_z++) {
    535                 for (loop_t = 0; loop_t < poly->nT; loop_t++) {
     533    for (loop_x = 0; loop_x < (1 + poly->COOL_4D_nX); loop_x++) {
     534        for (loop_y = 0; loop_y < (1 + poly->COOL_4D_nY); loop_y++) {
     535            for (loop_z = 0; loop_z < (1 + poly->COOL_4D_nZ); loop_z++) {
     536                for (loop_t = 0; loop_t < (1 + poly->COOL_4D_nZ); loop_t++) {
    536537                    if (poly->mask[loop_x][loop_y][loop_z][loop_t] == 0) {
    537538                        polySum += poly->coeff[loop_x][loop_y][loop_z][loop_t] *
    … …  
    610611/*****************************************************************************
    611612    This routine must allocate memory for the polynomial structures.
     613    XXX: replaces nterms variables with nOrder.  Lots of potential for bugs
     614         here.  Probably should create separately named private variables
     615         in these functions.
    612616 *****************************************************************************/
    613617psPolynomial1D* psPolynomial1DAlloc(unsigned int n,
    … …  
    623627
    624628    newPoly->type = type;
    625     newPoly->n = n;
     629    newPoly->COOL_1D_n = n;
    626630    newPoly->coeff = psAlloc(n * sizeof(psF64));
    627631    newPoly->coeffErr = psAlloc(n * sizeof(psF64));
    628     newPoly->mask = (psMaskType *)psAlloc(n * sizeof(psMaskType));
    629     for (i = 0; i < n; i++) {
     632    newPoly->mask = (psMaskType *)psAlloc((n + 1) * sizeof(psMaskType));
     633    for (i = 0; i < (n + 1); i++) {
    630634        newPoly->coeff[i] = 0.0;
    631635        newPoly->coeffErr[i] = 0.0;
    … …  
    651655
    652656    newPoly->type = type;
    653     newPoly->nX = nX;
    654     newPoly->nY = nY;
    655 
    656     newPoly->coeff = psAlloc(nX * sizeof(psF64 *));
    657     newPoly->coeffErr = psAlloc(nX * sizeof(psF64 *));
    658     newPoly->mask = (psMaskType **)psAlloc(nX * sizeof(psMaskType *));
    659     for (x = 0; x < nX; x++) {
    660         newPoly->coeff[x] = psAlloc(nY * sizeof(psF64));
    661         newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64));
    662         newPoly->mask[x] = (psMaskType *)psAlloc(nY * sizeof(psMaskType));
    663     }
    664     for (x = 0; x < nX; x++) {
    665         for (y = 0; y < nY; y++) {
     657    newPoly->COOL_2D_nX = nX;
     658    newPoly->COOL_2D_nY = nY;
     659
     660    newPoly->coeff = psAlloc((1 + nX) * sizeof(psF64 *));
     661    newPoly->coeffErr = psAlloc((1 + nX) * sizeof(psF64 *));
     662    newPoly->mask = (psMaskType **)psAlloc((1 + nX) * sizeof(psMaskType *));
     663    for (x = 0; x < (1 + nX); x++) {
     664        newPoly->coeff[x] = psAlloc((1 + nY) * sizeof(psF64));
     665        newPoly->coeffErr[x] = psAlloc((1 + nY) * sizeof(psF64));
     666        newPoly->mask[x] = (psMaskType *)psAlloc((1 + nY) * sizeof(psMaskType));
     667    }
     668    for (x = 0; x < (1 + nX); x++) {
     669        for (y = 0; y < (1 + nY); y++) {
    666670            newPoly->coeff[x][y] = 0.0;
    667671            newPoly->coeffErr[x][y] = 0.0;
    … …  
    691695
    692696    newPoly->type = type;
    693     newPoly->nX = nX;
    694     newPoly->nY = nY;
    695     newPoly->nZ = nZ;
    696 
    697     newPoly->coeff = psAlloc(nX * sizeof(psF64 **));
    698     newPoly->coeffErr = psAlloc(nX * sizeof(psF64 **));
    699     newPoly->mask = (psMaskType ***)psAlloc(nX * sizeof(psMaskType **));
    700     for (x = 0; x < nX; x++) {
    701         newPoly->coeff[x] = psAlloc(nY * sizeof(psF64 *));
    702         newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64 *));
    703         newPoly->mask[x] = (psMaskType **)psAlloc(nY * sizeof(psMaskType *));
    704         for (y = 0; y < nY; y++) {
    705             newPoly->coeff[x][y] = psAlloc(nZ * sizeof(psF64));
    706             newPoly->coeffErr[x][y] = psAlloc(nZ * sizeof(psF64));
    707             newPoly->mask[x][y] = (psMaskType *)psAlloc(nZ * sizeof(psMaskType));
    708         }
    709     }
    710     for (x = 0; x < nX; x++) {
    711         for (y = 0; y < nY; y++) {
    712             for (z = 0; z < nZ; z++) {
     697    newPoly->COOL_3D_nX = nX;
     698    newPoly->COOL_3D_nY = nY;
     699    newPoly->COOL_3D_nZ = nZ;
     700
     701    newPoly->coeff = psAlloc((nX + 1) * sizeof(psF64 **));
     702    newPoly->coeffErr = psAlloc((nX + 1) * sizeof(psF64 **));
     703    newPoly->mask = (psMaskType ***)psAlloc((nX + 1) * sizeof(psMaskType **));
     704    for (x = 0; x < (1 + nX); x++) {
     705        newPoly->coeff[x] = psAlloc((nY + 1) * sizeof(psF64 *));
     706        newPoly->coeffErr[x] = psAlloc((nY + 1) * sizeof(psF64 *));
     707        newPoly->mask[x] = (psMaskType **)psAlloc((nY + 1) * sizeof(psMaskType *));
     708        for (y = 0; y < (nY + 1); y++) {
     709            newPoly->coeff[x][y] = psAlloc((nZ + 1) * sizeof(psF64));
     710            newPoly->coeffErr[x][y] = psAlloc((nZ + 1) * sizeof(psF64));
     711            newPoly->mask[x][y] = (psMaskType *)psAlloc((nZ + 1) * sizeof(psMaskType));
     712        }
     713    }
     714    for (x = 0; x < (nX + 1); x++) {
     715        for (y = 0; y < (nY + 1); y++) {
     716            for (z = 0; z < (nZ + 1); z++) {
    713717                newPoly->coeff[x][y][z] = 0.0;
    714718                newPoly->coeffErr[x][y][z] = 0.0;
    … …  
    742746
    743747    newPoly->type = type;
    744     newPoly->nX = nX;
    745     newPoly->nY = nY;
    746     newPoly->nZ = nZ;
    747     newPoly->nT = nT;
    748 
    749     newPoly->coeff = psAlloc(nX * sizeof(psF64 ***));
    750     newPoly->coeffErr = psAlloc(nX * sizeof(psF64 ***));
    751     newPoly->mask = (psMaskType ****)psAlloc(nX * sizeof(psMaskType ***));
    752     for (x = 0; x < nX; x++) {
    753         newPoly->coeff[x] = psAlloc(nY * sizeof(psF64 **));
    754         newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64 **));
    755         newPoly->mask[x] = (psMaskType ***)psAlloc(nY * sizeof(psMaskType **));
    756         for (y = 0; y < nY; y++) {
    757             newPoly->coeff[x][y] = psAlloc(nZ * sizeof(psF64 *));
    758             newPoly->coeffErr[x][y] = psAlloc(nZ * sizeof(psF64 *));
    759             newPoly->mask[x][y] = (psMaskType **)psAlloc(nZ * sizeof(psMaskType *));
    760             for (z = 0; z < nZ; z++) {
    761                 newPoly->coeff[x][y][z] = psAlloc(nT * sizeof(psF64));
    762                 newPoly->coeffErr[x][y][z] = psAlloc(nT * sizeof(psF64));
    763                 newPoly->mask[x][y][z] = (psMaskType *)psAlloc(nT * sizeof(psMaskType));
    764             }
    765         }
    766     }
    767     for (x = 0; x < nX; x++) {
    768         for (y = 0; y < nY; y++) {
    769             for (z = 0; z < nZ; z++) {
    770                 for (t = 0; t < nT; t++) {
     748    newPoly->COOL_4D_nX = nX;
     749    newPoly->COOL_4D_nY = nY;
     750    newPoly->COOL_4D_nZ = nZ;
     751    newPoly->COOL_4D_nT = nT;
     752
     753    newPoly->coeff = psAlloc((nX + 1) * sizeof(psF64 ***));
     754    newPoly->coeffErr = psAlloc((nX + 1) * sizeof(psF64 ***));
     755    newPoly->mask = (psMaskType ****)psAlloc((nX + 1) * sizeof(psMaskType ***));
     756    for (x = 0; x < (nX + 1); x++) {
     757        newPoly->coeff[x] = psAlloc((nY + 1) * sizeof(psF64 **));
     758        newPoly->coeffErr[x] = psAlloc((nY + 1) * sizeof(psF64 **));
     759        newPoly->mask[x] = (psMaskType ***)psAlloc((nY + 1) * sizeof(psMaskType **));
     760        for (y = 0; y < (nY + 1); y++) {
     761            newPoly->coeff[x][y] = psAlloc((nZ + 1) * sizeof(psF64 *));
     762            newPoly->coeffErr[x][y] = psAlloc((nZ + 1) * sizeof(psF64 *));
     763            newPoly->mask[x][y] = (psMaskType **)psAlloc((nZ + 1) * sizeof(psMaskType *));
     764            for (z = 0; z < (nZ + 1); z++) {
     765                newPoly->coeff[x][y][z] = psAlloc((nT + 1) * sizeof(psF64));
     766                newPoly->coeffErr[x][y][z] = psAlloc((nT + 1) * sizeof(psF64));
     767                newPoly->mask[x][y][z] = (psMaskType *)psAlloc((nT + 1) * sizeof(psMaskType));
     768            }
     769        }
     770    }
     771    for (x = 0; x < (nX + 1); x++) {
     772        for (y = 0; y < (nY + 1); y++) {
     773            for (z = 0; z < (nZ + 1); z++) {
     774                for (t = 0; t < (nT + 1); t++) {
    771775                    newPoly->coeff[x][y][z][t] = 0.0;
    772776                    newPoly->coeffErr[x][y][z][t] = 0.0;
    … …  
    972976    }
    973977
    974     // Allocate output vector
     978    // Allocoutput vector
    975979    tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
    976980

  • trunk/psLib/src/math/psPolynomial.h

    r5066 r5090  
    1111 *  @author GLG, MHPCC
    1212 *
    13  *  @version $Revision: 1.56 $ $Name: not supported by cvs2svn $
    14  *  @date $Date: 2005-09-19 19:53:13 $
     13 *  @version $Revision: 1.57 $ $Name: not supported by cvs2svn $
     14 *  @date $Date: 2005-09-22 02:47:16 $
    1515 *
    1616 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    6464 */
    6565typedef enum {
    66     PS_POLYNOMIAL_ORD,                 ///< Ordinary Polynomial
    67     PS_POLYNOMIAL_CHEB                 ///< Chebyshev Polynomial
     66    PS_POLYNOMIAL_ORD,                  ///< Ordinary Polynomial
     67    PS_POLYNOMIAL_CHEB                  ///< Chebyshev Polynomial
    6868}
    6969psPolynomialType;
    7070
     71// XXX: These are incorrect names for the order of the polynomial.  We
     72// keep them here temporarily so we can later sed replace them with the
     73// correct names.
    7174/** One-dimensional polynomial */
    7275typedef struct
    7376{
    74     psPolynomialType type;             ///< Polynomial type
    75     unsigned int n;                    ///< Number of terms
    76     psF64 *coeff;                      ///< Coefficients
    77     psF64 *coeffErr;                   ///< Error in coefficients
    78     psMaskType *mask;                  ///< Coefficient mask
     77    psPolynomialType type;              ///< Polynomial type
     78    unsigned int COOL_1D_n;             ///< Polynomial order
     79    psF64 *coeff;                       ///< Coefficients
     80    psF64 *coeffErr;                    ///< Error in coefficients
     81    psMaskType *mask;                   ///< Coefficient mask
    7982}
    8083psPolynomial1D;
    8184
     85// XXX: These are incorrect names for the order of the polynomial.  We
     86// keep them here temporarily so we can later sed replace them with the
     87// correct names.
    8288/** Two-dimensional polynomial */
    8389typedef struct
    8490{
    85     psPolynomialType type;             ///< Polynomial type
    86     unsigned int nX;                   ///< Number of terms in x
    87     unsigned int nY;                   ///< Number of terms in y
    88     psF64 **coeff;                     ///< Coefficients
    89     psF64 **coeffErr;                  ///< Error in coefficients
    90     psMaskType **mask;                 ///< Coefficients mask
     91    psPolynomialType type;              ///< Polynomial type
     92    unsigned int COOL_2D_nX;            ///< Polynomial order in x
     93    unsigned int COOL_2D_nY;            ///< Polynomial order in y
     94    psF64 **coeff;                      ///< Coefficients
     95    psF64 **coeffErr;                   ///< Error in coefficients
     96    psMaskType **mask;                  ///< Coefficients mask
    9197}
    9298psPolynomial2D;
    9399
     100// XXX: These are incorrect names for the order of the polynomial.  We
     101// keep them here temporarily so we can later sed replace them with the
     102// correct names.
    94103/** Three-dimensional polynomial */
    95104typedef struct
    96105{
    97     psPolynomialType type;             ///< Polynomial type
    98     unsigned int nX;                   ///< Number of terms in x
    99     unsigned int nY;                   ///< Number of terms in y
    100     unsigned int nZ;                   ///< Number of terms in z
    101     psF64 ***coeff;                    ///< Coefficients
    102     psF64 ***coeffErr;                 ///< Error in coefficients
    103     psMaskType ***mask;                ///< Coefficients mask
     106    psPolynomialType type;              ///< Polynomial type
     107    unsigned int COOL_3D_nX;            ///< Polynomial order in x
     108    unsigned int COOL_3D_nY;            ///< Polynomial order in y
     109    unsigned int COOL_3D_nZ;            ///< Polynomial order in z
     110    psF64 ***coeff;                     ///< Coefficients
     111    psF64 ***coeffErr;                  ///< Error in coefficients
     112    psMaskType ***mask;                 ///< Coefficients mask
    104113}
    105114psPolynomial3D;
    106115
     116// XXX: These are incorrect names for the order of the polynomial.  We
     117// keep them here temporarily so we can later sed replace them with the
     118// correct names.
    107119/** Four-dimensional polynomial */
    108120typedef struct
    109121{
    110     psPolynomialType type;             ///< Polynomial type
    111     unsigned int nX;                   ///< Number of terms in x
    112     unsigned int nY;                   ///< Number of terms in y
    113     unsigned int nZ;                   ///< Number of terms in z
    114     unsigned int nT;                   ///< Number of terms in t
    115     psF64 ****coeff;                   ///< Coefficients
    116     psF64 ****coeffErr;                ///< Error in coefficients
    117     psMaskType ****mask;               ///< Coefficients mask
     122    psPolynomialType type;              ///< Polynomial type
     123    unsigned int COOL_4D_nX;            ///< Polynomial order in x
     124    unsigned int COOL_4D_nY;            ///< Polynomial order in y
     125    unsigned int COOL_4D_nZ;            ///< Polynomial order in z
     126    unsigned int COOL_4D_nT;            ///< Polynomial order in t
     127    psF64 ****coeff;                    ///< Coefficients
     128    psF64 ****coeffErr;                 ///< Error in coefficients
     129    psMaskType ****mask;                ///< Coefficients mask
    118130}
    119131psPolynomial4D;

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

    r5072 r5090  
    77*  splines.
    88*
    9 *  @version $Revision: 1.126 $ $Name: not supported by cvs2svn $
    10 *  @date $Date: 2005-09-20 02:43:53 $
     9*  @version $Revision: 1.127 $ $Name: not supported by cvs2svn $
     10*  @date $Date: 2005-09-22 02:47:16 $
    1111*
    1212*
    … …  
    383383    // If these are linear splines, which means their polynomials will have
    384384    // two coefficients, then we do the simple calculation.
    385     if (2 == (spline->spline[0])->n) {
     385    if (1 == (spline->spline[0])->COOL_1D_n) {
    386386        for (i=0;i<spline->n;i++) {
    387387            slope = (y32->data.F32[i+1] - y32->data.F32[i]) /
    … …  
    402402
    403403    // Check if these are cubic splines (n==4).  If not, psError.
    404     if (4 != (spline->spline[0])->n) {
     404    if (3 != (spline->spline[0])->COOL_1D_n) {
    405405        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
    406406                "Don't know how to generate %u-order splines.",
    407                 (spline->spline[0])->n-1);
     407                (spline->spline[0])->COOL_1D_n);
    408408        return(NULL);
    409409    }
    … …  
    745745    tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
    746746    for (unsigned int i=0;i<numSplines;i++) {
    747         (tmpSpline->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
     747        (tmpSpline->spline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD);
    748748    }
    749749
    … …  
    809809    tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *));
    810810    for (unsigned int i=0;i<numSplines;i++) {
    811         (tmpSpline->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD);
     811        (tmpSpline->spline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD);
    812812    }
    813813

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

    r4991 r5090  
    1717 *
    1818 *
    19  *  @version $Revision: 1.145 $ $Name: not supported by cvs2svn $
    20  *  @date $Date: 2005-09-11 22:18:40 $
     19 *  @version $Revision: 1.146 $ $Name: not supported by cvs2svn $
     20 *  @date $Date: 2005-09-22 02:47:16 $
    2121 *
    2222 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    15021502    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
    15031503    psVector *yErr = psVectorAlloc(3, PS_TYPE_F64);
     1504    // XXX: Why was this 2 when the alloc function specified number of terms?  Note: it's correct now.
    15041505    psPolynomial1D *myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);
    15051506
    … …  
    15521553
    15531554        // Determine the coefficients of the polynomial.
    1554         //        myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
    15551555        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, yErr, x);
    15561556        if (myPoly == NULL) {
    … …  
    18271827        }
    18281828
    1829         psPolynomial1D *tmpPoly = psPolynomial1DAlloc(3, PS_POLYNOMIAL_ORD);
     1829        psPolynomial1D *tmpPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);
    18301830        // XXX: What about the NULL x argument?
    18311831        tmpPoly = psVectorFitPolynomial1D(tmpPoly, NULL, 0, y, NULL, NULL);
Note: See TracChangeset for help on using the changeset viewer.