/** @file  psFunctions.c
 *
 *  @brief Contains basic function allocation, deallocation, and evaluation
 *         routines.
 *
 *  This file will hold the functions for allocated, freeing, and evaluating
 *  polynomials.  It also contains a Gaussian functions.
 *
 *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2004-07-15 23:52:34 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 */
/*****************************************************************************/
/*  INCLUDE FILES                                                            */
/*****************************************************************************/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <stdarg.h>
#include <float.h>
#include <math.h>

#include "psMemory.h"
#include "psVector.h"
#include "psTrace.h"
#include "psError.h"
#include "psAbort.h"
#include "psFunctions.h"

#include "float.h"
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

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);

/*****************************************************************************/
/*  FUNCTION IMPLEMENTATION - PUBLIC                                         */
/*****************************************************************************/

/*****************************************************************************
    Evaluate a non-normalized Gaussian with the given mean and sigma at the
    given coordianate.  Note that this is not a Gaussian deviate.  The
    evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f]
 *****************************************************************************/
float
psGaussian(float x,
           float mean,
           float stddev,
           int normal)
{
    float tmp = 0.0;

    if (normal == 1) {
        #ifdef DARWIN
        tmp = 1.0 / (float)sqrt(2.0 * M_PI * (stddev * stddev));
        #else

        tmp = 1.0 / sqrtf(2.0 * M_PI * (stddev * stddev));
        #endif

        return(tmp * exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
    } else {
        return(exp(-((x-mean) * (x-mean)) / (2.0 * stddev * stddev)));
    }
}

/*****************************************************************************
    psGaussianDev()
 This routine creates a psVector of the specified size and type F32
 and fills it with a random Gaussian distribution of numbers with
 the specified mean and sigma.  This routine makes use of the GSL
 routines for generating both uniformly distributed numbers and the
 Gaussian distribution as well.
 
 NOTE: There is no way to seed the random generator.
 *****************************************************************************/
psVector *psGaussianDev(float mean,
                        float sigma,
                        int Npts)
{
    psVector *gauss = NULL;
    const gsl_rng_type *T = NULL;
    gsl_rng *r = NULL;
    int i = 0;

    gauss = psVectorAlloc(Npts, PS_TYPE_F32);
    gauss->n = Npts;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc(T);

    for (i = 0; i < Npts; i++) {
        gauss->data.F32[i] = mean + gsl_ran_gaussian(r, sigma);
    }

    // NOTE: Should I free r as well?
    return(gauss);
}


/*****************************************************************************
    This routine must allocate memory for the polynomial structures.
 *****************************************************************************/
psPolynomial1D *psPolynomial1DAlloc(int n)
{
    int i = 0;
    psPolynomial1D *newPoly = NULL;

    newPoly = (psPolynomial1D *) psAlloc(sizeof(psPolynomial1D));
    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;i<n;i++) {
        newPoly->coeff[i] = 0.0;
        newPoly->coeffErr[i] = 0.0;
        newPoly->mask[i] = 0.0;
    }

    return(newPoly);
}

psPolynomial2D *psPolynomial2DAlloc(int nX, int nY)
{
    int x = 0;
    int y = 0;
    psPolynomial2D *newPoly = NULL;

    newPoly = (psPolynomial2D *) psAlloc(sizeof(psPolynomial2D));
    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;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;
            newPoly->mask[x][y] = 0;
        }
    }

    return(newPoly);
}

psPolynomial3D *psPolynomial3DAlloc(int nX, int nY, int nZ)
{
    int x = 0;
    int y = 0;
    int z = 0;
    psPolynomial3D *newPoly = NULL;

    newPoly = (psPolynomial3D *) psAlloc(sizeof(psPolynomial3D));
    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;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;
                newPoly->mask[x][y][z] = 0;
            }
        }
    }

    return(newPoly);
}

psPolynomial4D *psPolynomial4DAlloc(int nW, int nX, int nY, int nZ)
{
    int w = 0;
    int x = 0;
    int y = 0;
    int z = 0;
    psPolynomial4D *newPoly = NULL;

    newPoly = (psPolynomial4D *) psAlloc(sizeof(psPolynomial4D));
    p_psMemSetDeallocator(newPoly,(psFreeFcn)polynomial4DFree);
    newPoly->nW = nW;
    newPoly->nX = nX;
    newPoly->nY = nY;
    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;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;
                    newPoly->mask[w][x][y][z] = 0;
                }
            }
        }
    }

    return(newPoly);
}

static void polynomial1DFree(psPolynomial1D *myPoly)
{
    psFree(myPoly->coeff);
    psFree(myPoly->coeffErr);
    psFree(myPoly->mask);
}

static void polynomial2DFree(psPolynomial2D *myPoly)
{
    int x = 0;

    for (x=0;x<myPoly->nX;x++) {
        psFree(myPoly->coeff[x]);
        psFree(myPoly->coeffErr[x]);
        psFree(myPoly->mask[x]);
    }
    psFree(myPoly->coeff);
    psFree(myPoly->coeffErr);
    psFree(myPoly->mask);
}

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]);
            psFree(myPoly->mask[x][y]);
        }
        psFree(myPoly->coeff[x]);
        psFree(myPoly->coeffErr[x]);
        psFree(myPoly->mask[x]);
    }

    psFree(myPoly->coeff);
    psFree(myPoly->coeffErr);
    psFree(myPoly->mask);
}

static void polynomial4DFree(psPolynomial4D *myPoly)
{
    int w = 0;
    int x = 0;
    int y = 0;

    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]);
                psFree(myPoly->mask[w][x][y]);
            }
            psFree(myPoly->coeff[w][x]);
            psFree(myPoly->coeffErr[w][x]);
            psFree(myPoly->mask[w][x]);
        }
        psFree(myPoly->coeff[w]);
        psFree(myPoly->coeffErr[w]);
        psFree(myPoly->mask[w]);
    }

    psFree(myPoly->coeff);
    psFree(myPoly->coeffErr);
    psFree(myPoly->mask);
}

/*****************************************************************************
    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
 *****************************************************************************/
float
psEvalPolynomial1D(float x,
                   const psPolynomial1D *myPoly)
{
    int loop_x = 0;
    float polySum = 0.0;
    float xSum = 1.0;

    // NOTE: Do we want to flag this case?
    if (myPoly->n == 0) {
        return(1.0);
    }

    for (loop_x=0;loop_x<myPoly->n;loop_x++) {
        polySum+= xSum * myPoly->coeff[loop_x];
        xSum*=x;
    }

    return(polySum);
}


float
psEvalPolynomial2D(float x,
                   float y,
                   const psPolynomial2D *myPoly)
{
    int loop_x = 0;
    int loop_y = 0;
    float polySum = 0.0;
    float xSum = 1.0;
    float ySum = 1.0;

    for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
        ySum = xSum;
        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
psEvalPolynomial3D(float x,
                   float y,
                   float z,
                   const psPolynomial3D *myPoly)
{
    int loop_x = 0;
    int loop_y = 0;
    int loop_z = 0;
    float polySum = 0.0;
    float xSum = 1.0;
    float ySum = 1.0;
    float zSum = 1.0;

    for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
        ySum = xSum;
        for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
            zSum = ySum;
            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
psEvalPolynomial4D(float w,
                   float x,
                   float y,
                   float z,
                   const psPolynomial4D *myPoly)
{
    int loop_w = 0;
    int loop_x = 0;
    int loop_y = 0;
    int loop_z = 0;
    float polySum = 0.0;
    float wSum = 1.0;
    float xSum = 1.0;
    float ySum = 1.0;
    float zSum = 1.0;

    for (loop_w=0;loop_w<myPoly->nW;loop_w++) {
        xSum = wSum;
        for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
            ySum = xSum;
            for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
                zSum = ySum;
                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);
}




psDPolynomial1D *psDPolynomial1DAlloc(int n)
{
    int i = 0;
    psDPolynomial1D *newPoly = NULL;

    newPoly = (psDPolynomial1D *) psAlloc(sizeof(psDPolynomial1D));
    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;i<n;i++) {
        newPoly->coeff[i] = 0.0;
        newPoly->coeffErr[i] = 0.0;
        newPoly->mask[i] = 0.0;
    }

    return(newPoly);
}

psDPolynomial2D *psDPolynomial2DAlloc(int nX, int nY)
{
    int x = 0;
    int y = 0;
    psDPolynomial2D *newPoly = NULL;

    newPoly = (psDPolynomial2D *) psAlloc(sizeof(psDPolynomial2D));
    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;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;
            newPoly->mask[x][y] = 0;
        }
    }

    return(newPoly);
}

psDPolynomial3D *psDPolynomial3DAlloc(int nX, int nY, int nZ)
{
    int x = 0;
    int y = 0;
    int z = 0;
    psDPolynomial3D *newPoly = NULL;

    newPoly = (psDPolynomial3D *) psAlloc(sizeof(psDPolynomial3D));
    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;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;
                newPoly->mask[x][y][z] = 0;
            }
        }
    }

    return(newPoly);
}

psDPolynomial4D *psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ)
{
    int w = 0;
    int x = 0;
    int y = 0;
    int z = 0;
    psDPolynomial4D *newPoly = NULL;

    newPoly = (psDPolynomial4D *) psAlloc(sizeof(psDPolynomial4D));
    p_psMemSetDeallocator(newPoly,(psFreeFcn)dPolynomial4DFree);
    newPoly->nW = nW;
    newPoly->nX = nX;
    newPoly->nY = nY;
    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;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;
                    newPoly->mask[w][x][y][z] = 0;
                }
            }
        }
    }

    return(newPoly);
}

static void dPolynomial1DFree(psDPolynomial1D *myPoly)
{
    psFree(myPoly->coeff);
    psFree(myPoly->coeffErr);
    psFree(myPoly->mask);
}

static void dPolynomial2DFree(psDPolynomial2D *myPoly)
{
    int x = 0;

    for (x=0;x<myPoly->nX;x++) {
        psFree(myPoly->coeff[x]);
        psFree(myPoly->coeffErr[x]);
        psFree(myPoly->mask[x]);
    }
    psFree(myPoly->coeff);
    psFree(myPoly->coeffErr);
    psFree(myPoly->mask);
}

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]);
            psFree(myPoly->mask[x][y]);
        }
        psFree(myPoly->coeff[x]);
        psFree(myPoly->coeffErr[x]);
        psFree(myPoly->mask[x]);
    }

    psFree(myPoly->coeff);
    psFree(myPoly->coeffErr);
    psFree(myPoly->mask);
}

static void dPolynomial4DFree(psDPolynomial4D *myPoly)
{
    int w = 0;
    int x = 0;
    int y = 0;

    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]);
                psFree(myPoly->mask[w][x][y]);
            }
            psFree(myPoly->coeff[w][x]);
            psFree(myPoly->coeffErr[w][x]);
            psFree(myPoly->mask[w][x]);
        }
        psFree(myPoly->coeff[w]);
        psFree(myPoly->coeffErr[w]);
        psFree(myPoly->mask[w]);
    }

    psFree(myPoly->coeff);
    psFree(myPoly->coeffErr);
    psFree(myPoly->mask);
}

/*****************************************************************************
    Polynomial coefficients will be accessed in [w][x][y][z] fashion.
 *****************************************************************************/
double
psDEvalPolynomial1D(double x,
                    const psDPolynomial1D *myPoly)
{
    int loop_x = 0;
    double polySum = 0.0;
    double xSum = 1.0;

    // NOTE: Do we want to flag this case?
    if (myPoly->n == 0) {
        return(1.0);
    }

    for (loop_x=0;loop_x<myPoly->n;loop_x++) {
        polySum+= xSum * myPoly->coeff[loop_x];
        xSum*=x;
    }

    return(polySum);
}


double
psDEvalPolynomial2D(double x,
                    double y,
                    const psDPolynomial2D *myPoly)
{
    int loop_x = 0;
    int loop_y = 0;
    double polySum = 0.0;
    double xSum = 1.0;
    double ySum = 1.0;

    for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
        ySum = xSum;
        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
psDEvalPolynomial3D(double x,
                    double y,
                    double z,
                    const psDPolynomial3D *myPoly)
{
    int loop_x = 0;
    int loop_y = 0;
    int loop_z = 0;
    double polySum = 0.0;
    double xSum = 1.0;
    double ySum = 1.0;
    double zSum = 1.0;

    for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
        ySum = xSum;
        for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
            zSum = ySum;
            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
psDEvalPolynomial4D(double w,
                    double x,
                    double y,
                    double z,
                    const psDPolynomial4D *myPoly)
{
    int loop_w = 0;
    int loop_x = 0;
    int loop_y = 0;
    int loop_z = 0;
    double polySum = 0.0;
    double wSum = 1.0;
    double xSum = 1.0;
    double ySum = 1.0;
    double zSum = 1.0;

    for (loop_w=0;loop_w<myPoly->nW;loop_w++) {
        xSum = wSum;
        for (loop_x=0;loop_x<myPoly->nX;loop_x++) {
            ySum = xSum;
            for (loop_y=0;loop_y<myPoly->nY;loop_y++) {
                zSum = ySum;
                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);
}
