Index: trunk/psLib/src/math/psPolynomial.c =================================================================== --- trunk/psLib/src/math/psPolynomial.c (revision 1406) +++ trunk/psLib/src/math/psPolynomial.c (revision 1407) @@ -1,2 +1,3 @@ + /** @file psFunctions.c * @@ -7,11 +8,14 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-08-06 22:34:05 $ + * @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-08-07 00:06:06 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ -/*****************************************************************************/ + +/*****************************************************************************/ + /* INCLUDE FILES */ + /*****************************************************************************/ #include @@ -31,6 +35,9 @@ #include #include -/*****************************************************************************/ + +/*****************************************************************************/ + /* DEFINE STATEMENTS */ + /*****************************************************************************/ @@ -38,17 +45,21 @@ /*****************************************************************************/ + /* TYPE DEFINITIONS */ -/*****************************************************************************/ -static void polynomial1DFree(psPolynomial1D *myPoly); -static void polynomial2DFree(psPolynomial2D *myPoly); -static void polynomial3DFree(psPolynomial3D *myPoly); -static void polynomial4DFree(psPolynomial4D *myPoly); -static void dPolynomial1DFree(psDPolynomial1D *myPoly); -static void dPolynomial2DFree(psDPolynomial2D *myPoly); -static void dPolynomial3DFree(psDPolynomial3D *myPoly); -static void dPolynomial4DFree(psDPolynomial4D *myPoly); - -/*****************************************************************************/ + +/*****************************************************************************/ +static void polynomial1DFree(psPolynomial1D * myPoly); +static void polynomial2DFree(psPolynomial2D * myPoly); +static void polynomial3DFree(psPolynomial3D * myPoly); +static void polynomial4DFree(psPolynomial4D * myPoly); +static void dPolynomial1DFree(psDPolynomial1D * myPoly); +static void dPolynomial2DFree(psDPolynomial2D * myPoly); +static void dPolynomial3DFree(psDPolynomial3D * myPoly); +static void dPolynomial4DFree(psDPolynomial4D * myPoly); + +/*****************************************************************************/ + /* GLOBAL VARIABLES */ + /*****************************************************************************/ @@ -56,5 +67,7 @@ /*****************************************************************************/ + /* FILE STATIC VARIABLES */ + /*****************************************************************************/ @@ -62,5 +75,7 @@ /*****************************************************************************/ + /* FUNCTION IMPLEMENTATION - LOCAL */ + /*****************************************************************************/ @@ -68,5 +83,7 @@ /*****************************************************************************/ + /* FUNCTION IMPLEMENTATION - PUBLIC */ + /*****************************************************************************/ @@ -76,9 +93,5 @@ evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] *****************************************************************************/ -float -psGaussian(float x, - float mean, - float sigma, - bool normal) +float psGaussian(float x, float mean, float sigma, bool normal) { float tmp = 1.0; @@ -94,5 +107,5 @@ } - return(tmp * exp(-((x-mean) * (x-mean)) / (2.0 * sigma * sigma))); + return (tmp * exp(-((x - mean) * (x - mean)) / (2.0 * sigma * sigma))); } @@ -107,7 +120,5 @@ NOTE: There is no way to seed the random generator. *****************************************************************************/ -psVector *psGaussianDev(float mean, - float sigma, - int Npts) +psVector *psGaussianDev(float mean, float sigma, int Npts) { psVector *gauss = NULL; @@ -127,7 +138,6 @@ // NOTE: Should I free r as well? - return(gauss); -} - + return (gauss); +} /***************************************************************************** @@ -140,10 +150,10 @@ newPoly = (psPolynomial1D *) psAlloc(sizeof(psPolynomial1D)); - p_psMemSetDeallocator(newPoly,(psFreeFcn)polynomial1DFree); + p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial1DFree); newPoly->n = n; - newPoly->coeff = (float *) psAlloc(n * sizeof(float)); - newPoly->coeffErr = (float *) psAlloc(n * sizeof(float)); - newPoly->mask = (char *) psAlloc(n * sizeof(char)); - for (i=0;icoeff = (float *)psAlloc(n * sizeof(float)); + newPoly->coeffErr = (float *)psAlloc(n * sizeof(float)); + newPoly->mask = (char *)psAlloc(n * sizeof(char)); + for (i = 0; i < n; i++) { newPoly->coeff[i] = 0.0; newPoly->coeffErr[i] = 0.0; @@ -151,5 +161,5 @@ } - return(newPoly); + return (newPoly); } @@ -161,18 +171,18 @@ newPoly = (psPolynomial2D *) psAlloc(sizeof(psPolynomial2D)); - p_psMemSetDeallocator(newPoly,(psFreeFcn)polynomial2DFree); + p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial2DFree); newPoly->nX = nX; newPoly->nY = nY; - newPoly->coeff = (float **) psAlloc(nX * sizeof(float *)); - newPoly->coeffErr = (float **) psAlloc(nX * sizeof(float *)); - newPoly->mask = (char **) psAlloc(nX * sizeof(char *)); - for (x=0;xcoeff[x] = (float *) psAlloc(nY * sizeof(float)); - newPoly->coeffErr[x] = (float *) psAlloc(nY * sizeof(float)); - newPoly->mask[x] = (char *) psAlloc(nY * sizeof(char)); - } - for (x=0;xcoeff = (float **)psAlloc(nX * sizeof(float *)); + newPoly->coeffErr = (float **)psAlloc(nX * sizeof(float *)); + newPoly->mask = (char **)psAlloc(nX * sizeof(char *)); + for (x = 0; x < nX; x++) { + newPoly->coeff[x] = (float *)psAlloc(nY * sizeof(float)); + newPoly->coeffErr[x] = (float *)psAlloc(nY * sizeof(float)); + newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char)); + } + for (x = 0; x < nX; x++) { + for (y = 0; y < nY; y++) { newPoly->coeff[x][y] = 0.0; newPoly->coeffErr[x][y] = 0.0; @@ -181,5 +191,5 @@ } - return(newPoly); + return (newPoly); } @@ -192,25 +202,25 @@ newPoly = (psPolynomial3D *) psAlloc(sizeof(psPolynomial3D)); - p_psMemSetDeallocator(newPoly,(psFreeFcn)polynomial3DFree); + p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial3DFree); newPoly->nX = nX; newPoly->nY = nY; newPoly->nZ = nZ; - newPoly->coeff = (float ***) psAlloc(nX * sizeof(float **)); - newPoly->coeffErr = (float ***) psAlloc(nX * sizeof(float **)); - newPoly->mask = (char ***) psAlloc(nX * sizeof(char **)); - for (x=0;xcoeff[x] = (float **) psAlloc(nY * sizeof(float *)); - newPoly->coeffErr[x] = (float **) psAlloc(nY * sizeof(float *)); - newPoly->mask[x] = (char **) psAlloc(nY * sizeof(char *)); - for (y=0;ycoeff[x][y] = (float *) psAlloc(nZ * sizeof(float)); - newPoly->coeffErr[x][y] = (float *) psAlloc(nZ * sizeof(float)); - newPoly->mask[x][y] = (char *) psAlloc(nZ * sizeof(char)); - } - } - for (x=0;xcoeff = (float ***)psAlloc(nX * sizeof(float **)); + newPoly->coeffErr = (float ***)psAlloc(nX * sizeof(float **)); + newPoly->mask = (char ***)psAlloc(nX * sizeof(char **)); + for (x = 0; x < nX; x++) { + newPoly->coeff[x] = (float **)psAlloc(nY * sizeof(float *)); + newPoly->coeffErr[x] = (float **)psAlloc(nY * sizeof(float *)); + newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *)); + for (y = 0; y < nY; y++) { + newPoly->coeff[x][y] = (float *)psAlloc(nZ * sizeof(float)); + newPoly->coeffErr[x][y] = (float *)psAlloc(nZ * sizeof(float)); + newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char)); + } + } + for (x = 0; x < nX; x++) { + for (y = 0; y < nY; y++) { + for (z = 0; z < nZ; z++) { newPoly->coeff[x][y][z] = 0.0; newPoly->coeffErr[x][y][z] = 0.0; @@ -220,5 +230,5 @@ } - return(newPoly); + return (newPoly); } @@ -232,5 +242,5 @@ newPoly = (psPolynomial4D *) psAlloc(sizeof(psPolynomial4D)); - p_psMemSetDeallocator(newPoly,(psFreeFcn)polynomial4DFree); + p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial4DFree); newPoly->nW = nW; newPoly->nX = nX; @@ -238,26 +248,26 @@ newPoly->nZ = nZ; - newPoly->coeff = (float ****) psAlloc(nW * sizeof(float ***)); - newPoly->coeffErr = (float ****) psAlloc(nW * sizeof(float ***)); - newPoly->mask = (char ****) psAlloc(nW * sizeof(char ***)); - for (w=0;wcoeff[w] = (float ***) psAlloc(nX * sizeof(float **)); - newPoly->coeffErr[w] = (float ***) psAlloc(nX * sizeof(float **)); - newPoly->mask[w] = (char ***) psAlloc(nX * sizeof(char **)); - for (x=0;xcoeff[w][x] = (float **) psAlloc(nY * sizeof(float *)); - newPoly->coeffErr[w][x] = (float **) psAlloc(nY * sizeof(float *)); - newPoly->mask[w][x] = (char **) psAlloc(nY * sizeof(char *)); - for (y=0;ycoeff[w][x][y] = (float *) psAlloc(nZ * sizeof(float)); - newPoly->coeffErr[w][x][y] = (float *) psAlloc(nZ * sizeof(float)); - newPoly->mask[w][x][y] = (char *) psAlloc(nZ * sizeof(char)); - } - } - } - for (w=0;wcoeff = (float ****)psAlloc(nW * sizeof(float ***)); + newPoly->coeffErr = (float ****)psAlloc(nW * sizeof(float ***)); + newPoly->mask = (char ****)psAlloc(nW * sizeof(char ***)); + for (w = 0; w < nW; w++) { + newPoly->coeff[w] = (float ***)psAlloc(nX * sizeof(float **)); + newPoly->coeffErr[w] = (float ***)psAlloc(nX * sizeof(float **)); + newPoly->mask[w] = (char ***)psAlloc(nX * sizeof(char **)); + for (x = 0; x < nX; x++) { + newPoly->coeff[w][x] = (float **)psAlloc(nY * sizeof(float *)); + newPoly->coeffErr[w][x] = (float **)psAlloc(nY * sizeof(float *)); + newPoly->mask[w][x] = (char **)psAlloc(nY * sizeof(char *)); + for (y = 0; y < nY; y++) { + newPoly->coeff[w][x][y] = (float *)psAlloc(nZ * sizeof(float)); + newPoly->coeffErr[w][x][y] = (float *)psAlloc(nZ * sizeof(float)); + newPoly->mask[w][x][y] = (char *)psAlloc(nZ * sizeof(char)); + } + } + } + for (w = 0; w < nW; w++) { + for (x = 0; x < nX; x++) { + for (y = 0; y < nY; y++) { + for (z = 0; z < nZ; z++) { newPoly->coeff[w][x][y][z] = 0.0; newPoly->coeffErr[w][x][y][z] = 0.0; @@ -268,8 +278,8 @@ } - return(newPoly); -} - -static void polynomial1DFree(psPolynomial1D *myPoly) + return (newPoly); +} + +static void polynomial1DFree(psPolynomial1D * myPoly) { psFree(myPoly->coeff); @@ -278,9 +288,9 @@ } -static void polynomial2DFree(psPolynomial2D *myPoly) -{ - int x = 0; - - for (x=0;xnX;x++) { +static void polynomial2DFree(psPolynomial2D * myPoly) +{ + int x = 0; + + for (x = 0; x < myPoly->nX; x++) { psFree(myPoly->coeff[x]); psFree(myPoly->coeffErr[x]); @@ -292,11 +302,11 @@ } -static void polynomial3DFree(psPolynomial3D *myPoly) -{ - int x = 0; - int y = 0; - - for (x=0;xnX;x++) { - for (y=0;ynY;y++) { +static void polynomial3DFree(psPolynomial3D * myPoly) +{ + int x = 0; + int y = 0; + + for (x = 0; x < myPoly->nX; x++) { + for (y = 0; y < myPoly->nY; y++) { psFree(myPoly->coeff[x][y]); psFree(myPoly->coeffErr[x][y]); @@ -313,5 +323,5 @@ } -static void polynomial4DFree(psPolynomial4D *myPoly) +static void polynomial4DFree(psPolynomial4D * myPoly) { int w = 0; @@ -319,7 +329,7 @@ int y = 0; - for (w=0;wnW;w++) { - for (x=0;xnX;x++) { - for (y=0;ynY;y++) { + for (w = 0; w < myPoly->nW; w++) { + for (x = 0; x < myPoly->nX; x++) { + for (y = 0; y < myPoly->nY; y++) { psFree(myPoly->coeff[w][x][y]); psFree(myPoly->coeffErr[w][x][y]); @@ -343,7 +353,5 @@ Polynomial coefficients will be accessed in [w][x][y][z] fashion. *****************************************************************************/ -float -psPolynomial1DEval(float x, - const psPolynomial1D *myPoly) +float psPolynomial1DEval(float x, const psPolynomial1D * myPoly) { int loop_x = 0; @@ -355,8 +363,7 @@ } - // NOTE: Do we want to flag this case? if (myPoly->n == 0) { - return(1.0); + return (1.0); } @@ -365,17 +372,13 @@ } - for (loop_x=0;loop_xn;loop_x++) { - polySum+= xSum * myPoly->coeff[loop_x]; - xSum*=x; - } - - return(polySum); -} - - -float -psPolynomial2DEval(float x, - float y, - const psPolynomial2D *myPoly) + for (loop_x = 0; loop_x < myPoly->n; loop_x++) { + polySum += xSum * myPoly->coeff[loop_x]; + xSum *= x; + } + + return (polySum); +} + +float psPolynomial2DEval(float x, float y, const psPolynomial2D * myPoly) { int loop_x = 0; @@ -385,21 +388,17 @@ float ySum = 1.0; - for (loop_x=0;loop_xnX;loop_x++) { + for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { ySum = xSum; - for (loop_y=0;loop_ynY;loop_y++) { - polySum+= ySum * myPoly->coeff[loop_x][loop_y]; - ySum*=y; - } - xSum*=x; - } - - return(polySum); -} - -float -psPolynomial3DEval(float x, - float y, - float z, - const psPolynomial3D *myPoly) + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { + polySum += ySum * myPoly->coeff[loop_x][loop_y]; + ySum *= y; + } + xSum *= x; + } + + return (polySum); +} + +float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D * myPoly) { int loop_x = 0; @@ -411,26 +410,21 @@ float zSum = 1.0; - for (loop_x=0;loop_xnX;loop_x++) { + for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { ySum = xSum; - for (loop_y=0;loop_ynY;loop_y++) { + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { zSum = ySum; - for (loop_z=0;loop_znZ;loop_z++) { - polySum+= zSum * myPoly->coeff[loop_x][loop_y][loop_z]; - zSum*=z; - } - ySum*=y; - } - xSum*=x; - } - - return(polySum); -} - -float -psPolynomial4DEval(float w, - float x, - float y, - float z, - const psPolynomial4D *myPoly) + for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { + polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z]; + zSum *= z; + } + ySum *= y; + } + xSum *= x; + } + + return (polySum); +} + +float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D * myPoly) { int loop_w = 0; @@ -444,27 +438,23 @@ float zSum = 1.0; - for (loop_w=0;loop_wnW;loop_w++) { + for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { xSum = wSum; - for (loop_x=0;loop_xnX;loop_x++) { + for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { ySum = xSum; - for (loop_y=0;loop_ynY;loop_y++) { + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { zSum = ySum; - for (loop_z=0;loop_znZ;loop_z++) { - polySum+= zSum * - myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; - zSum*=z; + for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { + polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; + zSum *= z; } - ySum*=y; - } - xSum*=x; - } - wSum*=w; - } - - return(polySum); -} - - - + ySum *= y; + } + xSum *= x; + } + wSum *= w; + } + + return (polySum); +} psDPolynomial1D *psDPolynomial1DAlloc(int n) @@ -474,10 +464,10 @@ newPoly = (psDPolynomial1D *) psAlloc(sizeof(psDPolynomial1D)); - p_psMemSetDeallocator(newPoly,(psFreeFcn)dPolynomial1DFree); + p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial1DFree); newPoly->n = n; - newPoly->coeff = (double *) psAlloc(n * sizeof(double)); - newPoly->coeffErr = (double *) psAlloc(n * sizeof(double)); - newPoly->mask = (char *) psAlloc(n * sizeof(char)); - for (i=0;icoeff = (double *)psAlloc(n * sizeof(double)); + newPoly->coeffErr = (double *)psAlloc(n * sizeof(double)); + newPoly->mask = (char *)psAlloc(n * sizeof(char)); + for (i = 0; i < n; i++) { newPoly->coeff[i] = 0.0; newPoly->coeffErr[i] = 0.0; @@ -485,5 +475,5 @@ } - return(newPoly); + return (newPoly); } @@ -495,18 +485,18 @@ newPoly = (psDPolynomial2D *) psAlloc(sizeof(psDPolynomial2D)); - p_psMemSetDeallocator(newPoly,(psFreeFcn)dPolynomial2DFree); + p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial2DFree); newPoly->nX = nX; newPoly->nY = nY; - newPoly->coeff = (double **) psAlloc(nX * sizeof(double *)); - newPoly->coeffErr = (double **) psAlloc(nX * sizeof(double *)); - newPoly->mask = (char **) psAlloc(nX * sizeof(char *)); - for (x=0;xcoeff[x] = (double *) psAlloc(nY * sizeof(double)); - newPoly->coeffErr[x] = (double *) psAlloc(nY * sizeof(double)); - newPoly->mask[x] = (char *) psAlloc(nY * sizeof(char)); - } - for (x=0;xcoeff = (double **)psAlloc(nX * sizeof(double *)); + newPoly->coeffErr = (double **)psAlloc(nX * sizeof(double *)); + newPoly->mask = (char **)psAlloc(nX * sizeof(char *)); + for (x = 0; x < nX; x++) { + newPoly->coeff[x] = (double *)psAlloc(nY * sizeof(double)); + newPoly->coeffErr[x] = (double *)psAlloc(nY * sizeof(double)); + newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char)); + } + for (x = 0; x < nX; x++) { + for (y = 0; y < nY; y++) { newPoly->coeff[x][y] = 0.0; newPoly->coeffErr[x][y] = 0.0; @@ -515,5 +505,5 @@ } - return(newPoly); + return (newPoly); } @@ -526,25 +516,25 @@ newPoly = (psDPolynomial3D *) psAlloc(sizeof(psDPolynomial3D)); - p_psMemSetDeallocator(newPoly,(psFreeFcn)dPolynomial3DFree); + p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial3DFree); newPoly->nX = nX; newPoly->nY = nY; newPoly->nZ = nZ; - newPoly->coeff = (double ***) psAlloc(nX * sizeof(double **)); - newPoly->coeffErr = (double ***) psAlloc(nX * sizeof(double **)); - newPoly->mask = (char ***) psAlloc(nX * sizeof(char **)); - for (x=0;xcoeff[x] = (double **) psAlloc(nY * sizeof(double *)); - newPoly->coeffErr[x] = (double **) psAlloc(nY * sizeof(double *)); - newPoly->mask[x] = (char **) psAlloc(nY * sizeof(char *)); - for (y=0;ycoeff[x][y] = (double *) psAlloc(nZ * sizeof(double)); - newPoly->coeffErr[x][y] = (double *) psAlloc(nZ * sizeof(double)); - newPoly->mask[x][y] = (char *) psAlloc(nZ * sizeof(char)); - } - } - for (x=0;xcoeff = (double ***)psAlloc(nX * sizeof(double **)); + newPoly->coeffErr = (double ***)psAlloc(nX * sizeof(double **)); + newPoly->mask = (char ***)psAlloc(nX * sizeof(char **)); + for (x = 0; x < nX; x++) { + newPoly->coeff[x] = (double **)psAlloc(nY * sizeof(double *)); + newPoly->coeffErr[x] = (double **)psAlloc(nY * sizeof(double *)); + newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *)); + for (y = 0; y < nY; y++) { + newPoly->coeff[x][y] = (double *)psAlloc(nZ * sizeof(double)); + newPoly->coeffErr[x][y] = (double *)psAlloc(nZ * sizeof(double)); + newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char)); + } + } + for (x = 0; x < nX; x++) { + for (y = 0; y < nY; y++) { + for (z = 0; z < nZ; z++) { newPoly->coeff[x][y][z] = 0.0; newPoly->coeffErr[x][y][z] = 0.0; @@ -554,5 +544,5 @@ } - return(newPoly); + return (newPoly); } @@ -566,5 +556,5 @@ newPoly = (psDPolynomial4D *) psAlloc(sizeof(psDPolynomial4D)); - p_psMemSetDeallocator(newPoly,(psFreeFcn)dPolynomial4DFree); + p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial4DFree); newPoly->nW = nW; newPoly->nX = nX; @@ -572,26 +562,26 @@ newPoly->nZ = nZ; - newPoly->coeff = (double ****) psAlloc(nW * sizeof(double ***)); - newPoly->coeffErr = (double ****) psAlloc(nW * sizeof(double ***)); - newPoly->mask = (char ****) psAlloc(nW * sizeof(char ***)); - for (w=0;wcoeff[w] = (double ***) psAlloc(nX * sizeof(double **)); - newPoly->coeffErr[w] = (double ***) psAlloc(nX * sizeof(double **)); - newPoly->mask[w] = (char ***) psAlloc(nX * sizeof(char **)); - for (x=0;xcoeff[w][x] = (double **) psAlloc(nY * sizeof(double *)); - newPoly->coeffErr[w][x] = (double **) psAlloc(nY * sizeof(double *)); - newPoly->mask[w][x] = (char **) psAlloc(nY * sizeof(char *)); - for (y=0;ycoeff[w][x][y] = (double *) psAlloc(nZ * sizeof(double)); - newPoly->coeffErr[w][x][y] = (double *) psAlloc(nZ * sizeof(double)); - newPoly->mask[w][x][y] = (char *) psAlloc(nZ * sizeof(char)); - } - } - } - for (w=0;wcoeff = (double ****)psAlloc(nW * sizeof(double ***)); + newPoly->coeffErr = (double ****)psAlloc(nW * sizeof(double ***)); + newPoly->mask = (char ****)psAlloc(nW * sizeof(char ***)); + for (w = 0; w < nW; w++) { + newPoly->coeff[w] = (double ***)psAlloc(nX * sizeof(double **)); + newPoly->coeffErr[w] = (double ***)psAlloc(nX * sizeof(double **)); + newPoly->mask[w] = (char ***)psAlloc(nX * sizeof(char **)); + for (x = 0; x < nX; x++) { + newPoly->coeff[w][x] = (double **)psAlloc(nY * sizeof(double *)); + newPoly->coeffErr[w][x] = (double **)psAlloc(nY * sizeof(double *)); + newPoly->mask[w][x] = (char **)psAlloc(nY * sizeof(char *)); + for (y = 0; y < nY; y++) { + newPoly->coeff[w][x][y] = (double *)psAlloc(nZ * sizeof(double)); + newPoly->coeffErr[w][x][y] = (double *)psAlloc(nZ * sizeof(double)); + newPoly->mask[w][x][y] = (char *)psAlloc(nZ * sizeof(char)); + } + } + } + for (w = 0; w < nW; w++) { + for (x = 0; x < nX; x++) { + for (y = 0; y < nY; y++) { + for (z = 0; z < nZ; z++) { newPoly->coeff[w][x][y][z] = 0.0; newPoly->coeffErr[w][x][y][z] = 0.0; @@ -602,8 +592,8 @@ } - return(newPoly); -} - -static void dPolynomial1DFree(psDPolynomial1D *myPoly) + return (newPoly); +} + +static void dPolynomial1DFree(psDPolynomial1D * myPoly) { psFree(myPoly->coeff); @@ -612,9 +602,9 @@ } -static void dPolynomial2DFree(psDPolynomial2D *myPoly) -{ - int x = 0; - - for (x=0;xnX;x++) { +static void dPolynomial2DFree(psDPolynomial2D * myPoly) +{ + int x = 0; + + for (x = 0; x < myPoly->nX; x++) { psFree(myPoly->coeff[x]); psFree(myPoly->coeffErr[x]); @@ -626,11 +616,11 @@ } -static void dPolynomial3DFree(psDPolynomial3D *myPoly) -{ - int x = 0; - int y = 0; - - for (x=0;xnX;x++) { - for (y=0;ynY;y++) { +static void dPolynomial3DFree(psDPolynomial3D * myPoly) +{ + int x = 0; + int y = 0; + + for (x = 0; x < myPoly->nX; x++) { + for (y = 0; y < myPoly->nY; y++) { psFree(myPoly->coeff[x][y]); psFree(myPoly->coeffErr[x][y]); @@ -647,5 +637,5 @@ } -static void dPolynomial4DFree(psDPolynomial4D *myPoly) +static void dPolynomial4DFree(psDPolynomial4D * myPoly) { int w = 0; @@ -653,7 +643,7 @@ int y = 0; - for (w=0;wnW;w++) { - for (x=0;xnX;x++) { - for (y=0;ynY;y++) { + for (w = 0; w < myPoly->nW; w++) { + for (x = 0; x < myPoly->nX; x++) { + for (y = 0; y < myPoly->nY; y++) { psFree(myPoly->coeff[w][x][y]); psFree(myPoly->coeffErr[w][x][y]); @@ -677,7 +667,5 @@ Polynomial coefficients will be accessed in [w][x][y][z] fashion. *****************************************************************************/ -double -psDPolynomial1DEval(double x, - const psDPolynomial1D *myPoly) +double psDPolynomial1DEval(double x, const psDPolynomial1D * myPoly) { int loop_x = 0; @@ -687,20 +675,16 @@ // NOTE: Do we want to flag this case? if (myPoly->n == 0) { - return(1.0); - } - - for (loop_x=0;loop_xn;loop_x++) { - polySum+= xSum * myPoly->coeff[loop_x]; - xSum*=x; - } - - return(polySum); -} - - -double -psDPolynomial2DEval(double x, - double y, - const psDPolynomial2D *myPoly) + return (1.0); + } + + for (loop_x = 0; loop_x < myPoly->n; loop_x++) { + polySum += xSum * myPoly->coeff[loop_x]; + xSum *= x; + } + + return (polySum); +} + +double psDPolynomial2DEval(double x, double y, const psDPolynomial2D * myPoly) { int loop_x = 0; @@ -710,21 +694,17 @@ double ySum = 1.0; - for (loop_x=0;loop_xnX;loop_x++) { + for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { ySum = xSum; - for (loop_y=0;loop_ynY;loop_y++) { - polySum+= ySum * myPoly->coeff[loop_x][loop_y]; - ySum*=y; - } - xSum*=x; - } - - return(polySum); -} - -double -psDPolynomial3DEval(double x, - double y, - double z, - const psDPolynomial3D *myPoly) + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { + polySum += ySum * myPoly->coeff[loop_x][loop_y]; + ySum *= y; + } + xSum *= x; + } + + return (polySum); +} + +double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D * myPoly) { int loop_x = 0; @@ -736,26 +716,21 @@ double zSum = 1.0; - for (loop_x=0;loop_xnX;loop_x++) { + for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { ySum = xSum; - for (loop_y=0;loop_ynY;loop_y++) { + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { zSum = ySum; - for (loop_z=0;loop_znZ;loop_z++) { - polySum+= zSum * myPoly->coeff[loop_x][loop_y][loop_z]; - zSum*=z; - } - ySum*=y; - } - xSum*=x; - } - - return(polySum); -} - -double -psDPolynomial4DEval(double w, - double x, - double y, - double z, - const psDPolynomial4D *myPoly) + for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { + polySum += zSum * myPoly->coeff[loop_x][loop_y][loop_z]; + zSum *= z; + } + ySum *= y; + } + xSum *= x; + } + + return (polySum); +} + +double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D * myPoly) { int loop_w = 0; @@ -769,22 +744,21 @@ double zSum = 1.0; - for (loop_w=0;loop_wnW;loop_w++) { + for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { xSum = wSum; - for (loop_x=0;loop_xnX;loop_x++) { + for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { ySum = xSum; - for (loop_y=0;loop_ynY;loop_y++) { + for (loop_y = 0; loop_y < myPoly->nY; loop_y++) { zSum = ySum; - for (loop_z=0;loop_znZ;loop_z++) { - polySum+= zSum * - myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; - zSum*=z; + for (loop_z = 0; loop_z < myPoly->nZ; loop_z++) { + polySum += zSum * myPoly->coeff[loop_w][loop_x][loop_y][loop_z]; + zSum *= z; } - ySum*=y; - } - xSum*=x; - } - wSum*=w; - } - - return(polySum); -} + ySum *= y; + } + xSum *= x; + } + wSum *= w; + } + + return (polySum); +}