/** @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 GLG, MHPCC
*
*  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
*  @date $Date: 2004-12-29 22:15:51 $
*
*  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 psF32      value on the gaussian curve given the input parameters
 */
psF32 psGaussian(
    psF32 x,                           ///< Value at which to evaluate
    psF32 mean,                        ///< Mean for the Gaussian
    psF32 stddev,                      ///< Standard deviation for the Gaussian
    psBool 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* p_psGaussianDev(
    psF32 mean,                        ///< The mean of the Gaussian
    psF32 sigma,                       ///< The sigma of the Gaussian
    psS32 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
    psS32 n;                             ///< Number of terms
    psF32 *coeff;                      ///< Coefficients
    psF32 *coeffErr;                   ///< Error in coefficients
    psU8 *mask;                        ///< Coefficient mask
}
psPolynomial1D;

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

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

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


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

/** Allocates a 2-D polynomial structure
 *
 *  @return  psPolynomial2D*    new 2-D polynomial struct
 */
psPolynomial2D* psPolynomial2DAlloc(
    psS32 nX,                            ///< Number of terms in x
    psS32 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(
    psS32 nX,                            ///< Number of terms in x
    psS32 nY,                            ///< Number of terms in y
    psS32 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(
    psS32 nW,                            ///< Number of terms in w
    psS32 nX,                            ///< Number of terms in x
    psS32 nY,                            ///< Number of terms in y
    psS32 nZ,                            ///< Number of terms in z
    psPolynomialType type              ///< Polynomial Type
);

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

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

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

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

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

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

psVector *psPolynomial3DEvalVector(
    const psPolynomial3D *myPoly,   ///< Coefficients for the polynomial
    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
);

psVector *psPolynomial4DEvalVector(
    const psPolynomial4D *myPoly,   ///< Coefficients for the polynomial
    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
);

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

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

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

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

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

/** Double-precision four-dimensional polynomial */
typedef struct
{
    psPolynomialType type;             ///< Polynomial type
    psS32 nW;                            ///< Number of terms in w
    psS32 nX;                            ///< Number of terms in x
    psS32 nY;                            ///< Number of terms in y
    psS32 nZ;                            ///< Number of terms in z
    psF64 ****coeff;                  ///< Coefficients
    psF64 ****coeffErr;               ///< Error in coefficients
    psU8 ****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(
    psS32 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(
    psS32 nX,                            ///< Number of terms in x
    psS32 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(
    psS32 nX,                            ///< Number of terms in x
    psS32 nY,                            ///< Number of terms in y
    psS32 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(
    psS32 nW,                            ///< Number of terms in w
    psS32 nX,                            ///< Number of terms in x
    psS32 nY,                            ///< Number of terms in y
    psS32 nZ,                            ///< Number of terms in z
    psPolynomialType type              ///< Polynomial Type
);

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

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

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

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

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

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

psVector *psDPolynomial3DEvalVector(
    const psDPolynomial3D *myPoly,  ///< Coefficients for the polynomial
    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
);

psVector *psDPolynomial4DEvalVector(
    const psDPolynomial4D *myPoly,  ///< Coefficients for the polynomial
    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
);



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

psSpline1D *psSpline1DAlloc(psS32 n,
                            psS32 order,
                            psF32 min,
                            psF32 max);

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

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

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

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

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

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

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

#endif

