#ifndef __JMATH__JLEGENDRE__ #define __JMATH__JLEGENDRE__ #include #include #include #include #include"JMath/JZero.hh" #include"JMath/JMathSupportkit.hh" /** * \author mdejong */ namespace JMATH {} namespace JPP { using namespace JMATH; } namespace JMATH { /** * Base class for Legendre polynome. */ template struct JLegendre_t : public std::vector { typedef typename std::vector::const_iterator const_iterator; typedef typename std::vector::iterator iterator; typedef typename std::vector::const_reverse_iterator const_reverse_iterator; typedef typename std::vector::reverse_iterator reverse_iterator; /** * Default constructor. */ JLegendre_t() : xmin(-1.0), xmax(+1.0) {} /** * Constructor. * * \param xmin minimal abscissa value * \param xmax maximal abscissa value */ JLegendre_t(const double xmin, const double xmax) : xmin(xmin), xmax(xmax) {} /** * Virtual destructor. */ virtual ~JLegendre_t() {} /** * Function value. * * \param x abscissa value * \return function value */ virtual JOrdinate_t getValue(const double x) const = 0; /** * Function value. * * \param x abscissa value * \return function value */ JOrdinate_t operator()(const double x) const { return this->getValue(x); } /** * Get minimal abscissa value. * * \return minimal abscissa value */ double getXmin() const { return xmin; } /** * Get maximal abscissa value. * * \return maximal abscissa value */ double getXmax() const { return xmax; } /** * Get normalised abscissa value. * * \param x abscissa value * \return normalised abscissa value */ double getX(const double x) const { return 2.0 * (x - xmin) / (xmax - xmin) - 1.0; } /** * Reset. */ void reset() { for (iterator i = this->begin(); i != this->end(); ++i) { *i = zero; } } /** * Set abscissa values. * * \param xmin minimal abscissa value * \param xmax maximal abscissa value */ void set(const double xmin, const double xmax) { this->xmin = xmin; this->xmax = xmax; } /** * Set abscissa values. * * The template argument T refers to an iterator of a data structure which should have the following data member: * - double %first; // * which conforms with std::pair. * * \param __begin begin of data * \param __end end of data */ template void set(T __begin, T __end) { this->xmin = std::numeric_limits::max(); this->xmax = std::numeric_limits::lowest(); for (T i = __begin; i != __end; ++i) { if (i->first < this->xmin) { this->xmin = i->first; } if (i->first > this->xmax) { this->xmax = i->first; } } } /** * Read Legendre polynome from input. * * \param in input stream * \param object Legendre polynome * \return input stream */ friend inline std::istream& operator>>(std::istream& in, JLegendre_t& object) { object.clear(); for (JOrdinate_t x; in >> x; ) { object.push_back(x); } return in; } /** * Write Legendre polynome to output. * * \param out output stream * \param object Legendre polynome * \return output stream */ friend inline std::ostream& operator<<(std::ostream& out, const JLegendre_t& object) { for (typename JLegendre_t::const_iterator i = object.begin(); i != object.end(); ++i) { out << ' ' << *i; } return out; } protected: double xmin; //!< minimal abscissa double xmax; //!< maximal abscissa }; /** * Template definition for function evaluation of Legendre polynome. */ template struct JLegendre; /** * Template specialisation for function evaluation of of Legendre polynome for undefined number of degrees. */ template struct JLegendre : public JLegendre_t { /** * Default constructor. */ JLegendre() {} /** * Constructor. * * \param xmin minimal abscissa value * \param xmax maximal abscissa value */ JLegendre(const double xmin, const double xmax) { this->set(xmin, xmax); } /** * Function value. * * \param x abscissa value * \return function value */ virtual JOrdinate_t getValue(const double x) const override { const double z = this->getX(x); JOrdinate_t y = zero; for (size_t n = 0; n != this->size(); ++n) { y += (*this)[n] * legendre(n,z); } return y; } }; /** * Template specialisation for function evaluation of of Legendre polynome for defined number of degrees. */ template struct JLegendre : JLegendre { using JLegendre_t::set; /** * Default constructor. */ JLegendre() { this->JLegendre::resize(N+1, getZero()); } /** * Constructor. * * \param xmin minimal abscissa value * \param xmax maximal abscissa value */ JLegendre(const double xmin, const double xmax) { this->JLegendre::resize(N+1, getZero()); this->set(xmin, xmax); } /** * Constructor. * * The template argument T refers to an iterator of a data structure which should have the following data members: * - double %first; // * - JOrdinate_t %second; // * which conforms with std::pair. * * \param __begin begin of data * \param __end end of data */ template JLegendre(T __begin, T __end) { this->JLegendre::resize(N+1, getZero()); this->set(__begin, __end); } /** * Set Legendre polynome. * * The template argument T refers to an iterator of a data structure which should have the following data members: * - double %first; // * - JOrdinate_t %second; // * which conforms with std::pair. * * \param __begin begin of data * \param __end end of data */ template void set(T __begin, T __end) { this->reset(); this->JLegendre::set(__begin, __end); for (size_t n = 0; n != N + 1; ++n) { JOrdinate_t V = zero; double W = 0.0; for (T i = __begin; i != __end; ++i) { const double x = i->first; const JOrdinate_t& y = i->second; const double z = this->getX(x); const double w = legendre(n, z); V += w * (y - (*this)(x)); W += w * (w); } (*this)[n] = V / W; } } /** * Read Legendre polynome from input. * * \param in input stream * \param object Legendre polynome * \return input stream */ friend inline std::istream& operator>>(std::istream& in, JLegendre& object) { for (typename JLegendre::iterator i = object.begin(); i != object.end(); ++i) { in >> *i; } return in; } private: void clear(); void resize(); void erase(); void push_back(); void pop_back(); }; } #endif