/** @file psFunctions.h
*  \brief Standard Mathematical Functions.
*  \ingroup Stats
*
*  This file will hold the prototypes for procedures which allocate, free,
*  and evaluate various polynomials.  Those polynomial structures are also
*  defined here.
*
*  @ingroup Stats
*
*  @author Someone at IfA
*  @author George Gusciora, MHPCC
*
*  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
*  @date $Date: 2004-10-06 00:48:15 $
*
*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
*/

#if !defined(PS_FUNCTIONS_H)
#define PS_FUNCTIONS_H

#include <stdio.h>
#include <stdbool.h>
#include <float.h>
#include <math.h>

#include "psVector.h"
#include "psScalar.h"

/** \addtogroup Stats
 *  \{
 */

/** 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] 
 *
 *  @return float      value on the gaussian curve given the input parameters
 */
float psGaussian(
    float x,                           ///< Value at which to evaluate
    float mean,                        ///< Mean for the Gaussian
    float stddev,                      ///< Standard deviation for the Gaussian
    bool normal                        ///< Indicates whether result should be normalized
);

/** Produce a vector of random numbers from a Gaussian distribution with
 *  the specified mean and sigma 
 *  
 *  @return psVector*    vector of random numbers
 *  
 */
psVector* psGaussianDev(
    float mean,                        ///< The mean of the Gaussian
    float sigma,                       ///< The sigma of the Gaussian
    int Npts                           ///< The size of the vector
);

typedef enum {
    PS_POLYNOMIAL_ORD,                 ///< Ordinary Polynomial
    PS_POLYNOMIAL_CHEB                 ///< Chebyshev Polynomial
} psPolynomialType;

/** One-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    int n;                             ///< Number of terms
    float *coeff;                      ///< Coefficients
    float *coeffErr;                   ///< Error in coefficients
    char *mask;                        ///< Coefficient mask
}
psPolynomial1D;

/** Two-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    int nX;                            ///< Number of terms in x
    int nY;                            ///< Number of terms in y
    float **coeff;                     ///< Coefficients
    float **coeffErr;                  ///< Error in coefficients
    char **mask;                       ///< Coefficients mask
}
psPolynomial2D;

/** Three-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    int nX;                            ///< Number of terms in x
    int nY;                            ///< Number of terms in y
    int nZ;                            ///< Number of terms in z
    float ***coeff;                    ///< Coefficients
    float ***coeffErr;                 ///< Error in coefficients
    char ***mask;                      ///< Coefficients mask
}
psPolynomial3D;

/** Four-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    int nW;                            ///< Number of terms in w
    int nX;                            ///< Number of terms in x
    int nY;                            ///< Number of terms in y
    int nZ;                            ///< Number of terms in z
    float ****coeff;                   ///< Coefficients
    float ****coeffErr;                ///< Error in coefficients
    char ****mask;                     ///< Coefficients mask
}
psPolynomial4D;


/** Allocates a psPolynomial1D structure with n terms
 *
 *  @return  psPolynomial1D*    new 1-D polynomial struct
 */
psPolynomial1D* psPolynomial1DAlloc(
    int n,                              ///< Number of terms
    psPolynomialType type               ///< Polynomial Type
);

/** Allocates a 2-D polynomial structure
 *
 *  @return  psPolynomial2D*    new 2-D polynomial struct
 */
psPolynomial2D* psPolynomial2DAlloc(
    int nX,                            ///< Number of terms in x
    int nY,                            ///< Number of terms in y
    psPolynomialType type              ///< Polynomial Type
);

/** Allocates a 3-D polynomial structure
 *
 *  @return  psPolynomial3D*    new 3-D polynomial struct
 */
psPolynomial3D* psPolynomial3DAlloc(
    int nX,                            ///< Number of terms in x
    int nY,                            ///< Number of terms in y
    int nZ,                            ///< Number of terms in z
    psPolynomialType type              ///< Polynomial Type
);

/** Allocates a 4-D polynomial structure
 *
 *  @return  psPolynomial4D*    new 4-D polynomial struct
 */
psPolynomial4D* psPolynomial4DAlloc(
    int nW,                            ///< Number of terms in w
    int nX,                            ///< Number of terms in x
    int nY,                            ///< Number of terms in y
    int nZ,                            ///< Number of terms in z
    psPolynomialType type              ///< Polynomial Type
);

/** Evaluates a 1-D polynomial at specific coordinates.
 *
 *  @return float    result of polynomial at given location
 */
float psPolynomial1DEval(
    float x,                           ///< location at which to evaluate
    const psPolynomial1D* myPoly       ///< Coefficients for the polynomial
);

/** Evaluates a 2-D polynomial at specific coordinates.
 *
 *  @return float    result of polynomial at given location
 */
float psPolynomial2DEval(
    float x,                           ///< x location at which to evaluate
    float y,                           ///< y location at which to evaluate
    const psPolynomial2D* myPoly       ///< Coefficients for the polynomial
);

/** Evaluates a 3-D polynomial at specific coordinates.
 *
 *  @return float    result of polynomial at given location
 */
float psPolynomial3DEval(
    float x,                           ///< x location at which to evaluate
    float y,                           ///< y location at which to evaluate
    float z,                           ///< z location at which to evaluate
    const psPolynomial3D* myPoly       ///< Coefficients for the polynomial
);

/** Evaluates a 4-D polynomial at specific coordinates.
 *
 *  @return float    result of polynomial at given location
 */
float psPolynomial4DEval(
    float w,                           ///< w location at which to evaluate
    float x,                           ///< x location at which to evaluate
    float y,                           ///< y location at which to evaluate
    float z,                           ///< z location at which to evaluate
    const psPolynomial4D* myPoly       ///< Coefficients for the polynomial
);

psVector *psPolynomial1DEvalVector(const psVector *x,             ///< x locations at which to evaluate
                                   const psPolynomial1D *myPoly   ///< Coefficients for the polynomial
                                  );

psVector *psPolynomial2DEvalVector(const psVector *x,             ///< x locations at which to evaluate
                                   const psVector *y,             ///< y locations at which to evaluate
                                   const psPolynomial2D *myPoly   ///< Coefficients for the polynomial
                                  );

psVector *psPolynomial3DEvalVector(const psVector *x,             ///< x locations at which to evaluate
                                   const psVector *y,             ///< y locations at which to evaluate
                                   const psVector *z,             ///< z locations at which to evaluate
                                   const psPolynomial3D *myPoly   ///< Coefficients for the polynomial
                                  );

psVector *psPolynomial4DEvalVector(const psVector *w,             ///< w locations at which to evaluate
                                   const psVector *x,             ///< x locations at which to evaluate
                                   const psVector *y,             ///< y locations at which to evaluate
                                   const psVector *z,             ///< z locations at which to evaluate
                                   const psPolynomial4D *myPoly   ///< Coefficients for the polynomial
                                  );

/*****************************************************************************/

/* Double-precision polynomials, mainly for use in astrometry */

/** Double-precision one-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    int n;                             ///< Number of terms
    double *coeff;                     ///< Coefficients
    double *coeffErr;                  ///< Error in coefficients
    char *mask;                        ///< Coefficient mask
}
psDPolynomial1D;

/** Double-precision two-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    int nX;                            ///< Number of terms in x
    int nY;                            ///< Number of terms in y
    double **coeff;                    ///< Coefficients
    double **coeffErr;                 ///< Error in coefficients
    char **mask;                       ///< Coefficients mask
}
psDPolynomial2D;

/** Double-precision three-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    int nX;                            ///< Number of terms in x
    int nY;                            ///< Number of terms in y
    int nZ;                            ///< Number of terms in z
    double ***coeff;                   ///< Coefficients
    double ***coeffErr;                ///< Error in coefficients
    char ***mask;                      ///< Coefficient mask
}
psDPolynomial3D;

/** Double-precision four-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    int nW;                            ///< Number of terms in w
    int nX;                            ///< Number of terms in x
    int nY;                            ///< Number of terms in y
    int nZ;                            ///< Number of terms in z
    double ****coeff;                  ///< Coefficients
    double ****coeffErr;               ///< Error in coefficients
    char ****mask;                     ///< Coefficients mask
}
psDPolynomial4D;

/** Allocates a double-precision 1-D polynomial structure with n terms
 *
 *  @return  psPolynomial1D*    new double-precision 1-D polynomial struct
 */
psDPolynomial1D* psDPolynomial1DAlloc(
    int n,                             ///< Number of terms
    psPolynomialType type              ///< Polynomial Type
);

/** Allocates a double-precision 2-D polynomial structure
 *
 *  @return  psPolynomial2D*    new double-precision 2-D polynomial struct
 */
psDPolynomial2D* psDPolynomial2DAlloc(
    int nX,                            ///< Number of terms in x
    int nY,                            ///< Number of terms in y
    psPolynomialType type              ///< Polynomial Type
);

/** Allocates a double-precision 3-D polynomial structure
 *
 *  @return  psPolynomial3D*    new double-precision 3-D polynomial struct
 */
psDPolynomial3D* psDPolynomial3DAlloc(
    int nX,                            ///< Number of terms in x
    int nY,                            ///< Number of terms in y
    int nZ,                            ///< Number of terms in z
    psPolynomialType type              ///< Polynomial Type
);

/** Allocates a double-precision 4-D polynomial structure
 *
 *  @return  psPolynomial4D*    new double-precision 4-D polynomial struct
 */
psDPolynomial4D* psDPolynomial4DAlloc(
    int nW,                            ///< Number of terms in w
    int nX,                            ///< Number of terms in x
    int nY,                            ///< Number of terms in y
    int nZ,                            ///< Number of terms in z
    psPolynomialType type              ///< Polynomial Type
);

/** Evaluates a double-precision 1-D polynomial at specific coordinates.
 *
 *  @return float    result of polynomial at given location
 */
double psDPolynomial1DEval(
    double x,                          ///< Value at which to evaluate
    const psDPolynomial1D* myPoly      ///< Coefficients for the polynomial
);

/** Evaluates a double-precision 2-D polynomial at specific coordinates.
 *
 *  @return float    result of polynomial at given location
 */
double psDPolynomial2DEval(
    double x,                           ///< Value x at which to evaluate
    double y,                           ///< Value y at which to evaluate
    const psDPolynomial2D* myPoly       ///< Coefficients for the polynomial
);

/** Evaluates a double-precision 3-D polynomial at specific coordinates.
 *
 *  @return float    result of polynomial at given location
 */
double psDPolynomial3DEval(
    double x,                          ///< Value x at which to evaluate
    double y,                          ///< Value y at which to evaluate
    double z,                          ///< Value z at which to evaluate
    const psDPolynomial3D* myPoly      ///< Coefficients for the polynomial
);

/** Evaluates a double-precision 4-D polynomial at specific coordinates.
 *
 *  @return float    result of polynomial at given location
 */
double psDPolynomial4DEval(
    double w,                          ///< Value w at which to evaluate
    double x,                          ///< Value x at which to evaluate
    double y,                          ///< Value y at which to evaluate
    double z,                          ///< Value z at which to evaluate
    const psDPolynomial4D* myPoly      ///< Coefficients for the polynomial
);

psVector *psDPolynomial1DEvalVector(const psVector *x,             ///< x locations at which to evaluate
                                    const psDPolynomial1D *myPoly  ///< Coefficients for the polynomial
                                   );

psVector *psDPolynomial2DEvalVector(const psVector *x,             ///< x locations at which to evaluate
                                    const psVector *y,             ///< y locations at which to evaluate
                                    const psDPolynomial2D *myPoly  ///< Coefficients for the polynomial
                                   );

psVector *psDPolynomial3DEvalVector(const psVector *x,             ///< x locations at which to evaluate
                                    const psVector *y,             ///< y locations at which to evaluate
                                    const psVector *z,             ///< z locations at which to evaluate
                                    const psDPolynomial3D *myPoly  ///< Coefficients for the polynomial
                                   );

psVector *psDPolynomial4DEvalVector(const psVector *w,             ///< w locations at which to evaluate
                                    const psVector *x,             ///< x locations at which to evaluate
                                    const psVector *y,             ///< y locations at which to evaluate
                                    const psVector *z,             ///< z locations at which to evaluate
                                    const psDPolynomial4D *myPoly  ///< Coefficients for the polynomial
                                   );



typedef struct
{
    int n;                        ///< The number of spline polynomials
    psPolynomial1D **spline;      ///< An array of n pointers to the spline polynomials
    float *p_psDeriv2;            ///< For cubic splines, the second derivative at each domain point.  Size is n+1.
    float *domains;               ///< The boundaries between each spline piece.  Size is n+1.
}
psSpline1D;

psSpline1D *psSpline1DAlloc(int n,
                            int order,
                            float min,
                            float max);

int p_psSpline1DFree(psSpline1D *tmpSpline);

psSpline1D *psSpline1DAllocGeneric(const psVector *bounds,
                                   int order);

float psSpline1DEval(const psSpline1D *spline,
                     float x);

psVector *psSpline1DEvalVector(const psVector *x,
                               const psSpline1D *spline);

int p_psVectorBinDisectF32(float *bins,
                           int numBins,
                           float x);

int p_psVectorBinDisectS32(int *bins,
                           int numBins,
                           int x);

int p_psVectorBinDisect(psVector *bins,
                        psScalar *x);

psScalar *p_psVectorInterpolate(psVector *domain,
                                psVector *range,
                                int order,
                                psScalar *x);

float p_psNRSpline1DEval(psSpline1D *spline,
                         const psVector* restrict x,
                         const psVector* restrict y,
                         float X);

/* \} */// End of MathGroup Functions

#endif

