#if !defined(PS_FUNCTIONS_H)
#define PS_FUNCTIONS_H

/** \file psFunctions.h
 *  \brief Standard Mathematical Functions. 
 *  \ingroup MathGroup
 */

/** 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[ 1/(\sqrt(2\pi)\sigma) exp(-\frac{(x-mean)^2}{2\sigma^2}) \f]
 */
float
psGaussian(float x,			///< Value at which to evaluate
	   float mean,			///< Mean for the Gaussian
	   float sigma,			///< Standard deviation for the Gaussian
	   bool normal)			///< normalized Gaussian?
;

psVector *
psGaussianDev(float mean,		///< Mean for the Gaussian
	      float sigma,		///< Standard deviation for the Gaussian
	      int Npts)			///< Number of points generated
;

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

/** Type of polynomial */
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 *restrict coeff;		///< Coefficients
    float *restrict coeffErr;		///< Error in coefficients
    bool *restrict mask;		///< Coefficient mask
} psPolynomial1D;

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

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

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


/** Functions **************************************************************/
/** \addtogroup MathGroup Math Utilities
 *  \{
 */

/** Constructor */
psPolynomial1D *psPolynomial1DAlloc(int n ///< Number of terms
				    psPolynomialType type ///< Polynomial type
				    );
/** Constructor */
psPolynomial2D *psPolynomial2DAlloc(int nX, int nY ///< Number of terms in x and y
				    psPolynomialType type ///< Polynomial type
				    );
/** Constructor */
psPolynomial3D *psPolynomial3DAlloc(int nX, int nY, int nZ ///< Number of terms in x, y and z
				    psPolynomialType type ///< Polynomial type
				    );
/** Constructor */
psPolynomial4D *psPolynomial4DAlloc(int nW, int nX, int nY, int nZ ///< Number of terms in w, x, y and z
				    psPolynomialType type ///< Polynomial type
				    );

/** Destructor */
void psPolynomial1DFree(psPolynomial1D *restrict myPoly ///< Polynomial to destroy
    );
     
/** Destructor */
void psPolynomial2DFree(psPolynomial2D *restrict myPoly ///< Polynomial to destroy
    );
/** Destructor */
void psPolynomial3DFree(psPolynomial3D *restrict myPoly ///< Polynomial to destroy
    );
/** Destructor */
void psPolynomial4DFree(psPolynomial4D *restrict myPoly ///< Polynomial to destroy
    );

/** Evaluate 1D polynomial */
float
psPolynomial1DEval(float x,		///< Value at which to evaluate
		   const psPolynomial1D *myPoly ///< Coefficients for the polynomial
		   );

/** Evaluate 2D polynomial */
float
psPolynomial2DEval(float x,		///< Value x at which to evaluate
		   float y,		///< Value y at which to evaluate
		   const psPolynomial2D *myPoly ///< Coefficients for the polynomial
		   );

/** Evaluate 3D polynomial */
float
psPolynomial3DEval(float x,		///< Value x at which to evaluate
		   float y,		///< Value y at which to evaluate
		   float z,		///< Value z at which to evaluate
		   const psPolynomial3D *myPoly ///< Coefficients for the polynomial
		   );

/** Evaluate 4D polynomial */
float
psPolynomial4DEval(float w,		///< Value w at which to evaluate
		   float x,		///< Value x at which to evaluate
		   float y,		///< Value y at which to evaluate
		   float z,		///< Value z at which to evaluate
		   const psPolynomial4D *myPoly ///< Coefficients for the polynomial
		   );

/** Evaluate 1D polynomial: vector input */
psVector *
psPolynomial1DEvalVector(const psVector *x, ///< Value at which to evaluate
			 const psPolynomial1D *myPoly ///< Coefficients for the polynomial
			 );

/** Evaluate 2D polynomial: vector input */
psVector *
psPolynomial2DEvalVector(const psVector *x, ///< Value x at which to evaluate
			 const psVector *y, ///< Value y at which to evaluate
			 const psPolynomial2D *myPoly ///< Coefficients for the polynomial
			 );

/** Evaluate 3D polynomial: vector input */
psVector *
psPolynomial3DEvalVector(const psVector *x, ///< Value x at which to evaluate
			 const psVector *y, ///< Value y at which to evaluate
			 const psVector *z, ///< Value z at which to evaluate
			 const psPolynomial3D *myPoly ///< Coefficients for the polynomial
			 );

/** Evaluate 4D polynomial: vector input */
psVector *
psPolynomial4DEvalVector(const psVector *w, ///< Value w at which to evaluate
			 const psVector *x, ///< Value x at which to evaluate
			 const psVector *y, ///< Value y at which to evaluate
			 const psVector *z, ///< Value z at which to evaluate
			 const psPolynomial4D *myPoly ///< Coefficients for the polynomial
			 );


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

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

/* 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 *restrict coeff;		///< Coefficients
    double *restrict coeffErr;		///< Error in coefficients
    bool *restrict mask;		///< Coefficient mask
} psDPolynomial1D;

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

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

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

/** Functions **************************************************************/
/** \addtogroup MathGroup Math Utilities
 *  \{
 */

/** Constructor */
psDPolynomial1D *psDPolynomial1DAlloc(int n, ///< Number of terms
				      psPolynomialType type ///< Polynomial type
				      );
/** Constructor */
psDPolynomial2D *psDPolynomial2DAlloc(int nX, int nY ///< Number of terms in x and y
				      psPolynomialType type ///< Polynomial type
				      );
/** Constructor */
psDPolynomial3D *psDPolynomial3DAlloc(int nX, int nY, int nZ ///< Number of terms in x, y and z
				      psPolynomialType type ///< Polynomial type
				      );
/** Constructor */
psDPolynomial4D *psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ ///< Number of terms in w, x, y and z
				      psPolynomialType type ///< Polynomial type
				      );

/** Evaluate 1D polynomial (double precision) */
double
psDPolynomial1DEval(double x,		///< Value at which to evaluate
		    const psDPolynomial1D *myPoly ///< Coefficients for the polynomial
		    );

/** Evaluate 2D polynomial (double precision) */
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
		    );

/** Evaluate 3D polynomial (double precision) */
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
		    );

/** Evaluate 4D polynomial (double precision) */
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
		    );

/** Evaluate 1D polynomial: vector input */
psVector *
psDPolynomial1DEvalVector(const psVector *x, ///< Value at which to evaluate
			  const psDPolynomial1D *myPoly ///< Coefficients for the polynomial
			  );

/** Evaluate 2D polynomial: vector input */
psVector *
psDPolynomial2DEvalVector(const psVector *x, ///< Value x at which to evaluate
			  const psVector *y, ///< Value y at which to evaluate
			  const psDPolynomial2D *myPoly ///< Coefficients for the polynomial
			  );

/** Evaluate 3D polynomial: vector input */
psVector *
psDPolynomial3DEvalVector(const psVector *x, ///< Value x at which to evaluate
			  const psVector *y, ///< Value y at which to evaluate
			  const psVector *z, ///< Value z at which to evaluate
			  const psDPolynomial3D *myPoly ///< Coefficients for the polynomial
			  );

/** Evaluate 4D polynomial: vector input */
psVector *
psDPolynomial4DEvalVector(const psVector *w, ///< Value w at which to evaluate
			  const psVector *x, ///< Value x at which to evaluate
			  const psVector *y, ///< Value y at which to evaluate
			  const psVector *z, ///< Value z at which to evaluate
			  const psDPolynomial4D *myPoly ///< Coefficients for the polynomial
			  );

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

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

/* Splines */

/** A 1D spline */
typedef struct {
    int n;				///< Number of spline pieces
    psPolynomial1D **spline;		///< Array of n pointers to splines
    float *domains;			///< The boundaries between each spline piece.  Size is n+1.
} psSpline1D;

/** Functions **************************************************************/
/** \addtogroup MathGroup Math Utilities
 *  \{
 */

/** Constructors */
psSpline1D *psSpline1DAlloc(int n,	///< Number of spline pieces
			    int order,	///< Order of spline pieces
			    float min,	///< Minimum data value
			    float max	///< Maximum data value
			    );
psSpline1D *psSpline1DAllocGeneric(const psVector *bounds ///< Boundaries between each spline piece.  Number
							  ///< of spline pieces can hence be inferred.
				   int order	///< Order of spline pieces
				   );
/** Destructor */
void psSpline1DFree(psSpline1D *mySpline ///< Spline to destroy
		    );
/** Evaluator */
float psSpline1DEval(float x,		///< Coordinate at which to evaluate
		     const psSpline1D *spline ///< Spline to evaluate
		     );

/** Evaluator: vector version */
psVector *psSpline1DEvalVector(const psVector *x, ///< Coordinates at which to evaluate
			       const psSpline1D *spline ///< Spline to evaluate
			       );

#endif
