// This is a part of the Object Oriented Mathematical Library
// Copyright (C) 1999-2008 Renishaw plc
//
// This source file is the confidential property of Renishaw plc.
// Reproduction or transmission in whole or in part, in any form or
// by any means, electronic, mechanical or otherwise, is prohibited
// without the prior written consent of the copyright owner.
//
// objnumerical.h -- numerical methods


// ========================================================================================================================
#ifndef __OBJNUMERICAL_H__
#define __OBJNUMERICAL_H__

// ========================================================================================================================
#include "OOMathsExport.h"

#include "objarray.h"
#include "objalgeb.h"

////////////////////////////////////////////////////////////////////////////////
// CMaNumericalAnalysis - math library numerical methods (abstract base class)

#define INFINITY MaGetLarge()     // 1.0e+300

// ========================================================================================================================
class OOMATHS_API CMaNumericalAnalysis : public CMaObject
{
public:
// destructor
   virtual ~CMaNumericalAnalysis();

protected:
// constructor
   CMaNumericalAnalysis();
};

// ========================================================================================================================
////////////////////////////////////////////////////////////////////////////////
// CMaIntegral - math library integral over an interval [a, b] (concrete class)

class OOMATHS_API CMaIntegral : public CMaNumericalAnalysis
{
public:

// integral types
   enum EMaIntegralType
   {
      // proper integrals
      Proper,           // finite limits, nonsingular smooth (e.g. analytic) integrand

      // improper integrals
      SingularLower,    // finite limits, integrable singularity at lower limit
      SingularUpper,    // finite limits, integrable singularity at upper limit
      InfiniteLower,    // negative infinite lower limit and negative finite upper limit
      InfiniteUpper,    // positive finite lower limit and positive infinite upper limit
      InfiniteExp,      // infinite upper limit with a negative exponential integrand

      // Gauss quadratures based on following weights
      Chebyshev,        // Chebyshev weight function
      Jacobi,           // Jacobi weight function
      Legendre,         // Legendre weight function
      Laguerre,         // Laguerre weight function
      Hermite           // Hermite weight function
   };
   // note: on SingularLower and SingularUpper (integrable singularity)
   //    1. definition: integrand goes to infinite value for a particular x in interval [a, b];
   //          or integrand cannot be computed at this x, but it goes to a finite limiting value.
   //    2. implementation: singularity is allowed at only one of limits. One integral is to be
   //          split into several integrals to satisfy this requirement.
   //
   // note: on InfiniteLower and InfiniteUpper
   //    3. you need to break the integral into two pieces at some value to make the lower limit
   //       and upper limit have the same sign.

// operations
   double Calculate(double (*f)(double x),            // integrand function f(x)
                    double a = 0.0,                   // lower limit value
                    double b = 1.0,                   // upper limit value
                    EMaIntegralType type = Proper,    // integral type
                    double tolerance = MaGetTol(),         // expected precision
                    double alpha = 0.0,               // additional parameter for Gauss quadratures
                    double beta  = 0.0);              // additional parameter for Gauss quadratures

private:
// implementation
   EMaIntegralType m_nType;
   double m_A;
   double m_B;
};

// ========================================================================================================================
////////////////////////////////////////////////////////////////////////////////
// CMaIntegral2D - math library integral over a two-dimensional space (concrete class)

class OOMATHS_API CMaIntegral2D : public CMaNumericalAnalysis
{
public:
// operations
   double Calculate(double (*f)(double x, double y),  // integrand function f(x,y)
                    double (*y1)(double x),           // lower limit function y1(x)
                    double (*y2)(double x),           // upper limit function y2(x)
                    double x1,                        // lower limit value x1
                    double x2);                       // upper limit value x2
};

// ========================================================================================================================
////////////////////////////////////////////////////////////////////////////////
// CMaIntegral3D - math library integral over a three-dimensional space (concrete class)

class OOMATHS_API CMaIntegral3D : public CMaNumericalAnalysis
{
public:
// operations
   double Calculate(double (*f)(double x, double y, double z), // integrand function f(x,y,z)
                    double (*z1)(double x, double y), // lower limit function z1(x,y)
                    double (*z2)(double x, double y), // upper limit function z2(x,y)
                    double (*y1)(double x),           // lower limit function y1(x) 
                    double (*y2)(double x),           // upper limit function y2(x) 
                    double x1,                        // lower limit value x1
                    double x2);                       // upper limit value x2
};

// ========================================================================================================================
////////////////////////////////////////////////////////////////////////////////
// CMaInterpolation - math library interpolation or extrapolation (concrete class)

// formula: yi = f(xi), i = 0,1,..,n, with x0 < x1 < ... < xn
// find: y = f(x)

class OOMATHS_API CMaInterpolation : public CMaNumericalAnalysis
{
public:
// interpolation types
   enum EMaInterpolationType
   {
      Polynomial,
      Rational,     // quotient of polynomials
      CubicSpline
   };

// constructor
   CMaInterpolation(EMaInterpolationType type = Polynomial);

// destructor
   ~CMaInterpolation();

// operations
   // set data
   void SetDegree(int degree = 3);       // not for CubicSpline
   void SetData(const CMaDoubleArray& X, const CMaDoubleArray& Y);

   void SetFirstDerivatives(double y1 = INFINITY, double y2 = INFINITY);
   // note: on SetFirstDerivatives for CubicSpline
   //    specify first derivatives on function's boundaries;
   //    if one or two are infinite, then a natural cubic spline is indicated with
   //    second derivatives being zeros on the boundaries.

   // calculate
   double ValueAt(double x);
   double ErrorAt(double x);             // not for CubicSpline

   // coefficient of interpolating polynomial
   CMaPolynomial ToPolynomial() const;   // for polynomial interpolation only

private:
// implementation
   int m_nD;
   double m_x, m_error, m_y1, m_y2;
   EMaInterpolationType m_nType;
   CMaDoubleArray m_aX, m_aY;      // tabulated data
   CMaDoubleArray m_aY2;           // second derivatives of interpolating function
};

// ========================================================================================================================
////////////////////////////////////////////////////////////////////////////////
// CMaInterpolation2D - math library interpolation in two dimensions (concrete class)

// formula: zij = f(xi, yj), i = 0,1,...,m and j = 0,1,..,n,
//          with x0 < x1 < ... < xm and y0 < y1 < ... < yn
// find: z = f(x, y), with x0 <= x <= xm and y0 <= y <= yn

class OOMATHS_API CMaInterpolation2D : public CMaNumericalAnalysis
{
public:
// interpolation types
   enum EMaInterpolationType
   {
      BiPolynomial,
      BiCubicSpline   // natural spline
   };

// constructor
   CMaInterpolation2D(EMaInterpolationType type = BiPolynomial);

// destructor
   ~CMaInterpolation2D();

// operations
   // set data
   void SetDegrees(int xDegree = 3, int yDegree = 3);   // for BiPolynomial
   void SetData(const CMaMatrix& Z, const CMaDoubleArray& X, const CMaDoubleArray& Y);

   // calculate
   double ValueAt(double x, double y);
   double ErrorAt(double x, double y);                  // for BiPolynomial

private:
// implementation
   int m_nDx, m_nDy;
   double m_x, m_y, m_error;
   EMaInterpolationType m_nType;
   CMaDoubleArray m_aX, m_aY;
   CMaMatrix  m_aZ;               // tabulated data
   CMaMatrix* m_aZ2;              // second derivatives of interpolating function
};

// ========================================================================================================================

#endif // __OBJNUMERICAL_H__
// ========================================================================================================================
