IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.