Index: trunk/psLib/src/math/psMinimize.c =================================================================== --- trunk/psLib/src/math/psMinimize.c (revision 5074) +++ trunk/psLib/src/math/psMinimize.c (revision 5090) @@ -10,6 +10,6 @@ * @author EAM, IfA * - * @version $Revision: 1.138 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-13 21:43:47 $ + * @version $Revision: 1.139 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-22 02:47:16 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -1459,5 +1459,5 @@ // XXX: these ASSERTS are redundant. PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->n, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_1D_n, NULL); PS_ASSERT_VECTOR_NON_NULL(y, NULL); PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); @@ -1523,5 +1523,5 @@ fac = 2.0/((psF32) n); // XXX: is this loop bound correct? - for (j=0;jn;j++) { + for (j=0;jCOOL_1D_n+1;j++) { sum = 0.0; for (k=0;kn, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_1D_n, NULL); PS_ASSERT_VECTOR_NON_NULL(f, NULL); PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL); @@ -1597,5 +1597,5 @@ } - nTerm = myPoly->n; + nTerm = 1 + myPoly->COOL_1D_n; nOrder = nTerm - 1; @@ -1697,5 +1697,5 @@ PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_INT_NONNEGATIVE(poly->n, NULL); + PS_ASSERT_INT_NONNEGATIVE(poly->COOL_1D_n, NULL); PS_ASSERT_VECTOR_NON_NULL(f, NULL); PS_ASSERT_VECTOR_NON_EMPTY(f, NULL); @@ -1837,6 +1837,6 @@ // These ASSERTS are redundant. PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_2D_nX, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_2D_nY, NULL); PS_ASSERT_VECTOR_NON_NULL(f, NULL); @@ -1865,6 +1865,6 @@ // XXX:Watch for changes to the psPolys: nTerm != nOrder. - psS32 nXterm = myPoly->nX; - psS32 nYterm = myPoly->nY; + psS32 nXterm = 1 + myPoly->COOL_2D_nX; + psS32 nYterm = 1 + myPoly->COOL_2D_nY; nTerm = nXterm * nYterm; @@ -2264,7 +2264,7 @@ // These ASSERTS are redundant. PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_3D_nX, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_3D_nY, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_3D_nZ, NULL); PS_ASSERT_VECTOR_NON_NULL(f, NULL); @@ -2527,8 +2527,8 @@ // These ASSERTS are redundant. PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nT, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_4D_nX, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_4D_nY, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_4D_nZ, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->COOL_4D_nT, NULL); PS_ASSERT_VECTOR_NON_NULL(f, NULL); PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL); Index: trunk/psLib/src/math/psPolynomial.c =================================================================== --- trunk/psLib/src/math/psPolynomial.c (revision 5074) +++ trunk/psLib/src/math/psPolynomial.c (revision 5090) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * -* @version $Revision: 1.122 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-09-20 20:48:35 $ +* @version $Revision: 1.123 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-09-22 02:47:16 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -95,5 +95,5 @@ unsigned int x = 0; - for (x = 0; x < poly->nX; x++) { + for (x = 0; x < poly->COOL_2D_nX; x++) { psFree(poly->coeff[x]); psFree(poly->coeffErr[x]); @@ -110,6 +110,6 @@ unsigned int y = 0; - for (x = 0; x < poly->nX; x++) { - for (y = 0; y < poly->nY; y++) { + for (x = 0; x < poly->COOL_3D_nX; x++) { + for (y = 0; y < poly->COOL_3D_nY; y++) { psFree(poly->coeff[x][y]); psFree(poly->coeffErr[x][y]); @@ -132,7 +132,7 @@ unsigned int z = 0; - for (x = 0; x < poly->nX; x++) { - for (y = 0; y < poly->nY; y++) { - for (z = 0; z < poly->nZ; z++) { + for (x = 0; x < poly->COOL_4D_nX; x++) { + for (y = 0; y < poly->COOL_4D_nY; y++) { + for (z = 0; z < poly->COOL_4D_nZ; z++) { psFree(poly->coeff[x][y][z]); psFree(poly->coeffErr[x][y][z]); @@ -170,5 +170,5 @@ chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *)); for (unsigned int i = 0; i < maxChebyPoly; i++) { - chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD); + chebPolys[i] = psPolynomial1DAlloc(i, PS_POLYNOMIAL_ORD); } @@ -182,8 +182,8 @@ for (psS32 i = 2; i < maxChebyPoly; i++) { - for (psS32 j = 0; j < chebPolys[i - 1]->n; j++) { + for (psS32 j = 0; j < chebPolys[i - 1]->COOL_1D_n; j++) { chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j]; } - for (psS32 j = 0; j < chebPolys[i - 2]->n; j++) { + for (psS32 j = 0; j < chebPolys[i - 2]->COOL_1D_n; j++) { chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j]; } @@ -210,11 +210,11 @@ "---- Calling ordPolynomial1DEval(%lf)\n", x); psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 4, - "Polynomial order is %u\n", poly->n); - for (loop_x = 0; loop_x < poly->n; loop_x++) { + "Polynomial order is %u\n", poly->COOL_1D_n); + for (loop_x = 0; loop_x < poly->COOL_1D_n+1; loop_x++) { psTrace(".psLib.dataManip.psPolynomial.ordPolynomial1DEval", 4, "Polynomial coeff[%u] is %lf\n", loop_x, poly->coeff[loop_x]); } - for (loop_x = 0; loop_x < poly->n; loop_x++) { + for (loop_x = 0; loop_x < poly->COOL_1D_n+1; loop_x++) { if (poly->mask[loop_x] == 0) { // XXX: If you set the tracelevel to 10 here, and later set the tracelevel to @@ -238,10 +238,11 @@ PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0); // XXX: Create a macro for this in psConstants.h - if (poly->n < 1) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error: Chebyshev polynomial as order %u.", poly->n); + if (poly->COOL_1D_n < 1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error: Chebyshev polynomial is order %u.", poly->COOL_1D_n); return(NAN); } psVector *d; - unsigned int n = poly->n; + // XXX: n should be nTerms here (for clarity). Or get rid of the variable. + unsigned int n = 1 + poly->COOL_1D_n; unsigned int i; psF64 tmp = 0.0; @@ -302,9 +303,9 @@ psPolynomial1D **chebPolys = NULL; - n = poly->n; + n = 1 + poly->COOL_1D_n; chebPolys = createChebyshevPolys(n); tmp = 0.0; - for (i=0;in;i++) { + for (i=0;i<(1 + poly->COOL_1D_n);i++) { tmp+= (poly->coeff[i] * psPolynomial1DEval(x, chebPolys[i])); } @@ -328,7 +329,7 @@ psF64 ySum = 1.0; - for (loop_x = 0; loop_x < poly->nX; loop_x++) { + for (loop_x = 0; loop_x < (1 + poly->COOL_2D_nX); loop_x++) { ySum = xSum; - for (loop_y = 0; loop_y < poly->nY; loop_y++) { + for (loop_y = 0; loop_y < (1 + poly->COOL_2D_nY); loop_y++) { if (poly->mask[loop_x][loop_y] == 0) { polySum += ySum * poly->coeff[loop_x][loop_y]; @@ -359,12 +360,12 @@ // Determine how many Chebyshev polynomials // are needed, then create them. - maxChebyPoly = poly->nX; - if (poly->nY > maxChebyPoly) { - maxChebyPoly = poly->nY; + maxChebyPoly = 1 + poly->COOL_2D_nX; + if ((1 + poly->COOL_2D_nY) > maxChebyPoly) { + maxChebyPoly = (1 + poly->COOL_2D_nY); } chebPolys = createChebyshevPolys(maxChebyPoly); - for (loop_x = 0; loop_x < poly->nX; loop_x++) { - for (loop_y = 0; loop_y < poly->nY; loop_y++) { + for (loop_x = 0; loop_x < (1 + poly->COOL_2D_nX); loop_x++) { + for (loop_y = 0; loop_y < (1 + poly->COOL_2D_nY); loop_y++) { if (poly->mask[loop_x][loop_y] == 0) { polySum += poly->coeff[loop_x][loop_y] * @@ -394,9 +395,9 @@ psF64 zSum = 1.0; - for (loop_x = 0; loop_x < poly->nX; loop_x++) { + for (loop_x = 0; loop_x < (1 + poly->COOL_3D_nX); loop_x++) { ySum = xSum; - for (loop_y = 0; loop_y < poly->nY; loop_y++) { + for (loop_y = 0; loop_y < (1 + poly->COOL_3D_nY); loop_y++) { zSum = ySum; - for (loop_z = 0; loop_z < poly->nZ; loop_z++) { + for (loop_z = 0; loop_z < (1 + poly->COOL_3D_nZ); loop_z++) { if (poly->mask[loop_x][loop_y][loop_z] == 0) { polySum += zSum * poly->coeff[loop_x][loop_y][loop_z]; @@ -430,16 +431,16 @@ // Determine how many Chebyshev polynomials // are needed, then create them. - maxChebyPoly = poly->nX; - if (poly->nY > maxChebyPoly) { - maxChebyPoly = poly->nY; - } - if (poly->nZ > maxChebyPoly) { - maxChebyPoly = poly->nZ; - } - chebPolys = createChebyshevPolys(maxChebyPoly); - - for (loop_x = 0; loop_x < poly->nX; loop_x++) { - for (loop_y = 0; loop_y < poly->nY; loop_y++) { - for (loop_z = 0; loop_z < poly->nZ; loop_z++) { + maxChebyPoly = poly->COOL_3D_nX; + if (poly->COOL_3D_nY > maxChebyPoly) { + maxChebyPoly = poly->COOL_3D_nY; + } + if (poly->COOL_3D_nZ > maxChebyPoly) { + maxChebyPoly = poly->COOL_3D_nZ; + } + chebPolys = createChebyshevPolys(maxChebyPoly + 1); + + for (loop_x = 0; loop_x < (1 + poly->COOL_3D_nX); loop_x++) { + for (loop_y = 0; loop_y < (1 + poly->COOL_3D_nY); loop_y++) { + for (loop_z = 0; loop_z < (1 + poly->COOL_3D_nZ); loop_z++) { if (poly->mask[loop_x][loop_y][loop_z] == 0) { polySum += poly->coeff[loop_x][loop_y][loop_z] * @@ -475,11 +476,11 @@ psF64 tSum = 1.0; - for (loop_x = 0; loop_x < poly->nX; loop_x++) { + for (loop_x = 0; loop_x < (1 + poly->COOL_4D_nX); loop_x++) { ySum = xSum; - for (loop_y = 0; loop_y < poly->nY; loop_y++) { + for (loop_y = 0; loop_y < (1 + poly->COOL_4D_nY); loop_y++) { zSum = ySum; - for (loop_z = 0; loop_z < poly->nZ; loop_z++) { + for (loop_z = 0; loop_z < (1 + poly->COOL_4D_nZ); loop_z++) { tSum = zSum; - for (loop_t = 0; loop_t < poly->nT; loop_t++) { + for (loop_t = 0; loop_t < (1 + poly->COOL_4D_nT); loop_t++) { if (poly->mask[loop_x][loop_y][loop_z][loop_t] == 0) { polySum += tSum * poly->coeff[loop_x][loop_y][loop_z][loop_t]; @@ -518,20 +519,20 @@ // Determine how many Chebyshev polynomials // are needed, then create them. - maxChebyPoly = poly->nX; - if (poly->nY > maxChebyPoly) { - maxChebyPoly = poly->nY; - } - if (poly->nZ > maxChebyPoly) { - maxChebyPoly = poly->nZ; - } - if (poly->nT > maxChebyPoly) { - maxChebyPoly = poly->nT; + maxChebyPoly = poly->COOL_4D_nX; + if (poly->COOL_4D_nY > maxChebyPoly) { + maxChebyPoly = poly->COOL_4D_nY; + } + if (poly->COOL_4D_nZ > maxChebyPoly) { + maxChebyPoly = poly->COOL_4D_nZ; + } + if (poly->COOL_4D_nT > maxChebyPoly) { + maxChebyPoly = poly->COOL_4D_nT; } chebPolys = createChebyshevPolys(maxChebyPoly); - for (loop_x = 0; loop_x < poly->nX; loop_x++) { - for (loop_y = 0; loop_y < poly->nY; loop_y++) { - for (loop_z = 0; loop_z < poly->nZ; loop_z++) { - for (loop_t = 0; loop_t < poly->nT; loop_t++) { + for (loop_x = 0; loop_x < (1 + poly->COOL_4D_nX); loop_x++) { + for (loop_y = 0; loop_y < (1 + poly->COOL_4D_nY); loop_y++) { + for (loop_z = 0; loop_z < (1 + poly->COOL_4D_nZ); loop_z++) { + for (loop_t = 0; loop_t < (1 + poly->COOL_4D_nZ); loop_t++) { if (poly->mask[loop_x][loop_y][loop_z][loop_t] == 0) { polySum += poly->coeff[loop_x][loop_y][loop_z][loop_t] * @@ -610,4 +611,7 @@ /***************************************************************************** This routine must allocate memory for the polynomial structures. + XXX: replaces nterms variables with nOrder. Lots of potential for bugs + here. Probably should create separately named private variables + in these functions. *****************************************************************************/ psPolynomial1D* psPolynomial1DAlloc(unsigned int n, @@ -623,9 +627,9 @@ newPoly->type = type; - newPoly->n = n; + newPoly->COOL_1D_n = n; newPoly->coeff = psAlloc(n * sizeof(psF64)); newPoly->coeffErr = psAlloc(n * sizeof(psF64)); - newPoly->mask = (psMaskType *)psAlloc(n * sizeof(psMaskType)); - for (i = 0; i < n; i++) { + newPoly->mask = (psMaskType *)psAlloc((n + 1) * sizeof(psMaskType)); + for (i = 0; i < (n + 1); i++) { newPoly->coeff[i] = 0.0; newPoly->coeffErr[i] = 0.0; @@ -651,17 +655,17 @@ newPoly->type = type; - newPoly->nX = nX; - newPoly->nY = nY; - - newPoly->coeff = psAlloc(nX * sizeof(psF64 *)); - newPoly->coeffErr = psAlloc(nX * sizeof(psF64 *)); - newPoly->mask = (psMaskType **)psAlloc(nX * sizeof(psMaskType *)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = psAlloc(nY * sizeof(psF64)); - newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64)); - newPoly->mask[x] = (psMaskType *)psAlloc(nY * sizeof(psMaskType)); - } - for (x = 0; x < nX; x++) { - for (y = 0; y < nY; y++) { + newPoly->COOL_2D_nX = nX; + newPoly->COOL_2D_nY = nY; + + newPoly->coeff = psAlloc((1 + nX) * sizeof(psF64 *)); + newPoly->coeffErr = psAlloc((1 + nX) * sizeof(psF64 *)); + newPoly->mask = (psMaskType **)psAlloc((1 + nX) * sizeof(psMaskType *)); + for (x = 0; x < (1 + nX); x++) { + newPoly->coeff[x] = psAlloc((1 + nY) * sizeof(psF64)); + newPoly->coeffErr[x] = psAlloc((1 + nY) * sizeof(psF64)); + newPoly->mask[x] = (psMaskType *)psAlloc((1 + nY) * sizeof(psMaskType)); + } + for (x = 0; x < (1 + nX); x++) { + for (y = 0; y < (1 + nY); y++) { newPoly->coeff[x][y] = 0.0; newPoly->coeffErr[x][y] = 0.0; @@ -691,24 +695,24 @@ newPoly->type = type; - newPoly->nX = nX; - newPoly->nY = nY; - newPoly->nZ = nZ; - - newPoly->coeff = psAlloc(nX * sizeof(psF64 **)); - newPoly->coeffErr = psAlloc(nX * sizeof(psF64 **)); - newPoly->mask = (psMaskType ***)psAlloc(nX * sizeof(psMaskType **)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = psAlloc(nY * sizeof(psF64 *)); - newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64 *)); - newPoly->mask[x] = (psMaskType **)psAlloc(nY * sizeof(psMaskType *)); - for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = psAlloc(nZ * sizeof(psF64)); - newPoly->coeffErr[x][y] = psAlloc(nZ * sizeof(psF64)); - newPoly->mask[x][y] = (psMaskType *)psAlloc(nZ * sizeof(psMaskType)); - } - } - for (x = 0; x < nX; x++) { - for (y = 0; y < nY; y++) { - for (z = 0; z < nZ; z++) { + newPoly->COOL_3D_nX = nX; + newPoly->COOL_3D_nY = nY; + newPoly->COOL_3D_nZ = nZ; + + newPoly->coeff = psAlloc((nX + 1) * sizeof(psF64 **)); + newPoly->coeffErr = psAlloc((nX + 1) * sizeof(psF64 **)); + newPoly->mask = (psMaskType ***)psAlloc((nX + 1) * sizeof(psMaskType **)); + for (x = 0; x < (1 + nX); x++) { + newPoly->coeff[x] = psAlloc((nY + 1) * sizeof(psF64 *)); + newPoly->coeffErr[x] = psAlloc((nY + 1) * sizeof(psF64 *)); + newPoly->mask[x] = (psMaskType **)psAlloc((nY + 1) * sizeof(psMaskType *)); + for (y = 0; y < (nY + 1); y++) { + newPoly->coeff[x][y] = psAlloc((nZ + 1) * sizeof(psF64)); + newPoly->coeffErr[x][y] = psAlloc((nZ + 1) * sizeof(psF64)); + newPoly->mask[x][y] = (psMaskType *)psAlloc((nZ + 1) * sizeof(psMaskType)); + } + } + for (x = 0; x < (nX + 1); x++) { + for (y = 0; y < (nY + 1); y++) { + for (z = 0; z < (nZ + 1); z++) { newPoly->coeff[x][y][z] = 0.0; newPoly->coeffErr[x][y][z] = 0.0; @@ -742,31 +746,31 @@ newPoly->type = type; - newPoly->nX = nX; - newPoly->nY = nY; - newPoly->nZ = nZ; - newPoly->nT = nT; - - newPoly->coeff = psAlloc(nX * sizeof(psF64 ***)); - newPoly->coeffErr = psAlloc(nX * sizeof(psF64 ***)); - newPoly->mask = (psMaskType ****)psAlloc(nX * sizeof(psMaskType ***)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = psAlloc(nY * sizeof(psF64 **)); - newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64 **)); - newPoly->mask[x] = (psMaskType ***)psAlloc(nY * sizeof(psMaskType **)); - for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = psAlloc(nZ * sizeof(psF64 *)); - newPoly->coeffErr[x][y] = psAlloc(nZ * sizeof(psF64 *)); - newPoly->mask[x][y] = (psMaskType **)psAlloc(nZ * sizeof(psMaskType *)); - for (z = 0; z < nZ; z++) { - newPoly->coeff[x][y][z] = psAlloc(nT * sizeof(psF64)); - newPoly->coeffErr[x][y][z] = psAlloc(nT * sizeof(psF64)); - newPoly->mask[x][y][z] = (psMaskType *)psAlloc(nT * sizeof(psMaskType)); - } - } - } - for (x = 0; x < nX; x++) { - for (y = 0; y < nY; y++) { - for (z = 0; z < nZ; z++) { - for (t = 0; t < nT; t++) { + newPoly->COOL_4D_nX = nX; + newPoly->COOL_4D_nY = nY; + newPoly->COOL_4D_nZ = nZ; + newPoly->COOL_4D_nT = nT; + + newPoly->coeff = psAlloc((nX + 1) * sizeof(psF64 ***)); + newPoly->coeffErr = psAlloc((nX + 1) * sizeof(psF64 ***)); + newPoly->mask = (psMaskType ****)psAlloc((nX + 1) * sizeof(psMaskType ***)); + for (x = 0; x < (nX + 1); x++) { + newPoly->coeff[x] = psAlloc((nY + 1) * sizeof(psF64 **)); + newPoly->coeffErr[x] = psAlloc((nY + 1) * sizeof(psF64 **)); + newPoly->mask[x] = (psMaskType ***)psAlloc((nY + 1) * sizeof(psMaskType **)); + for (y = 0; y < (nY + 1); y++) { + newPoly->coeff[x][y] = psAlloc((nZ + 1) * sizeof(psF64 *)); + newPoly->coeffErr[x][y] = psAlloc((nZ + 1) * sizeof(psF64 *)); + newPoly->mask[x][y] = (psMaskType **)psAlloc((nZ + 1) * sizeof(psMaskType *)); + for (z = 0; z < (nZ + 1); z++) { + newPoly->coeff[x][y][z] = psAlloc((nT + 1) * sizeof(psF64)); + newPoly->coeffErr[x][y][z] = psAlloc((nT + 1) * sizeof(psF64)); + newPoly->mask[x][y][z] = (psMaskType *)psAlloc((nT + 1) * sizeof(psMaskType)); + } + } + } + for (x = 0; x < (nX + 1); x++) { + for (y = 0; y < (nY + 1); y++) { + for (z = 0; z < (nZ + 1); z++) { + for (t = 0; t < (nT + 1); t++) { newPoly->coeff[x][y][z][t] = 0.0; newPoly->coeffErr[x][y][z][t] = 0.0; @@ -972,5 +976,5 @@ } - // Allocate output vector + // Allocoutput vector tmp = psVectorAlloc(vecLen, PS_TYPE_F64); Index: trunk/psLib/src/math/psPolynomial.h =================================================================== --- trunk/psLib/src/math/psPolynomial.h (revision 5074) +++ trunk/psLib/src/math/psPolynomial.h (revision 5090) @@ -11,6 +11,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.56 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-19 19:53:13 $ + * @version $Revision: 1.57 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-22 02:47:16 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -64,56 +64,68 @@ */ typedef enum { - PS_POLYNOMIAL_ORD, ///< Ordinary Polynomial - PS_POLYNOMIAL_CHEB ///< Chebyshev Polynomial + PS_POLYNOMIAL_ORD, ///< Ordinary Polynomial + PS_POLYNOMIAL_CHEB ///< Chebyshev Polynomial } psPolynomialType; +// XXX: These are incorrect names for the order of the polynomial. We +// keep them here temporarily so we can later sed replace them with the +// correct names. /** One-dimensional polynomial */ typedef struct { - psPolynomialType type; ///< Polynomial type - unsigned int n; ///< Number of terms - psF64 *coeff; ///< Coefficients - psF64 *coeffErr; ///< Error in coefficients - psMaskType *mask; ///< Coefficient mask + psPolynomialType type; ///< Polynomial type + unsigned int COOL_1D_n; ///< Polynomial order + psF64 *coeff; ///< Coefficients + psF64 *coeffErr; ///< Error in coefficients + psMaskType *mask; ///< Coefficient mask } psPolynomial1D; +// XXX: These are incorrect names for the order of the polynomial. We +// keep them here temporarily so we can later sed replace them with the +// correct names. /** Two-dimensional polynomial */ typedef struct { - psPolynomialType type; ///< Polynomial type - unsigned int nX; ///< Number of terms in x - unsigned int nY; ///< Number of terms in y - psF64 **coeff; ///< Coefficients - psF64 **coeffErr; ///< Error in coefficients - psMaskType **mask; ///< Coefficients mask + psPolynomialType type; ///< Polynomial type + unsigned int COOL_2D_nX; ///< Polynomial order in x + unsigned int COOL_2D_nY; ///< Polynomial order in y + psF64 **coeff; ///< Coefficients + psF64 **coeffErr; ///< Error in coefficients + psMaskType **mask; ///< Coefficients mask } psPolynomial2D; +// XXX: These are incorrect names for the order of the polynomial. We +// keep them here temporarily so we can later sed replace them with the +// correct names. /** Three-dimensional polynomial */ typedef struct { - psPolynomialType type; ///< Polynomial type - unsigned int nX; ///< Number of terms in x - unsigned int nY; ///< Number of terms in y - unsigned int nZ; ///< Number of terms in z - psF64 ***coeff; ///< Coefficients - psF64 ***coeffErr; ///< Error in coefficients - psMaskType ***mask; ///< Coefficients mask + psPolynomialType type; ///< Polynomial type + unsigned int COOL_3D_nX; ///< Polynomial order in x + unsigned int COOL_3D_nY; ///< Polynomial order in y + unsigned int COOL_3D_nZ; ///< Polynomial order in z + psF64 ***coeff; ///< Coefficients + psF64 ***coeffErr; ///< Error in coefficients + psMaskType ***mask; ///< Coefficients mask } psPolynomial3D; +// XXX: These are incorrect names for the order of the polynomial. We +// keep them here temporarily so we can later sed replace them with the +// correct names. /** Four-dimensional polynomial */ typedef struct { - psPolynomialType type; ///< Polynomial type - unsigned int nX; ///< Number of terms in x - unsigned int nY; ///< Number of terms in y - unsigned int nZ; ///< Number of terms in z - unsigned int nT; ///< Number of terms in t - psF64 ****coeff; ///< Coefficients - psF64 ****coeffErr; ///< Error in coefficients - psMaskType ****mask; ///< Coefficients mask + psPolynomialType type; ///< Polynomial type + unsigned int COOL_4D_nX; ///< Polynomial order in x + unsigned int COOL_4D_nY; ///< Polynomial order in y + unsigned int COOL_4D_nZ; ///< Polynomial order in z + unsigned int COOL_4D_nT; ///< Polynomial order in t + psF64 ****coeff; ///< Coefficients + psF64 ****coeffErr; ///< Error in coefficients + psMaskType ****mask; ///< Coefficients mask } psPolynomial4D; Index: trunk/psLib/src/math/psSpline.c =================================================================== --- trunk/psLib/src/math/psSpline.c (revision 5074) +++ trunk/psLib/src/math/psSpline.c (revision 5090) @@ -7,6 +7,6 @@ * splines. * -* @version $Revision: 1.126 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-09-20 02:43:53 $ +* @version $Revision: 1.127 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-09-22 02:47:16 $ * * @@ -383,5 +383,5 @@ // If these are linear splines, which means their polynomials will have // two coefficients, then we do the simple calculation. - if (2 == (spline->spline[0])->n) { + if (1 == (spline->spline[0])->COOL_1D_n) { for (i=0;in;i++) { slope = (y32->data.F32[i+1] - y32->data.F32[i]) / @@ -402,8 +402,8 @@ // Check if these are cubic splines (n==4). If not, psError. - if (4 != (spline->spline[0])->n) { + if (3 != (spline->spline[0])->COOL_1D_n) { psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Don't know how to generate %u-order splines.", - (spline->spline[0])->n-1); + (spline->spline[0])->COOL_1D_n); return(NULL); } @@ -745,5 +745,5 @@ tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); for (unsigned int i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); + (tmpSpline->spline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); } @@ -809,5 +809,5 @@ tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); for (unsigned int i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); + (tmpSpline->spline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); } Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 5074) +++ trunk/psLib/src/math/psStats.c (revision 5090) @@ -17,6 +17,6 @@ * * - * @version $Revision: 1.145 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-11 22:18:40 $ + * @version $Revision: 1.146 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-22 02:47:16 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -1502,4 +1502,5 @@ psVector *y = psVectorAlloc(3, PS_TYPE_F64); psVector *yErr = psVectorAlloc(3, PS_TYPE_F64); + // XXX: Why was this 2 when the alloc function specified number of terms? Note: it's correct now. psPolynomial1D *myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); @@ -1552,5 +1553,4 @@ // Determine the coefficients of the polynomial. - // myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, yErr, x); if (myPoly == NULL) { @@ -1827,5 +1827,5 @@ } - psPolynomial1D *tmpPoly = psPolynomial1DAlloc(3, PS_POLYNOMIAL_ORD); + psPolynomial1D *tmpPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); // XXX: What about the NULL x argument? tmpPoly = psVectorFitPolynomial1D(tmpPoly, NULL, 0, y, NULL, NULL);