From 893785fe95656edeeacfade1865e5dad3902cb6d Mon Sep 17 00:00:00 2001 From: jowr Date: Fri, 6 Jun 2014 14:07:57 +0200 Subject: [PATCH] More matrix conversions --- include/MatrixMath.h | 108 ++ include/PolyMath.h | 1165 +++++++++-------- src/PolyMath.cpp | 2627 ++++++++++++++++++++------------------- src/Tests/eigenTest.cpp | 27 + 4 files changed, 2130 insertions(+), 1797 deletions(-) diff --git a/include/MatrixMath.h b/include/MatrixMath.h index dfc1acdc..37d31244 100644 --- a/include/MatrixMath.h +++ b/include/MatrixMath.h @@ -10,6 +10,8 @@ #include #include "float.h" +#include + /// A wrapper around std::vector /** This wrapper makes the standard vector multi-dimensional. * A useful thing even though we might not need it that @@ -26,6 +28,111 @@ template struct VectorNd<0,T> { namespace CoolProp{ +/// Convert vectors and matrices +/** Conversion functions for the different kinds of object-like + * parameters. This might be obsolete at a later stage, but now + * it is still needed. + */ +/// @param coefficients matrix containing the ordered coefficients +template std::vector > convert(const Eigen::Matrix &coefficients, std::vector > &result){ + // Eigen uses columns as major axis, this might be faster than the row iteration. + // However, the 2D vector stores things differently, no idea what is faster... + //size_t nRows = coefficients.rows(), nCols = coefficients.cols(); + //std::vector > result(R, std::vector(C, 0)); + std::vector col = std::vector(C, 0); + // extend vector if necessary + result.resize(R, col); + for (size_t i = 0; i < R; ++i) { + result[i].resize(C, 0); + for (size_t j = 0; j < C; ++j) { + result[i][j] = coefficients(i,j); + } + } + return result; +} + +/// @param coefficients matrix containing the ordered coefficients +template void convert(const std::vector > &coefficients, Eigen::Matrix &result){ + size_t nRows = num_rows(coefficients), nCols = num_cols(coefficients); + //Eigen::MatrixBase result(nRows,nCols); + for (size_t i = 0; i < nCols; ++i) { + for (size_t j = 0; j < nRows; ++j) { + result(j,i) = coefficients[j][i]; + } + } + //return result; +} + + + +///// @param coefficients matrix containing the ordered coefficients +//template Eigen::Matrix convert(const std::vector > &coefficients){ +// size_t nRows = num_rows(coefficients), nCols = num_cols(coefficients); +// Eigen::MatrixBase result(nRows,nCols); +// for (size_t i = 0; i < nCols; ++i) { +// for (size_t j = 0; j < nRows; ++j) { +// result(j,i) = coefficients[j][i]; +// } +// } +// return result; +//} +// +///// @param coefficients matrix containing the ordered coefficients +//template void convert(const std::vector > &coefficients, Eigen::Matrix &result){ +// size_t nRows = num_rows(coefficients), nCols = num_cols(coefficients); +// //Eigen::MatrixBase result(nRows,nCols); +// for (size_t i = 0; i < nCols; ++i) { +// for (size_t j = 0; j < nRows; ++j) { +// result(j,i) = coefficients[j][i]; +// } +// } +// //return result; +//} + +// +//template void convert(const std::vector > &coefficients, Eigen::MatrixBase &result){ +// size_t nRows = num_rows(coefficients), nCols = num_cols(coefficients); +// //Eigen::MatrixBase result; +// //if ((R!=nRows) || (C!=nCols)) +// result.resize(nRows,nCols); +// for (size_t i = 0; i < nCols; ++i) { +// for (size_t j = 0; j < nRows; ++j) { +// result(j,i) = coefficients[j][i]; +// } +// } +// //return result; +//} + +//template +//inline void func1(MatrixBase &out_mat ){ +// // Do something then return a matrix +// out_mat = ... +//} + +//template +//Eigen::Matrix +//Multiply(const Eigen::MatrixBase& p1, +// const Eigen::MatrixBase& p2) +//{ +// return (p1 * p2); +//} +// +// +//template +//Eigen::Matrix +//Multiply(const Eigen::MatrixBase& p1, +// const Eigen::MatrixBase& p2) +//{ +// return (p1 * p2); +//} + + + + + + + + /// Templates for printing numbers, vectors and matrices static const char* stdFmt = "%7.3f"; @@ -153,6 +260,7 @@ template std::string vec_to_string(const std::vector > & //template std::string vec_to_string( std::vector const& a, const char *fmt); //template std::string vec_to_string(std::vector > const& A, const char *fmt); + // Forward definitions template std::size_t num_rows (std::vector > const& in); template std::size_t max_cols (std::vector > const& in); diff --git a/include/PolyMath.h b/include/PolyMath.h index 649b587f..06fe4260 100644 --- a/include/PolyMath.h +++ b/include/PolyMath.h @@ -16,30 +16,124 @@ namespace CoolProp{ +///// The base class for all Polynomials +//class Polynomial1D{ +// +//private: +// Eigen::VectorXd coefficients; +// +//public: +// /// Constructors +// Polynomial1D(); +// Polynomial1D(const Eigen::VectorXd &coefficients){ +// this->setCoefficients(coefficients); +// } +// Polynomial1D(const std::vector &coefficients){ +// this->setCoefficients(coefficients); +// } +// +// /// Destructor. No implementation +// virtual ~Polynomial1D(){}; +// +//public: +// /// Set the coefficient vector. +// /// @param coefficients vector containing the ordered coefficients +// bool setCoefficients(const Eigen::VectorXd &coefficients); +// /// @param coefficients vector containing the ordered coefficients +// bool setCoefficients(const std::vector &coefficients); +// +// /// Basic checks for coefficient vectors. +// /** Starts with only the first coefficient dimension +// * and checks the vector length against parameter n. */ +// /// @param n unsigned integer value that represents the desired degree of the polynomial +// bool checkCoefficients(const unsigned int n); +// /// @param coefficients vector containing the ordered coefficients +// /// @param n unsigned integer value that represents the desired degree of the polynomial +// bool checkCoefficients(const Eigen::VectorXd &coefficients, const unsigned int n); +// /// @param coefficients vector containing the ordered coefficients +// /// @param n unsigned integer value that represents the desired degree of the polynomial +// bool checkCoefficients(const std::vector &coefficients, const unsigned int n); +// +//protected: +// /// Integration functions +// /** Integrating coefficients for polynomials is done by dividing the +// * original coefficients by (i+1) and elevating the order by 1 +// * through adding a zero as first coefficient. +// * Some reslicing needs to be applied to integrate along the x-axis. +// * In the brine/solution equations, reordering of the parameters +// * avoids this expensive operation. However, it is included for the +// * sake of completeness. +// */ +// Eigen::VectorXd integrateCoeffs(const Eigen::VectorXd &coefficients); +// std::vector integrateCoeffs(const std::vector &coefficients); +// +// /// Derivative coefficients calculation +// /** Deriving coefficients for polynomials is done by multiplying the +// * original coefficients with i and lowering the order by 1. +// * +// * It is not really deprecated, but untested and therefore a warning +// * is issued. Please check this method before you use it. +// */ +// Eigen::VectorXd deriveCoeffs(const Eigen::VectorXd &coefficients); +// std::vector deriveCoeffs(const std::vector &coefficients); +// +//public: +// /// The core functions to evaluate the polynomial +// /** It is here we implement the different special +// * functions that allow us to specify certain +// * types of polynomials. +// * The derivative might bee needed during the +// * solution process of the solver. It could also +// * be a protected function... +// */ +// double evaluate(const double &x_in); +// double derivative(const double &x_in); +// double solve(const double &y_in); +//}; +/// The base class for all Polynomials +class Polynomial2D { - - - -/// The base class for Polynomials -class BasePolynomial{ +private: + Eigen::MatrixXd coefficients; public: - // Constructor - BasePolynomial(); - // Destructor. No implementation - virtual ~BasePolynomial(){}; + /// Constructors + Polynomial2D(); + Polynomial2D(const Eigen::MatrixXd &coefficients){ + this->setCoefficients(coefficients); + } + Polynomial2D(const std::vector > &coefficients){ + this->setCoefficients(coefficients); + } + + /// Destructor. No implementation + virtual ~Polynomial2D(){}; public: + /// Set the coefficient matrix. + /// @param coefficients matrix containing the ordered coefficients + bool setCoefficients(const Eigen::MatrixXd &coefficients); + bool setCoefficients(const std::vector > &coefficients); + /// Basic checks for coefficient vectors. /** Starts with only the first coefficient dimension - * and checks the vector length against parameter n. */ - bool checkCoefficients(const Eigen::VectorXd &coefficients, const unsigned int n); + * and checks the matrix size against the parameters rows and columns. */ + /// @param rows unsigned integer value that represents the desired degree of the polynomial in the 1st dimension + /// @param columns unsigned integer value that represents the desired degree of the polynomial in the 2nd dimension + bool checkCoefficients(const unsigned int rows, const unsigned int columns); + /// @param coefficients matrix containing the ordered coefficients + /// @param rows unsigned integer value that represents the desired degree of the polynomial in the 1st dimension + /// @param columns unsigned integer value that represents the desired degree of the polynomial in the 2nd dimension bool checkCoefficients(const Eigen::MatrixXd &coefficients, const unsigned int rows, const unsigned int columns); - bool checkCoefficients(const std::vector &coefficients, const unsigned int n); - bool checkCoefficients(const std::vector< std::vector > &coefficients, const unsigned int rows, const unsigned int columns); + /// @param coefficients matrix containing the ordered coefficients + /// @param rows unsigned integer value that represents the desired degree of the polynomial in the 1st dimension + /// @param columns unsigned integer value that represents the desired degree of the polynomial in the 2nd dimension + bool checkCoefficients(const std::vector > &coefficients, const unsigned int rows, const unsigned int columns); +protected: + /// Integration functions /** Integrating coefficients for polynomials is done by dividing the * original coefficients by (i+1) and elevating the order by 1 * through adding a zero as first coefficient. @@ -48,504 +142,587 @@ public: * avoids this expensive operation. However, it is included for the * sake of completeness. */ - std::vector integrateCoeffs(const std::vector &coefficients); - std::vector< std::vector > integrateCoeffs(const std::vector< std::vector > &coefficients, bool axis); + /// @param coefficients matrix containing the ordered coefficients + /// @param axis unsigned integer value that represents the desired direction of integration + Eigen::MatrixXd integrateCoeffs(const Eigen::MatrixXd &coefficients, unsigned int axis = 1); + /// @param coefficients matrix containing the ordered coefficients + /// @param axis unsigned integer value that represents the desired direction of integration + std::vector > integrateCoeffs(const std::vector > &coefficients, unsigned int axis = 1); + /// Derivative coefficients calculation /** Deriving coefficients for polynomials is done by multiplying the * original coefficients with i and lowering the order by 1. * * It is not really deprecated, but untested and therefore a warning * is issued. Please check this method before you use it. */ - std::vector deriveCoeffs(const std::vector &coefficients); - std::vector< std::vector > deriveCoeffs(const std::vector< std::vector > &coefficients, unsigned int axis); - -private: - /** The core of the polynomial wrappers are the different - * implementations that follow below. In case there are - * new calculation schemes available, please do not delete - * the implementations, but mark them as deprecated. - * The old functions are good for debugging since the - * structure is easier to read than the backward Horner-scheme - * or the recursive Horner-scheme. - */ - - /// Simple polynomial function generator. <- Deprecated due to poor performance, use Horner-scheme instead - /** Base function to produce n-th order polynomials - * based on the length of the coefficient vector. - * Starts with only the first coefficient at x^0. */ - DEPRECATED(double simplePolynomial(const std::vector &coefficients, double x)); - DEPRECATED(double simplePolynomial(const std::vector > &coefficients, double x, double y)); - - /// Simple integrated polynomial function generator. - /** Base function to produce integrals of n-th order polynomials based on - * the length of the coefficient vector. - * Starts with only the first coefficient at x^0 */ - ///Indefinite integral in x-direction - double simplePolynomialInt(const std::vector &coefficients, double x); - ///Indefinite integral in y-direction only - double simplePolynomialInt(const std::vector > &coefficients, double x, double y); - - /// Simple integrated polynomial function generator divided by independent variable. - /** Base function to produce integrals of n-th order - * polynomials based on the length of the coefficient - * vector. Starts with only the first coefficient at x^0 */ - ///Indefinite integral of a polynomial divided by its independent variable - double simpleFracInt(const std::vector &coefficients, double x); - ///Indefinite integral of a polynomial divided by its 2nd independent variable - double simpleFracInt(const std::vector > &coefficients, double x, double y); - - /** Simple integrated centred(!) polynomial function generator divided by independent variable. - * We need to rewrite some of the functions in order to - * use central fit. Having a central temperature xbase - * allows for a better fit, but requires a different - * formulation of the fracInt function group. Other - * functions are not affected. - * Starts with only the first coefficient at x^0 */ - ///Helper function to calculate the D vector: - double factorial(double nValue); - double binom(double nValue, double nValue2); - std::vector fracIntCentralDvector(int m, double x, double xbase); - ///Indefinite integral of a centred polynomial divided by its independent variable - double fracIntCentral(const std::vector &coefficients, double x, double xbase); - - /// Horner function generator implementations - /** Represent polynomials according to Horner's scheme. - * This avoids unnecessary multiplication and thus - * speeds up calculation. - */ - double baseHorner(const std::vector &coefficients, double x); - double baseHorner(const std::vector< std::vector > &coefficients, double x, double y); - ///Indefinite integral in x-direction - double baseHornerInt(const std::vector &coefficients, double x); - ///Indefinite integral in y-direction only - double baseHornerInt(const std::vector > &coefficients, double x, double y); - ///Indefinite integral of a polynomial divided by its independent variable - double baseHornerFracInt(const std::vector &coefficients, double x); - ///Indefinite integral of a polynomial divided by its 2nd independent variable - double baseHornerFracInt(const std::vector > &coefficients, double x, double y); - - /** Alternatives - * Simple functions that heavily rely on other parts of this file. - * We still need to check which combinations yield the best - * performance. - */ - ///Derivative in x-direction - double deriveIn2Steps(const std::vector &coefficients, double x); // TODO: Check results! - ///Derivative in terms of x(axis=true) or y(axis=false). - double deriveIn2Steps(const std::vector< std::vector > &coefficients, double x, double y, bool axis); // TODO: Check results! - ///Indefinite integral in x-direction - double integrateIn2Steps(const std::vector &coefficients, double x); - ///Indefinite integral in terms of x(axis=true) or y(axis=false). - double integrateIn2Steps(const std::vector< std::vector > &coefficients, double x, double y, bool axis); - ///Indefinite integral in x-direction of a polynomial divided by its independent variable - double fracIntIn2Steps(const std::vector &coefficients, double x); - ///Indefinite integral in y-direction of a polynomial divided by its 2nd independent variable - double fracIntIn2Steps(const std::vector > &coefficients, double x, double y); - ///Indefinite integral of a centred polynomial divided by its 2nd independent variable - double fracIntCentral2Steps(const std::vector > &coefficients, double x, double y, double ybase); + /// @param coefficients matrix containing the ordered coefficients + /// @param axis unsigned integer value that represents the desired direction of derivation + Eigen::MatrixXd deriveCoeffs(const Eigen::MatrixXd &coefficients, unsigned int axis = 1); + /// @param coefficients matrix containing the ordered coefficients + /// @param axis unsigned integer value that represents the desired direction of derivation + std::vector > deriveCoeffs(const std::vector > &coefficients, unsigned int axis = 1); public: - /** Here we define the functions that should be used by the - * respective implementations. Please do no use any other - * method since this would break the purpose of this interface. - * Note that the functions below are supposed to be aliases - * to implementations declared elsewhere in this file. + /// The core functions to evaluate the polynomial + /** It is here we implement the different special + * functions that allow us to specify certain + * types of polynomials. + * The derivative might bee needed during the + * solution process of the solver. It could also + * be a protected function... */ - - /** Everything related to the normal polynomials goes in this - * section, holds all the functions for evaluating polynomials. - */ - /// Evaluates a one-dimensional polynomial for the given coefficients - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input - virtual inline double polyval(const std::vector &coefficients, double x){ - return baseHorner(coefficients,x); - } - - /// Evaluates a two-dimensional polynomial for the given coefficients - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param y double value that represents the current input in the 2nd dimension - virtual inline double polyval(const std::vector< std::vector > &coefficients, double x, double y){ - return baseHorner(coefficients,x,y); - } - - - /** Everything related to the integrated polynomials goes in this - * section, holds all the functions for evaluating polynomials. - */ - /// Evaluates the indefinite integral of a one-dimensional polynomial - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input - virtual inline double polyint(const std::vector &coefficients, double x){ - return baseHornerInt(coefficients,x); - } - - /// Evaluates the indefinite integral of a two-dimensional polynomial along the 2nd axis (y) - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param y double value that represents the current input in the 2nd dimension - virtual inline double polyint(const std::vector< std::vector > &coefficients, double x, double y){ - return baseHornerInt(coefficients,x,y); - } - - - /** Everything related to the derived polynomials goes in this - * section, holds all the functions for evaluating polynomials. - */ - /// Evaluates the derivative of a one-dimensional polynomial - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input - virtual inline double polyder(const std::vector &coefficients, double x){ - return deriveIn2Steps(coefficients,x); - } - - /// Evaluates the derivative of a two-dimensional polynomial along the 2nd axis (y) - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param y double value that represents the current input in the 2nd dimension - virtual inline double polyder(const std::vector< std::vector > &coefficients, double x, double y){ - return deriveIn2Steps(coefficients,x,y,false); - } - - - /** Everything related to the polynomials divided by one variable goes in this - * section, holds all the functions for evaluating polynomials. - */ - /// Evaluates the indefinite integral of a one-dimensional polynomial divided by its independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current position - virtual inline double polyfracval(const std::vector &coefficients, double x){ - return baseHorner(coefficients,x)/x; - } - - /// Evaluates the indefinite integral of a two-dimensional polynomial divided by its 2nd independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param y double value that represents the current input in the 2nd dimension - virtual inline double polyfracval(const std::vector< std::vector > &coefficients, double x, double y){ - return baseHorner(coefficients,x,y)/y; - } - - - /** Everything related to the integrated polynomials divided by one variable goes in this - * section, holds all the functions for solving polynomials. - */ - /// Evaluates the indefinite integral of a one-dimensional polynomial divided by its independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current position - virtual inline double polyfracint(const std::vector &coefficients, double x){ - return baseHornerFracInt(coefficients,x); - } - - /// Evaluates the indefinite integral of a two-dimensional polynomial divided by its 2nd independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param y double value that represents the current input in the 2nd dimension - virtual inline double polyfracint(const std::vector< std::vector > &coefficients, double x, double y){ - return baseHornerFracInt(coefficients,x,y); - } - - /// Evaluates the indefinite integral of a centred one-dimensional polynomial divided by its independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current position - /// @param xbase central temperature for fitted function - virtual inline double polyfracintcentral(const std::vector &coefficients, double x, double xbase){ - return fracIntCentral(coefficients,x,xbase); - } - - /// Evaluates the indefinite integral of a centred two-dimensional polynomial divided by its 2nd independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param y double value that represents the current input in the 2nd dimension - /// @param ybase central temperature for fitted function - virtual inline double polyfracintcentral(const std::vector< std::vector > &coefficients, double x, double y, double ybase){ - return fracIntCentral2Steps(coefficients,x,y,ybase); - } - - - /** Everything related to the derived polynomials divided by one variable goes in this - * section, holds all the functions for solving polynomials. - */ - /// Evaluates the derivative of a one-dimensional polynomial divided by its independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current position - virtual inline double polyfracder(const std::vector &coefficients, double x){ - throw CoolProp::NotImplementedError("Derivatives of polynomials divided by their independent variable have not been implemented."); // TODO: Implement polyfracder1D - } - - /// Evaluates the derivative of a two-dimensional polynomial divided by its 2nd independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param y double value that represents the current input in the 2nd dimension - virtual inline double polyfracder(const std::vector< std::vector > &coefficients, double x, double y){ - throw CoolProp::NotImplementedError("Derivatives of polynomials divided by their independent variable have not been implemented."); // TODO: Implement polyfracder2D - } - - /// Evaluates the derivative of a centred one-dimensional polynomial divided by its independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current position - /// @param xbase central temperature for fitted function - virtual inline double polyfracdercentral(const std::vector &coefficients, double x, double xbase){ - throw CoolProp::NotImplementedError("Derivatives of polynomials divided by their independent variable have not been implemented."); // TODO: Implement polyfracdercentral1D - } - - /// Evaluates the derivative of a centred two-dimensional polynomial divided by its 2nd independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param y double value that represents the current input in the 2nd dimension - /// @param ybase central temperature for fitted function - virtual inline double polyfracdercentral(const std::vector< std::vector > &coefficients, double x, double y, double ybase){ - throw CoolProp::NotImplementedError("Derivatives of polynomials divided by their independent variable have not been implemented."); // TODO: Implement polyfracdercentral2D - } + /// @param x_in double value that represents the current input in the 1st dimension + /// @param y_in double value that represents the current input in the 2nd dimension + double evaluate(const double &x_in, const double &y_in); + /// @param x_in double value that represents the current input in the 1st dimension + /// @param y_in double value that represents the current input in the 2nd dimension + /// @param axis unsigned integer value that represents the axis to solve for (1=x, 2=y) + double derivative(const double &x_in, const double &y_in, unsigned int axis = 1); + /// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension) + /// @param z_in double value that represents the current output in the 3rd dimension + /// @param axis unsigned integer value that represents the axis to solve for (1=x, 2=y) + double solve(const double &in, const double &z_in, unsigned int axis = 1); }; -/** Implements the function wrapper interface and can be - * used by the solvers. - * TODO: Make multidimensional - */ -class PolyResidual : public FuncWrapper1D { -protected: - enum dims {i1D, i2D}; - /// Object that evaluates the equation - BasePolynomial poly; - /// Current output value - double output, firstDim; - int dim; - std::vector< std::vector > coefficients; -private: - PolyResidual(); -public: - PolyResidual(const std::vector &coefficients, double y); - PolyResidual(const std::vector< std::vector > &coefficients, double x, double z); - virtual ~PolyResidual(){}; - bool is2D(){return (this->dim==i2D);}; - virtual double call(double x); - virtual double deriv(double x); -}; -class PolyIntResidual : public PolyResidual { -public: - PolyIntResidual(const std::vector &coefficients, double y):PolyResidual(coefficients, y){}; - PolyIntResidual(const std::vector< std::vector > &coefficients, double x, double z):PolyResidual(coefficients, x, z){}; - virtual double call(double x); - virtual double deriv(double x); -}; -class PolyFracIntResidual : public PolyResidual { -public: - PolyFracIntResidual(const std::vector &coefficients, double y):PolyResidual(coefficients, y){}; - PolyFracIntResidual(const std::vector< std::vector > &coefficients, double x, double z):PolyResidual(coefficients, x, z){}; - virtual double call(double x); - virtual double deriv(double x); -}; -class PolyFracIntCentralResidual : public PolyResidual { -protected: - double baseVal; -public: - PolyFracIntCentralResidual(const std::vector &coefficients, double y, double xBase):PolyResidual(coefficients, y){this->baseVal = xBase;}; - PolyFracIntCentralResidual(const std::vector< std::vector > &coefficients, double x, double z, double yBase): PolyResidual(coefficients, x, z){this->baseVal = yBase;}; - virtual double call(double x); - virtual double deriv(double x); -}; -class PolyDerResidual : public PolyResidual { -public: - PolyDerResidual(const std::vector &coefficients, double y):PolyResidual(coefficients, y){}; - PolyDerResidual(const std::vector< std::vector > &coefficients, double x, double z):PolyResidual(coefficients, x, z){}; - virtual double call(double x); - virtual double deriv(double x); -}; - - -/** Implements the same public functions as the - * but solves the polynomial for the given value - * instead of evaluating it. - * TODO: This class does not check for bijective - * polynomials and is therefore a little - * fragile. - */ -class PolynomialSolver : public BasePolynomial{ -private: - enum solvers {iNewton, iBrent}; - int uses; - double guess, min, max; - double macheps, tol; - int maxiter; - -public: - // Constructor - PolynomialSolver(); - // Destructor. No implementation - virtual ~PolynomialSolver(){}; - -public: - /** Here we redefine the functions that solve the polynomials. - * These implementations all use the base class to evaluate - * the polynomial during the solution process. - */ - - /** Everything related to the normal polynomials goes in this - * section, holds all the functions for solving polynomials. - */ - /// Solves a one-dimensional polynomial for the given coefficients - /// @param coefficients vector containing the ordered coefficients - /// @param y double value that represents the current input - virtual double polyval(const std::vector &coefficients, double y); - - /// Solves a two-dimensional polynomial for the given coefficients - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param z double value that represents the current output - virtual double polyval(const std::vector< std::vector > &coefficients, double x, double z); - - - /** Everything related to the integrated polynomials goes in this - * section, holds all the functions for solving polynomials. - */ - /// Solves the indefinite integral of a one-dimensional polynomial - /// @param coefficients vector containing the ordered coefficients - /// @param y double value that represents the current output - virtual double polyint(const std::vector &coefficients, double y); - - /// Solves the indefinite integral of a two-dimensional polynomial along the 2nd axis (y) - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param z double value that represents the current output - virtual double polyint(const std::vector< std::vector > &coefficients, double x, double z); - - - /** Everything related to the derived polynomials goes in this - * section, holds all the functions for solving polynomials. - */ - /// Solves the derivative of a one-dimensional polynomial - /// @param coefficients vector containing the ordered coefficients - /// @param y double value that represents the current output - virtual double polyder(const std::vector &coefficients, double y); - - /// Solves the derivative of a two-dimensional polynomial along the 2nd axis (y) - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param z double value that represents the current output - virtual double polyder(const std::vector< std::vector > &coefficients, double x, double z); - - - /** Everything related to the polynomials divided by one variable goes in this - * section, holds all the functions for solving polynomials. - */ - /// Solves the indefinite integral of a one-dimensional polynomial divided by its independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param y double value that represents the current output - virtual double polyfracval(const std::vector &coefficients, double y); - - /// Solves the indefinite integral of a two-dimensional polynomial divided by its 2nd independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param z double value that represents the current output - virtual double polyfracval(const std::vector< std::vector > &coefficients, double x, double z); - - - /** Everything related to the integrated polynomials divided by one variable goes in this - * section, holds all the functions for solving polynomials. - */ - /// Solves the indefinite integral of a one-dimensional polynomial divided by its independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param y double value that represents the current output - virtual double polyfracint(const std::vector &coefficients, double y); - - /// Solves the indefinite integral of a two-dimensional polynomial divided by its 2nd independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param z double value that represents the current output - virtual double polyfracint(const std::vector< std::vector > &coefficients, double x, double z); - - /// Solves the indefinite integral of a centred one-dimensional polynomial divided by its independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param y double value that represents the current output - /// @param xbase central x-value for fitted function - virtual double polyfracintcentral(const std::vector &coefficients, double y, double xbase); - - /// Solves the indefinite integral of a centred two-dimensional polynomial divided by its 2nd independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param z double value that represents the current output - /// @param ybase central y-value for fitted function - virtual double polyfracintcentral(const std::vector< std::vector > &coefficients, double x, double z, double ybase); - - - /** Everything related to the derived polynomials divided by one variable goes in this - * section, holds all the functions for solving polynomials. - */ - /// Solves the derivative of a one-dimensional polynomial divided by its independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param y double value that represents the current output - virtual double polyfracder(const std::vector &coefficients, double y); - - /// Solves the derivative of a two-dimensional polynomial divided by its 2nd independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param z double value that represents the current output - virtual double polyfracder(const std::vector< std::vector > &coefficients, double x, double z); - - /// Solves the derivative of a centred one-dimensional polynomial divided by its independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param y double value that represents the current output - /// @param xbase central x-value for fitted function - virtual double polyfracdercentral(const std::vector &coefficients, double y, double xbase); - - /// Solves the derivative of a centred two-dimensional polynomial divided by its 2nd independent variable - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param z double value that represents the current output - /// @param ybase central y-value for fitted function - virtual double polyfracdercentral(const std::vector< std::vector > &coefficients, double x, double z, double ybase); - - - /** Set the solvers and updates either the guess values or the - * boundaries for the variable to solve for. - */ - /// Sets the guess value for the Newton solver and enables it. - /// @param guess double value that represents the guess value - virtual void setGuess(double guess); - /// Sets the limits for the Brent solver and enables it. - /// @param min double value that represents the lower boundary - /// @param max double value that represents the upper boundary - virtual void setLimits(double min, double max); - /// Solves the equations based on previously defined parameters. - /// @param min double value that represents the lower boundary - /// @param max double value that represents the upper boundary - virtual double solve(PolyResidual &res); -}; - - -/// The base class for exponential functions -class BaseExponential{ - -protected: - BasePolynomial poly; - bool POLYMATH_DEBUG; - -public: - BaseExponential(); - virtual ~BaseExponential(){}; - -public: - /// Evaluates an exponential function for the given coefficients - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input - /// @param n int value that determines the kind of exponential function - double expval(const std::vector &coefficients, double x, int n); - - /// Evaluates an exponential function for the given coefficients - /// @param coefficients vector containing the ordered coefficients - /// @param x double value that represents the current input in the 1st dimension - /// @param y double value that represents the current input in the 2nd dimension - /// @param n int value that determines the kind of exponential function - double expval(const std::vector< std::vector > &coefficients, double x, double y, int n); -}; +// +// +// +// +// +// +// +// +// +///// The base class for Polynomials +//class BasePolynomial{ +// +//public: +// // Constructor +// BasePolynomial(); +// // Destructor. No implementation +// virtual ~BasePolynomial(){}; +// +//public: +// /// Basic checks for coefficient vectors. +// /** Starts with only the first coefficient dimension +// * and checks the vector length against parameter n. */ +// bool checkCoefficients(const Eigen::VectorXd &coefficients, const unsigned int n); +// bool checkCoefficients(const Eigen::MatrixXd &coefficients, const unsigned int rows, const unsigned int columns); +// bool checkCoefficients(const std::vector &coefficients, const unsigned int n); +// bool checkCoefficients(const std::vector< std::vector > &coefficients, const unsigned int rows, const unsigned int columns); +// +// /** Integrating coefficients for polynomials is done by dividing the +// * original coefficients by (i+1) and elevating the order by 1 +// * through adding a zero as first coefficient. +// * Some reslicing needs to be applied to integrate along the x-axis. +// * In the brine/solution equations, reordering of the parameters +// * avoids this expensive operation. However, it is included for the +// * sake of completeness. +// */ +// std::vector integrateCoeffs(const std::vector &coefficients); +// std::vector< std::vector > integrateCoeffs(const std::vector< std::vector > &coefficients, bool axis); +// +// /** Deriving coefficients for polynomials is done by multiplying the +// * original coefficients with i and lowering the order by 1. +// * +// * It is not really deprecated, but untested and therefore a warning +// * is issued. Please check this method before you use it. +// */ +// std::vector deriveCoeffs(const std::vector &coefficients); +// std::vector< std::vector > deriveCoeffs(const std::vector< std::vector > &coefficients, unsigned int axis); +// +//private: +// /** The core of the polynomial wrappers are the different +// * implementations that follow below. In case there are +// * new calculation schemes available, please do not delete +// * the implementations, but mark them as deprecated. +// * The old functions are good for debugging since the +// * structure is easier to read than the backward Horner-scheme +// * or the recursive Horner-scheme. +// */ +// +// /// Simple polynomial function generator. <- Deprecated due to poor performance, use Horner-scheme instead +// /** Base function to produce n-th order polynomials +// * based on the length of the coefficient vector. +// * Starts with only the first coefficient at x^0. */ +// DEPRECATED(double simplePolynomial(const std::vector &coefficients, double x)); +// DEPRECATED(double simplePolynomial(const std::vector > &coefficients, double x, double y)); +// +// /// Simple integrated polynomial function generator. +// /** Base function to produce integrals of n-th order polynomials based on +// * the length of the coefficient vector. +// * Starts with only the first coefficient at x^0 */ +// ///Indefinite integral in x-direction +// double simplePolynomialInt(const std::vector &coefficients, double x); +// ///Indefinite integral in y-direction only +// double simplePolynomialInt(const std::vector > &coefficients, double x, double y); +// +// /// Simple integrated polynomial function generator divided by independent variable. +// /** Base function to produce integrals of n-th order +// * polynomials based on the length of the coefficient +// * vector. Starts with only the first coefficient at x^0 */ +// ///Indefinite integral of a polynomial divided by its independent variable +// double simpleFracInt(const std::vector &coefficients, double x); +// ///Indefinite integral of a polynomial divided by its 2nd independent variable +// double simpleFracInt(const std::vector > &coefficients, double x, double y); +// +// /** Simple integrated centred(!) polynomial function generator divided by independent variable. +// * We need to rewrite some of the functions in order to +// * use central fit. Having a central temperature xbase +// * allows for a better fit, but requires a different +// * formulation of the fracInt function group. Other +// * functions are not affected. +// * Starts with only the first coefficient at x^0 */ +// ///Helper function to calculate the D vector: +// double factorial(double nValue); +// double binom(double nValue, double nValue2); +// std::vector fracIntCentralDvector(int m, double x, double xbase); +// ///Indefinite integral of a centred polynomial divided by its independent variable +// double fracIntCentral(const std::vector &coefficients, double x, double xbase); +// +// /// Horner function generator implementations +// /** Represent polynomials according to Horner's scheme. +// * This avoids unnecessary multiplication and thus +// * speeds up calculation. +// */ +// double baseHorner(const std::vector &coefficients, double x); +// double baseHorner(const std::vector< std::vector > &coefficients, double x, double y); +// ///Indefinite integral in x-direction +// double baseHornerInt(const std::vector &coefficients, double x); +// ///Indefinite integral in y-direction only +// double baseHornerInt(const std::vector > &coefficients, double x, double y); +// ///Indefinite integral of a polynomial divided by its independent variable +// double baseHornerFracInt(const std::vector &coefficients, double x); +// ///Indefinite integral of a polynomial divided by its 2nd independent variable +// double baseHornerFracInt(const std::vector > &coefficients, double x, double y); +// +// /** Alternatives +// * Simple functions that heavily rely on other parts of this file. +// * We still need to check which combinations yield the best +// * performance. +// */ +// ///Derivative in x-direction +// double deriveIn2Steps(const std::vector &coefficients, double x); // TODO: Check results! +// ///Derivative in terms of x(axis=true) or y(axis=false). +// double deriveIn2Steps(const std::vector< std::vector > &coefficients, double x, double y, bool axis); // TODO: Check results! +// ///Indefinite integral in x-direction +// double integrateIn2Steps(const std::vector &coefficients, double x); +// ///Indefinite integral in terms of x(axis=true) or y(axis=false). +// double integrateIn2Steps(const std::vector< std::vector > &coefficients, double x, double y, bool axis); +// ///Indefinite integral in x-direction of a polynomial divided by its independent variable +// double fracIntIn2Steps(const std::vector &coefficients, double x); +// ///Indefinite integral in y-direction of a polynomial divided by its 2nd independent variable +// double fracIntIn2Steps(const std::vector > &coefficients, double x, double y); +// ///Indefinite integral of a centred polynomial divided by its 2nd independent variable +// double fracIntCentral2Steps(const std::vector > &coefficients, double x, double y, double ybase); +// +//public: +// /** Here we define the functions that should be used by the +// * respective implementations. Please do no use any other +// * method since this would break the purpose of this interface. +// * Note that the functions below are supposed to be aliases +// * to implementations declared elsewhere in this file. +// */ +// +// /** Everything related to the normal polynomials goes in this +// * section, holds all the functions for evaluating polynomials. +// */ +// /// Evaluates a one-dimensional polynomial for the given coefficients +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input +// virtual inline double polyval(const std::vector &coefficients, double x){ +// return baseHorner(coefficients,x); +// } +// +// /// Evaluates a two-dimensional polynomial for the given coefficients +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param y double value that represents the current input in the 2nd dimension +// virtual inline double polyval(const std::vector< std::vector > &coefficients, double x, double y){ +// return baseHorner(coefficients,x,y); +// } +// +// +// /** Everything related to the integrated polynomials goes in this +// * section, holds all the functions for evaluating polynomials. +// */ +// /// Evaluates the indefinite integral of a one-dimensional polynomial +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input +// virtual inline double polyint(const std::vector &coefficients, double x){ +// return baseHornerInt(coefficients,x); +// } +// +// /// Evaluates the indefinite integral of a two-dimensional polynomial along the 2nd axis (y) +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param y double value that represents the current input in the 2nd dimension +// virtual inline double polyint(const std::vector< std::vector > &coefficients, double x, double y){ +// return baseHornerInt(coefficients,x,y); +// } +// +// +// /** Everything related to the derived polynomials goes in this +// * section, holds all the functions for evaluating polynomials. +// */ +// /// Evaluates the derivative of a one-dimensional polynomial +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input +// virtual inline double polyder(const std::vector &coefficients, double x){ +// return deriveIn2Steps(coefficients,x); +// } +// +// /// Evaluates the derivative of a two-dimensional polynomial along the 2nd axis (y) +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param y double value that represents the current input in the 2nd dimension +// virtual inline double polyder(const std::vector< std::vector > &coefficients, double x, double y){ +// return deriveIn2Steps(coefficients,x,y,false); +// } +// +// +// /** Everything related to the polynomials divided by one variable goes in this +// * section, holds all the functions for evaluating polynomials. +// */ +// /// Evaluates the indefinite integral of a one-dimensional polynomial divided by its independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current position +// virtual inline double polyfracval(const std::vector &coefficients, double x){ +// return baseHorner(coefficients,x)/x; +// } +// +// /// Evaluates the indefinite integral of a two-dimensional polynomial divided by its 2nd independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param y double value that represents the current input in the 2nd dimension +// virtual inline double polyfracval(const std::vector< std::vector > &coefficients, double x, double y){ +// return baseHorner(coefficients,x,y)/y; +// } +// +// +// /** Everything related to the integrated polynomials divided by one variable goes in this +// * section, holds all the functions for solving polynomials. +// */ +// /// Evaluates the indefinite integral of a one-dimensional polynomial divided by its independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current position +// virtual inline double polyfracint(const std::vector &coefficients, double x){ +// return baseHornerFracInt(coefficients,x); +// } +// +// /// Evaluates the indefinite integral of a two-dimensional polynomial divided by its 2nd independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param y double value that represents the current input in the 2nd dimension +// virtual inline double polyfracint(const std::vector< std::vector > &coefficients, double x, double y){ +// return baseHornerFracInt(coefficients,x,y); +// } +// +// /// Evaluates the indefinite integral of a centred one-dimensional polynomial divided by its independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current position +// /// @param xbase central temperature for fitted function +// virtual inline double polyfracintcentral(const std::vector &coefficients, double x, double xbase){ +// return fracIntCentral(coefficients,x,xbase); +// } +// +// /// Evaluates the indefinite integral of a centred two-dimensional polynomial divided by its 2nd independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param y double value that represents the current input in the 2nd dimension +// /// @param ybase central temperature for fitted function +// virtual inline double polyfracintcentral(const std::vector< std::vector > &coefficients, double x, double y, double ybase){ +// return fracIntCentral2Steps(coefficients,x,y,ybase); +// } +// +// +// /** Everything related to the derived polynomials divided by one variable goes in this +// * section, holds all the functions for solving polynomials. +// */ +// /// Evaluates the derivative of a one-dimensional polynomial divided by its independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current position +// virtual inline double polyfracder(const std::vector &coefficients, double x){ +// throw CoolProp::NotImplementedError("Derivatives of polynomials divided by their independent variable have not been implemented."); // TODO: Implement polyfracder1D +// } +// +// /// Evaluates the derivative of a two-dimensional polynomial divided by its 2nd independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param y double value that represents the current input in the 2nd dimension +// virtual inline double polyfracder(const std::vector< std::vector > &coefficients, double x, double y){ +// throw CoolProp::NotImplementedError("Derivatives of polynomials divided by their independent variable have not been implemented."); // TODO: Implement polyfracder2D +// } +// +// /// Evaluates the derivative of a centred one-dimensional polynomial divided by its independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current position +// /// @param xbase central temperature for fitted function +// virtual inline double polyfracdercentral(const std::vector &coefficients, double x, double xbase){ +// throw CoolProp::NotImplementedError("Derivatives of polynomials divided by their independent variable have not been implemented."); // TODO: Implement polyfracdercentral1D +// } +// +// /// Evaluates the derivative of a centred two-dimensional polynomial divided by its 2nd independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param y double value that represents the current input in the 2nd dimension +// /// @param ybase central temperature for fitted function +// virtual inline double polyfracdercentral(const std::vector< std::vector > &coefficients, double x, double y, double ybase){ +// throw CoolProp::NotImplementedError("Derivatives of polynomials divided by their independent variable have not been implemented."); // TODO: Implement polyfracdercentral2D +// } +//}; +// +// +// +// +///** Implements the function wrapper interface and can be +// * used by the solvers. +// * TODO: Make multidimensional +// */ +//class PolyResidual : public FuncWrapper1D { +//protected: +// enum dims {i1D, i2D}; +// /// Object that evaluates the equation +// BasePolynomial poly; +// /// Current output value +// double output, firstDim; +// int dim; +// std::vector< std::vector > coefficients; +//private: +// PolyResidual(); +//public: +// PolyResidual(const std::vector &coefficients, double y); +// PolyResidual(const std::vector< std::vector > &coefficients, double x, double z); +// virtual ~PolyResidual(){}; +// bool is2D(){return (this->dim==i2D);}; +// virtual double call(double x); +// virtual double deriv(double x); +//}; +//class PolyIntResidual : public PolyResidual { +//public: +// PolyIntResidual(const std::vector &coefficients, double y):PolyResidual(coefficients, y){}; +// PolyIntResidual(const std::vector< std::vector > &coefficients, double x, double z):PolyResidual(coefficients, x, z){}; +// virtual double call(double x); +// virtual double deriv(double x); +//}; +//class PolyFracIntResidual : public PolyResidual { +//public: +// PolyFracIntResidual(const std::vector &coefficients, double y):PolyResidual(coefficients, y){}; +// PolyFracIntResidual(const std::vector< std::vector > &coefficients, double x, double z):PolyResidual(coefficients, x, z){}; +// virtual double call(double x); +// virtual double deriv(double x); +//}; +//class PolyFracIntCentralResidual : public PolyResidual { +//protected: +// double baseVal; +//public: +// PolyFracIntCentralResidual(const std::vector &coefficients, double y, double xBase):PolyResidual(coefficients, y){this->baseVal = xBase;}; +// PolyFracIntCentralResidual(const std::vector< std::vector > &coefficients, double x, double z, double yBase): PolyResidual(coefficients, x, z){this->baseVal = yBase;}; +// virtual double call(double x); +// virtual double deriv(double x); +//}; +//class PolyDerResidual : public PolyResidual { +//public: +// PolyDerResidual(const std::vector &coefficients, double y):PolyResidual(coefficients, y){}; +// PolyDerResidual(const std::vector< std::vector > &coefficients, double x, double z):PolyResidual(coefficients, x, z){}; +// virtual double call(double x); +// virtual double deriv(double x); +//}; +// +// +// +// +///** Implements the same public functions as the +// * but solves the polynomial for the given value +// * instead of evaluating it. +// * TODO: This class does not check for bijective +// * polynomials and is therefore a little +// * fragile. +// */ +//class PolynomialSolver : public BasePolynomial{ +//private: +// enum solvers {iNewton, iBrent}; +// int uses; +// double guess, min, max; +// double macheps, tol; +// int maxiter; +// +//public: +// // Constructor +// PolynomialSolver(); +// // Destructor. No implementation +// virtual ~PolynomialSolver(){}; +// +//public: +// /** Here we redefine the functions that solve the polynomials. +// * These implementations all use the base class to evaluate +// * the polynomial during the solution process. +// */ +// +// /** Everything related to the normal polynomials goes in this +// * section, holds all the functions for solving polynomials. +// */ +// /// Solves a one-dimensional polynomial for the given coefficients +// /// @param coefficients vector containing the ordered coefficients +// /// @param y double value that represents the current input +// virtual double polyval(const std::vector &coefficients, double y); +// +// /// Solves a two-dimensional polynomial for the given coefficients +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param z double value that represents the current output +// virtual double polyval(const std::vector< std::vector > &coefficients, double x, double z); +// +// +// /** Everything related to the integrated polynomials goes in this +// * section, holds all the functions for solving polynomials. +// */ +// /// Solves the indefinite integral of a one-dimensional polynomial +// /// @param coefficients vector containing the ordered coefficients +// /// @param y double value that represents the current output +// virtual double polyint(const std::vector &coefficients, double y); +// +// /// Solves the indefinite integral of a two-dimensional polynomial along the 2nd axis (y) +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param z double value that represents the current output +// virtual double polyint(const std::vector< std::vector > &coefficients, double x, double z); +// +// +// /** Everything related to the derived polynomials goes in this +// * section, holds all the functions for solving polynomials. +// */ +// /// Solves the derivative of a one-dimensional polynomial +// /// @param coefficients vector containing the ordered coefficients +// /// @param y double value that represents the current output +// virtual double polyder(const std::vector &coefficients, double y); +// +// /// Solves the derivative of a two-dimensional polynomial along the 2nd axis (y) +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param z double value that represents the current output +// virtual double polyder(const std::vector< std::vector > &coefficients, double x, double z); +// +// +// /** Everything related to the polynomials divided by one variable goes in this +// * section, holds all the functions for solving polynomials. +// */ +// /// Solves the indefinite integral of a one-dimensional polynomial divided by its independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param y double value that represents the current output +// virtual double polyfracval(const std::vector &coefficients, double y); +// +// /// Solves the indefinite integral of a two-dimensional polynomial divided by its 2nd independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param z double value that represents the current output +// virtual double polyfracval(const std::vector< std::vector > &coefficients, double x, double z); +// +// +// /** Everything related to the integrated polynomials divided by one variable goes in this +// * section, holds all the functions for solving polynomials. +// */ +// /// Solves the indefinite integral of a one-dimensional polynomial divided by its independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param y double value that represents the current output +// virtual double polyfracint(const std::vector &coefficients, double y); +// +// /// Solves the indefinite integral of a two-dimensional polynomial divided by its 2nd independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param z double value that represents the current output +// virtual double polyfracint(const std::vector< std::vector > &coefficients, double x, double z); +// +// /// Solves the indefinite integral of a centred one-dimensional polynomial divided by its independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param y double value that represents the current output +// /// @param xbase central x-value for fitted function +// virtual double polyfracintcentral(const std::vector &coefficients, double y, double xbase); +// +// /// Solves the indefinite integral of a centred two-dimensional polynomial divided by its 2nd independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param z double value that represents the current output +// /// @param ybase central y-value for fitted function +// virtual double polyfracintcentral(const std::vector< std::vector > &coefficients, double x, double z, double ybase); +// +// +// /** Everything related to the derived polynomials divided by one variable goes in this +// * section, holds all the functions for solving polynomials. +// */ +// /// Solves the derivative of a one-dimensional polynomial divided by its independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param y double value that represents the current output +// virtual double polyfracder(const std::vector &coefficients, double y); +// +// /// Solves the derivative of a two-dimensional polynomial divided by its 2nd independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param z double value that represents the current output +// virtual double polyfracder(const std::vector< std::vector > &coefficients, double x, double z); +// +// /// Solves the derivative of a centred one-dimensional polynomial divided by its independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param y double value that represents the current output +// /// @param xbase central x-value for fitted function +// virtual double polyfracdercentral(const std::vector &coefficients, double y, double xbase); +// +// /// Solves the derivative of a centred two-dimensional polynomial divided by its 2nd independent variable +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param z double value that represents the current output +// /// @param ybase central y-value for fitted function +// virtual double polyfracdercentral(const std::vector< std::vector > &coefficients, double x, double z, double ybase); +// +// +// /** Set the solvers and updates either the guess values or the +// * boundaries for the variable to solve for. +// */ +// /// Sets the guess value for the Newton solver and enables it. +// /// @param guess double value that represents the guess value +// virtual void setGuess(double guess); +// /// Sets the limits for the Brent solver and enables it. +// /// @param min double value that represents the lower boundary +// /// @param max double value that represents the upper boundary +// virtual void setLimits(double min, double max); +// /// Solves the equations based on previously defined parameters. +// /// @param min double value that represents the lower boundary +// /// @param max double value that represents the upper boundary +// virtual double solve(PolyResidual &res); +//}; +// +// +///// The base class for exponential functions +//class BaseExponential{ +// +//protected: +// BasePolynomial poly; +// bool POLYMATH_DEBUG; +// +//public: +// BaseExponential(); +// virtual ~BaseExponential(){}; +// +//public: +// /// Evaluates an exponential function for the given coefficients +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input +// /// @param n int value that determines the kind of exponential function +// double expval(const std::vector &coefficients, double x, int n); +// +// /// Evaluates an exponential function for the given coefficients +// /// @param coefficients vector containing the ordered coefficients +// /// @param x double value that represents the current input in the 1st dimension +// /// @param y double value that represents the current input in the 2nd dimension +// /// @param n int value that determines the kind of exponential function +// double expval(const std::vector< std::vector > &coefficients, double x, double y, int n); +//}; }; /* namespace CoolProp */ diff --git a/src/PolyMath.cpp b/src/PolyMath.cpp index 7b109054..86bc3281 100644 --- a/src/PolyMath.cpp +++ b/src/PolyMath.cpp @@ -13,1235 +13,1041 @@ #include "Solvers.h" +#include + namespace CoolProp{ -BasePolynomial::BasePolynomial(){ - this->POLYMATH_DEBUG = false; +/// Constructors for the base class for all Polynomials +//Polynomial1D::Polynomial1D(); +bool Polynomial2D::setCoefficients(const Eigen::MatrixXd &coefficients){ + this.coefficients = coefficients; + return this.coefficients == coefficients; +} +bool Polynomial2D::setCoefficients(const std::vector< std::vector > &coefficients){ + return this->setCoefficients(convert(coefficients)); +} } -/// Basic checks for coefficient vectors. -/** Starts with only the first coefficient dimension - * and checks the vector length against parameter n. */ -bool BasePolynomial::checkCoefficients(const std::vector &coefficients, unsigned int n){ - if (coefficients.size() == n){ - return true; - } else { - throw ValueError(format("The number of coefficients %d does not match with %d. ",coefficients.size(),n)); - } - return false; -} -bool BasePolynomial::checkCoefficients(std::vector< std::vector > const& coefficients, unsigned int rows, unsigned int columns){ - if (coefficients.size() == rows){ - bool result = true; - for(unsigned int i=0; i BasePolynomial::integrateCoeffs(std::vector const& coefficients){ - std::vector newCoefficients; - unsigned int sizeX = coefficients.size(); - if (sizeX<1) throw ValueError(format("You have to provide coefficients, a vector length of %d is not a valid. ",sizeX)); - // pushing a zero elevates the order by 1 - newCoefficients.push_back(0.0); - for(unsigned int i=0; i > BasePolynomial::integrateCoeffs(std::vector< std::vector > const& coefficients, bool axis){ - std::vector< std::vector > newCoefficients; - unsigned int sizeX = coefficients.size(); - if (sizeX<1) throw ValueError(format("You have to provide coefficients, a vector length of %d is not a valid. ",sizeX)); - - if (axis==true){ - std::vector< std::vector > tmpCoefficients; - tmpCoefficients = transpose(coefficients); - unsigned int sizeY = tmpCoefficients.size(); - for(unsigned int i=0; i BasePolynomial::deriveCoeffs(std::vector const& coefficients){ - std::vector newCoefficients; - unsigned int sizeX = coefficients.size(); - if (sizeX<1) throw ValueError(format("You have to provide coefficients, a vector length of %d is not a valid. ",sizeX)); - // skipping the first element lowers the order - for(unsigned int i=1; i > BasePolynomial::deriveCoeffs(const std::vector< std::vector > &coefficients, unsigned int axis){ - std::vector< std::vector > newCoefficients; - unsigned int sizeX = coefficients.size(); - if (sizeX<1) throw ValueError(format("You have to provide coefficients, a vector length of %d is not a valid. ",sizeX)); - if (axis==0){ - std::vector< std::vector > tmpCoefficients; - tmpCoefficients = transpose(coefficients); - unsigned int sizeY = tmpCoefficients.size(); - for(unsigned int i=0; i const& coefficients, double x){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running simplePolynomial(std::vector, " << x << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - double result = 0.; - for(unsigned int i=0; iPOLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} -double BasePolynomial::simplePolynomial(std::vector > const& coefficients, double x, double y){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running simplePolynomial(std::vector, " << x << ", " << y << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - double result = 0; - for(unsigned int i=0; iPOLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - - -/// Simple integrated polynomial function generator. -/** Base function to produce integrals of n-th order - * polynomials based on the length of the coefficient - * vector. - * Starts with only the first coefficient at x^0 */ -///Indefinite integral in x-direction -double BasePolynomial::simplePolynomialInt(std::vector const& coefficients, double x){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running simplePolynomialInt(std::vector, " << x << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - double result = 0.; - for(unsigned int i=0; iPOLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - -///Indefinite integral in y-direction only -double BasePolynomial::simplePolynomialInt(std::vector > const& coefficients, double x, double y){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running simplePolynomialInt(std::vector, " << x << ", " << y << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - double result = 0.; - for(unsigned int i=0; iPOLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - - -/// Simple integrated polynomial function generator divided by independent variable. -/** Base function to produce integrals of n-th order - * polynomials based on the length of the coefficient - * vector. - * Starts with only the first coefficient at x^0 */ -double BasePolynomial::simpleFracInt(std::vector const& coefficients, double x){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running simpleFracInt(std::vector, " << x << "): "; - } - double result = coefficients[0] * log(x); - if (coefficients.size() > 1) { - for (unsigned int i=1; iPOLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - -double BasePolynomial::simpleFracInt(std::vector< std::vector > const& coefficients, double x, double y){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running simpleFracInt(std::vector, " << x << ", " << y << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - double result = 0; - for (unsigned int i=0; iPOLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - - -/** Simple integrated centred(!) polynomial function generator divided by independent variable. - * We need to rewrite some of the functions in order to - * use central fit. Having a central temperature xbase - * allows for a better fit, but requires a different - * formulation of the fracInt function group. Other - * functions are not affected. - * Starts with only the first coefficient at x^0 */ - -///Helper functions to calculate binomial coefficients: http://rosettacode.org/wiki/Evaluate_binomial_coefficients#C.2B.2B -//double BasePolynomial::factorial(double nValue){ -// double result = nValue; -// double result_next; -// double pc = nValue; -// do { -// result_next = result*(pc-1); -// result = result_next; -// pc--; -// } while(pc>2); -// nValue = result; -// return nValue; +//namespace CoolProp{ +// +//BasePolynomial::BasePolynomial(){ +// this->POLYMATH_DEBUG = false; //} -//double BasePolynomial::factorial(double nValue){ -// if (nValue == 0) return (1); -// else return (nValue * factorial(nValue - 1)); +// +// +///// Basic checks for coefficient vectors. +///** Starts with only the first coefficient dimension +// * and checks the vector length against parameter n. */ +//bool BasePolynomial::checkCoefficients(const std::vector &coefficients, unsigned int n){ +// if (coefficients.size() == n){ +// return true; +// } else { +// throw ValueError(format("The number of coefficients %d does not match with %d. ",coefficients.size(),n)); +// } +// return false; //} -double BasePolynomial::factorial(double nValue){ - double value = 1; - for(int i = 2; i <= nValue; i++){ - value = value * i; - } - return value; -} - -double BasePolynomial::binom(double nValue, double nValue2){ - double result; - if(nValue2 == 1) return nValue; - result = (factorial(nValue)) / (factorial(nValue2)*factorial((nValue - nValue2))); - nValue2 = result; - return nValue2; -} - -///Helper functions to calculate the D vector: -std::vector BasePolynomial::fracIntCentralDvector(int m, double x, double xbase){ - std::vector D; - double tmp; - if (m<1) throw ValueError(format("You have to provide coefficients, a vector length of %d is not a valid. ",m)); - for (int j=0; j const& coefficients, double x, double xbase){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running fracIntCentral(std::vector, " << x << ", " << xbase << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - int m = coefficients.size(); - std::vector D = fracIntCentralDvector(m, x, xbase); - double result = 0; - for(int j=0; jPOLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - - -/// Horner function generator implementations -/** Represent polynomials according to Horner's scheme. - * This avoids unnecessary multiplication and thus - * speeds up calculation. - */ -double BasePolynomial::baseHorner(std::vector const& coefficients, double x){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running baseHorner(std::vector, " << x << "): "; - } - double result = 0; - for(int i=coefficients.size()-1; i>=0; i--) { - result *= x; - result += coefficients[i]; - } - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - -double BasePolynomial::baseHorner(std::vector< std::vector > const& coefficients, double x, double y){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running baseHorner(std::vector, " << x << ", " << y << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - double result = 0; - for(int i=coefficients.size()-1; i>=0; i--) { - result *= x; - result += baseHorner(coefficients[i], y); - } - this->POLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - -///Indefinite integral in x-direction -double BasePolynomial::baseHornerInt(std::vector const& coefficients, double x){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running baseHornerInt(std::vector, " << x << "): "; - } - double result = 0; - for(int i=coefficients.size()-1; i>=0; i--) { - result *= x; - result += coefficients[i]/(i+1.); - } - result = result * x; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - -///Indefinite integral in y-direction only -double BasePolynomial::baseHornerInt(std::vector > const& coefficients, double x, double y){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running baseHornerInt(std::vector, " << x << ", " << y << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - double result = 0; - for(int i=coefficients.size()-1; i>=0; i--) { - result *= x; - result += baseHornerInt(coefficients[i], y); - } - this->POLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - -///Indefinite integral in x-direction of a polynomial divided by its independent variable -double BasePolynomial::baseHornerFracInt(std::vector const& coefficients, double x){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running baseHornerFra(std::vector, " << x << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - double result = 0; - if (coefficients.size() > 1) { - for(int i=coefficients.size()-1; i>=1; i--) { - result *= x; - result += coefficients[i]/(i); - } - result *= x; - } - result += coefficients[0] * log(x); - this->POLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - -///Indefinite integral in y-direction of a polynomial divided by its 2nd independent variable -double BasePolynomial::baseHornerFracInt(std::vector > const& coefficients, double x, double y){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running baseHornerFra(std::vector, " << x << ", " << y << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - - double result = 0; - for(int i=coefficients.size()-1; i>=0; i--) { - result *= x; - result += baseHornerFracInt(coefficients[i], y); - } - - this->POLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - - -/** Alternatives - * Simple functions that heavily rely on other parts of this file. - * We still need to check which combinations yield the best - * performance. - */ -///Derivative in x-direction -double BasePolynomial::deriveIn2Steps(std::vector const& coefficients, double x){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running deriveIn2Steps(std::vector, " << x << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - double result = polyval(deriveCoeffs(coefficients),x); - this->POLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - -///Derivative in terms of x(axis=true) or y(axis=false). -double BasePolynomial::deriveIn2Steps(std::vector< std::vector > const& coefficients, double x, double y, bool axis){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running deriveIn2Steps(std::vector, " << x << ", " << y << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - double result = polyval(deriveCoeffs(coefficients,axis),x,y); - this->POLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - -///Indefinite integral in x-direction -double BasePolynomial::integrateIn2Steps(std::vector const& coefficients, double x){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running integrateIn2Steps(std::vector, " << x << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - double result = polyval(integrateCoeffs(coefficients),x); - this->POLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - -///Indefinite integral in terms of x(axis=true) or y(axis=false). -double BasePolynomial::integrateIn2Steps(std::vector< std::vector > const& coefficients, double x, double y, bool axis){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running integrateIn2Steps(std::vector, " << x << ", " << y << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - double result = polyval(integrateCoeffs(coefficients,axis),x,y); - this->POLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - -///Indefinite integral in x-direction of a polynomial divided by its independent variable -double BasePolynomial::fracIntIn2Steps(std::vector const& coefficients, double x){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running fracIntIn2Steps(std::vector, " << x << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - double result = coefficients[0] * log(x); - if (coefficients.size() > 1) { - std::vector newCoeffs(coefficients.begin() + 1, coefficients.end()); - result += polyint(newCoeffs,x); - } - this->POLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - -///Indefinite integral in y-direction of a polynomial divided by its 2nd independent variable -double BasePolynomial::fracIntIn2Steps(std::vector > const& coefficients, double x, double y){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running fracIntIn2Steps(std::vector, " << x << ", " << y << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - std::vector newCoeffs; - for (unsigned int i=0; iPOLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - -///Indefinite integral in y-direction of a centred polynomial divided by its 2nd independent variable -double BasePolynomial::fracIntCentral2Steps(std::vector > const& coefficients, double x, double y, double ybase){ - if (this->POLYMATH_DEBUG) { - std::cout << "Running fracIntCentral2Steps(std::vector, " << x << ", " << y << ", " << ybase << "): "; - } - bool db = this->POLYMATH_DEBUG; - this->POLYMATH_DEBUG = false; - std::vector newCoeffs; - for (unsigned int i=0; iPOLYMATH_DEBUG = db; - if (this->POLYMATH_DEBUG) { - std::cout << result << std::endl; - } - return result; -} - - - - -/** Implements the function wrapper interface and can be - * used by the solvers. This is only an example and you should - * use local redefinitions of the class. - * TODO: Make multidimensional - */ -PolyResidual::PolyResidual(){ - this->dim = -1; -} - -PolyResidual::PolyResidual(const std::vector &coefficients, double y){ - this->output = y; - this->firstDim = 0; - this->coefficients.clear(); - this->coefficients.push_back(coefficients); - this->dim = i1D; -} - -PolyResidual::PolyResidual(const std::vector< std::vector > &coefficients, double x, double z){ - this->output = z; - this->firstDim = x; - this->coefficients = coefficients; - this->dim = i2D; -} - -double PolyResidual::call(double x){ - double polyRes = -1; - switch (this->dim) { - case i1D: - polyRes = this->poly.polyval(this->coefficients[0], x); - break; - case i2D: - polyRes = this->poly.polyval(this->coefficients, this->firstDim, x); - break; - default: - throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); - } - return polyRes - this->output; -} - -double PolyResidual::deriv(double x){ - double polyRes = -1; - switch (this->dim) { - case i1D: - polyRes = this->poly.polyder(this->coefficients[0], x); - break; - case i2D: - polyRes = this->poly.polyder(this->coefficients, this->firstDim, x); - break; - default: - throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); - } - return polyRes; -} - -double PolyIntResidual::call(double x){ - double polyRes = -1; - switch (this->dim) { - case i1D: - polyRes = this->poly.polyint(this->coefficients[0], x); - break; - case i2D: - polyRes = this->poly.polyint(this->coefficients, this->firstDim, x); - break; - default: - throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); - } - return polyRes - this->output; -} - -double PolyIntResidual::deriv(double x){ - double polyRes = -1; - switch (this->dim) { - case i1D: - polyRes = this->poly.polyval(this->coefficients[0], x); - break; - case i2D: - polyRes = this->poly.polyval(this->coefficients, this->firstDim, x); - break; - default: - throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); - } - return polyRes; -} - -double PolyFracIntResidual::call(double x){ - double polyRes = -1; - switch (this->dim) { - case i1D: - polyRes = this->poly.polyfracint(this->coefficients[0], x); - break; - case i2D: - polyRes = this->poly.polyfracint(this->coefficients, this->firstDim, x); - break; - default: - throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); - } - return polyRes - this->output; -} - -double PolyFracIntResidual::deriv(double x){ - double polyRes = -1; - switch (this->dim) { - case i1D: - polyRes = this->poly.polyfracval(this->coefficients[0], x); - break; - case i2D: - polyRes = this->poly.polyfracval(this->coefficients, this->firstDim, x); - break; - default: - throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); - } - return polyRes; -} - -double PolyFracIntCentralResidual::call(double x){ - double polyRes = -1; - switch (this->dim) { - case i1D: - polyRes = this->poly.polyfracintcentral(this->coefficients[0], x, this->baseVal); - break; - case i2D: - polyRes = this->poly.polyfracintcentral(this->coefficients, this->firstDim, x, this->baseVal); - break; - default: - throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); - } - return polyRes - this->output; -} - -double PolyFracIntCentralResidual::deriv(double x){ - throw CoolProp::NotImplementedError("Derivative of a polynomial frac int is not defined."); -} - -double PolyDerResidual::call(double x){ - double polyRes = -1; - switch (this->dim) { - case i1D: - polyRes = this->poly.polyder(this->coefficients[0], x); - break; - case i2D: - polyRes = this->poly.polyder(this->coefficients, this->firstDim, x); - break; - default: - throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); - } - return polyRes - this->output; -} - -double PolyDerResidual::deriv(double x){ - throw CoolProp::NotImplementedError("2nd derivative of a polynomial is not defined."); -} - - - - -/** Implements the same public functions as the BasePolynomial - * but solves the polynomial for the given value - * instead of evaluating it. - * TODO: This class does not check for bijective - * polynomials and is therefore a little - * fragile. - */ -PolynomialSolver::PolynomialSolver(){ - this->POLYMATH_DEBUG = false; - this->macheps = DBL_EPSILON; - this->tol = DBL_EPSILON*1e3; - this->maxiter = 50; -} - -/** Everything related to the normal polynomials goes in this - * section, holds all the functions for solving polynomials. - */ -/// Solves a one-dimensional polynomial for the given coefficients -/// @param coefficients vector containing the ordered coefficients -/// @param y double value that represents the current input -double PolynomialSolver::polyval(const std::vector &coefficients, double y) { - PolyResidual residual = PolyResidual(coefficients, y); - return this->solve(residual); -} - -/// Solves a two-dimensional polynomial for the given coefficients -/// @param coefficients vector containing the ordered coefficients -/// @param x double value that represents the current input in the 1st dimension -/// @param z double value that represents the current output -double PolynomialSolver::polyval(const std::vector< std::vector > &coefficients, double x, double z){ - PolyResidual residual = PolyResidual(coefficients, x, z); - return this->solve(residual); -} - - -/** Everything related to the integrated polynomials goes in this - * section, holds all the functions for solving polynomials. - */ -/// Solves the indefinite integral of a one-dimensional polynomial -/// @param coefficients vector containing the ordered coefficients -/// @param y double value that represents the current output -double PolynomialSolver::polyint(const std::vector &coefficients, double y){ - PolyIntResidual residual = PolyIntResidual(coefficients, y); - return this->solve(residual); -} - -/// Solves the indefinite integral of a two-dimensional polynomial along the 2nd axis (y) -/// @param coefficients vector containing the ordered coefficients -/// @param x double value that represents the current input in the 1st dimension -/// @param z double value that represents the current output -double PolynomialSolver::polyint(const std::vector< std::vector > &coefficients, double x, double z){ - PolyIntResidual residual = PolyIntResidual(coefficients, x, z); - return this->solve(residual); -} - - -/** Everything related to the derived polynomials goes in this - * section, holds all the functions for solving polynomials. - */ -/// Solves the derivative of a one-dimensional polynomial -/// @param coefficients vector containing the ordered coefficients -/// @param y double value that represents the current output -double PolynomialSolver::polyder(const std::vector &coefficients, double y){ - PolyDerResidual residual = PolyDerResidual(coefficients, y); - return this->solve(residual); -} - -/// Solves the derivative of a two-dimensional polynomial along the 2nd axis (y) -/// @param coefficients vector containing the ordered coefficients -/// @param x double value that represents the current input in the 1st dimension -/// @param z double value that represents the current output -double PolynomialSolver::polyder(const std::vector< std::vector > &coefficients, double x, double z){ - PolyDerResidual residual = PolyDerResidual(coefficients, x, z); - return this->solve(residual); -} - - -/** Everything related to the polynomials divided by one variable goes in this - * section, holds all the functions for solving polynomials. - */ -/// Solves the indefinite integral of a one-dimensional polynomial divided by its independent variable -/// @param coefficients vector containing the ordered coefficients -/// @param y double value that represents the current output -double PolynomialSolver::polyfracval(const std::vector &coefficients, double y){ - throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function -} - -/// Solves the indefinite integral of a two-dimensional polynomial divided by its 2nd independent variable -/// @param coefficients vector containing the ordered coefficients -/// @param x double value that represents the current input in the 1st dimension -/// @param z double value that represents the current output -double PolynomialSolver::polyfracval(const std::vector< std::vector > &coefficients, double x, double z){ - throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function -} - - -/** Everything related to the integrated polynomials divided by one variable goes in this - * section, holds all the functions for solving polynomials. - */ -/// Solves the indefinite integral of a one-dimensional polynomial divided by its independent variable -/// @param coefficients vector containing the ordered coefficients -/// @param y double value that represents the current output -double PolynomialSolver::polyfracint(const std::vector &coefficients, double y){ - PolyFracIntResidual residual = PolyFracIntResidual(coefficients, y); - return this->solve(residual); -} - -/// Solves the indefinite integral of a two-dimensional polynomial divided by its 2nd independent variable -/// @param coefficients vector containing the ordered coefficients -/// @param x double value that represents the current input in the 1st dimension -/// @param z double value that represents the current output -double PolynomialSolver::polyfracint(const std::vector< std::vector > &coefficients, double x, double z){ - PolyFracIntResidual residual = PolyFracIntResidual(coefficients, x, z); - return this->solve(residual); -} - -/// Solves the indefinite integral of a centred one-dimensional polynomial divided by its independent variable -/// @param coefficients vector containing the ordered coefficients -/// @param y double value that represents the current output -/// @param xbase central x-value for fitted function -double PolynomialSolver::polyfracintcentral(const std::vector &coefficients, double y, double xbase){ - PolyFracIntCentralResidual residual = PolyFracIntCentralResidual(coefficients, y, xbase); - return this->solve(residual); -} - -/// Solves the indefinite integral of a centred two-dimensional polynomial divided by its 2nd independent variable -/// @param coefficients vector containing the ordered coefficients -/// @param x double value that represents the current input in the 1st dimension -/// @param z double value that represents the current output -/// @param ybase central y-value for fitted function -double PolynomialSolver::polyfracintcentral(const std::vector< std::vector > &coefficients, double x, double z, double ybase){ - PolyFracIntCentralResidual residual = PolyFracIntCentralResidual(coefficients, x, z, ybase); - return this->solve(residual); -} - - -/** Everything related to the derived polynomials divided by one variable goes in this - * section, holds all the functions for solving polynomials. - */ -/// Solves the derivative of a one-dimensional polynomial divided by its independent variable -/// @param coefficients vector containing the ordered coefficients -/// @param y double value that represents the current output -double PolynomialSolver::polyfracder(const std::vector &coefficients, double y){ - throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function -} - -/// Solves the derivative of a two-dimensional polynomial divided by its 2nd independent variable -/// @param coefficients vector containing the ordered coefficients -/// @param x double value that represents the current input in the 1st dimension -/// @param z double value that represents the current output -double PolynomialSolver::polyfracder(const std::vector< std::vector > &coefficients, double x, double z){ - throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function -} - -/// Solves the derivative of a centred one-dimensional polynomial divided by its independent variable -/// @param coefficients vector containing the ordered coefficients -/// @param y double value that represents the current output -/// @param xbase central x-value for fitted function -double PolynomialSolver::polyfracdercentral(const std::vector &coefficients, double y, double xbase){ - throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function -} - -/// Solves the derivative of a centred two-dimensional polynomial divided by its 2nd independent variable -/// @param coefficients vector containing the ordered coefficients -/// @param x double value that represents the current input in the 1st dimension -/// @param z double value that represents the current output -/// @param ybase central y-value for fitted function -double PolynomialSolver::polyfracdercentral(const std::vector< std::vector > &coefficients, double x, double z, double ybase){ - throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function -} - - -/** Set the solvers and updates either the guess values or the - * boundaries for the variable to solve for. - */ -/// Sets the guess value for the Newton solver and enables it. -/// @param guess double value that represents the guess value -void PolynomialSolver::setGuess(double guess){ - this->uses = iNewton; - this->guess = guess; - this->min = -1; - this->max = -1; -} -/// Sets the limits for the Brent solver and enables it. -/// @param min double value that represents the lower boundary -/// @param max double value that represents the upper boundary -void PolynomialSolver::setLimits(double min, double max){ - this->uses = iBrent; - this->guess = -1; - this->min = min; - this->max = max; -} - -/// Solves the equations based on previously defined parameters. -/// @param min double value that represents the lower boundary -/// @param max double value that represents the upper boundary -double PolynomialSolver::solve(PolyResidual &res){ - std::string errstring; - double result = -1.0; - switch (this->uses) { - case iNewton: ///< Newton solver with derivative and guess value - if (res.is2D()) { - throw CoolProp::NotImplementedError("The Newton solver is not suitable for 2D polynomials, yet."); - } - result = Newton(res, this->guess, this->tol, this->maxiter, errstring); - break; - - case iBrent: ///< Brent solver with bounds - result = Brent(res, this->min, this->max, this->macheps, this->tol, this->maxiter, errstring); - break; - - default: - throw CoolProp::NotImplementedError("This solver has not been implemented or you forgot to select a solver..."); - } - return result; -} - - -/** Here we define the functions that should be to evaluate exponential - * functions. Not really polynomials, I know... - */ - -BaseExponential::BaseExponential(){ - this->POLYMATH_DEBUG = false; -// this->poly = new BaseExponential(); -} -// -//BaseExponential::~BaseExponential(){ -// delete this->poly; +//bool BasePolynomial::checkCoefficients(std::vector< std::vector > const& coefficients, unsigned int rows, unsigned int columns){ +// if (coefficients.size() == rows){ +// bool result = true; +// for(unsigned int i=0; i &coefficients, double x, int n){ - double result = 0.; - if (n==1) { - this->poly.checkCoefficients(coefficients,3); - result = exp(coefficients[0]/(x+coefficients[1]) - coefficients[2]); - } else if (n==2) { - result = exp(this->poly.polyval(coefficients, x)); - } else { - throw ValueError(format("There is no function defined for this input (%d). ",n)); - } - return result; -} - -/// Evaluates an exponential function for the given coefficients -/// @param coefficients vector containing the ordered coefficients -/// @param x double value that represents the current input in the 1st dimension -/// @param y double value that represents the current input in the 2nd dimension -/// @param n int value that determines the kind of exponential function -double BaseExponential::expval(const std::vector< std::vector > &coefficients, double x, double y, int n){ - double result = 0.; - if (n==2) { - result = exp(this->poly.polyval(coefficients, x, y)); - } else { - throw ValueError(format("There is no function defined for this input (%d). ",n)); - } - return result; -} - - -} - - -#ifdef ENABLE_CATCH -#include -#include "catch.hpp" - -class PolynomialConsistencyFixture { -public: - CoolProp::BasePolynomial poly; - CoolProp::PolynomialSolver solver; -// enum dims {i1D, i2D}; -// double firstDim; -// int dim; -// std::vector< std::vector > coefficients; // -// void setInputs(const std::vector &coefficients){ -// this->firstDim = 0; -// this->coefficients.clear(); -// this->coefficients.push_back(coefficients); -// this->dim = i1D; -// } // -// void setInputs(const std::vector< std::vector > &coefficients, double x){ -// this->firstDim = x; -// this->coefficients = coefficients; -// this->dim = i2D; -// } -}; - - -TEST_CASE("Internal consistency checks with PolyMath objects","[PolyMath]") -{ - CoolProp::BasePolynomial poly; - CoolProp::PolynomialSolver solver; - - /// Test case for "SylthermXLT" by "Dow Chemicals" - std::vector cHeat; - cHeat.clear(); - cHeat.push_back(+1.1562261074E+03); - cHeat.push_back(+2.0994549103E+00); - cHeat.push_back(+7.7175381057E-07); - cHeat.push_back(-3.7008444051E-20); - - double deltaT = 0.1; - double Tmin = 273.15- 50; - double Tmax = 273.15+250; - double Tinc = 15; - - double val1,val2,val3,val4; - - SECTION("DerFromVal1D") { - for (double T = Tmin; T BasePolynomial::integrateCoeffs(std::vector const& coefficients){ +// std::vector newCoefficients; +// unsigned int sizeX = coefficients.size(); +// if (sizeX<1) throw ValueError(format("You have to provide coefficients, a vector length of %d is not a valid. ",sizeX)); +// // pushing a zero elevates the order by 1 +// newCoefficients.push_back(0.0); +// for(unsigned int i=0; i > BasePolynomial::integrateCoeffs(std::vector< std::vector > const& coefficients, bool axis){ +// std::vector< std::vector > newCoefficients; +// unsigned int sizeX = coefficients.size(); +// if (sizeX<1) throw ValueError(format("You have to provide coefficients, a vector length of %d is not a valid. ",sizeX)); // -// val1 = poly.polyfracint(cHeat, T); -// solver.setGuess(T+100); -// val2 = solver.polyfracint(cHeat, val1); -// CAPTURE(T); -// CAPTURE(val1); -// CAPTURE(val2); -// CHECK(fabs(T-val2) < 1e-1); - } - } - SECTION("Solve1DBrent") { - for (double T = Tmin; T > cHeat2D; - cHeat2D.clear(); - cHeat2D.push_back(cHeat); - cHeat2D.push_back(cHeat); - cHeat2D.push_back(cHeat); - - //setInputs(cHeat2D, 0.3); - - SECTION("DerFromVal2D") { - for (double T = Tmin; T > tmpCoefficients; +// tmpCoefficients = transpose(coefficients); +// unsigned int sizeY = tmpCoefficients.size(); +// for(unsigned int i=0; i BasePolynomial::deriveCoeffs(std::vector const& coefficients){ +// std::vector newCoefficients; +// unsigned int sizeX = coefficients.size(); +// if (sizeX<1) throw ValueError(format("You have to provide coefficients, a vector length of %d is not a valid. ",sizeX)); +// // skipping the first element lowers the order +// for(unsigned int i=1; i > BasePolynomial::deriveCoeffs(const std::vector< std::vector > &coefficients, unsigned int axis){ +// std::vector< std::vector > newCoefficients; +// unsigned int sizeX = coefficients.size(); +// if (sizeX<1) throw ValueError(format("You have to provide coefficients, a vector length of %d is not a valid. ",sizeX)); +// +// if (axis==0){ +// std::vector< std::vector > tmpCoefficients; +// tmpCoefficients = transpose(coefficients); +// unsigned int sizeY = tmpCoefficients.size(); +// for(unsigned int i=0; i const& coefficients, double x){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running simplePolynomial(std::vector, " << x << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// double result = 0.; +// for(unsigned int i=0; iPOLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +//double BasePolynomial::simplePolynomial(std::vector > const& coefficients, double x, double y){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running simplePolynomial(std::vector, " << x << ", " << y << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// double result = 0; +// for(unsigned int i=0; iPOLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +// +///// Simple integrated polynomial function generator. +///** Base function to produce integrals of n-th order +// * polynomials based on the length of the coefficient +// * vector. +// * Starts with only the first coefficient at x^0 */ +/////Indefinite integral in x-direction +//double BasePolynomial::simplePolynomialInt(std::vector const& coefficients, double x){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running simplePolynomialInt(std::vector, " << x << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// double result = 0.; +// for(unsigned int i=0; iPOLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +/////Indefinite integral in y-direction only +//double BasePolynomial::simplePolynomialInt(std::vector > const& coefficients, double x, double y){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running simplePolynomialInt(std::vector, " << x << ", " << y << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// double result = 0.; +// for(unsigned int i=0; iPOLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +// +///// Simple integrated polynomial function generator divided by independent variable. +///** Base function to produce integrals of n-th order +// * polynomials based on the length of the coefficient +// * vector. +// * Starts with only the first coefficient at x^0 */ +//double BasePolynomial::simpleFracInt(std::vector const& coefficients, double x){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running simpleFracInt(std::vector, " << x << "): "; +// } +// double result = coefficients[0] * log(x); +// if (coefficients.size() > 1) { +// for (unsigned int i=1; iPOLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +//double BasePolynomial::simpleFracInt(std::vector< std::vector > const& coefficients, double x, double y){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running simpleFracInt(std::vector, " << x << ", " << y << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// double result = 0; +// for (unsigned int i=0; iPOLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +// +///** Simple integrated centred(!) polynomial function generator divided by independent variable. +// * We need to rewrite some of the functions in order to +// * use central fit. Having a central temperature xbase +// * allows for a better fit, but requires a different +// * formulation of the fracInt function group. Other +// * functions are not affected. +// * Starts with only the first coefficient at x^0 */ +// +/////Helper functions to calculate binomial coefficients: http://rosettacode.org/wiki/Evaluate_binomial_coefficients#C.2B.2B +////double BasePolynomial::factorial(double nValue){ +//// double result = nValue; +//// double result_next; +//// double pc = nValue; +//// do { +//// result_next = result*(pc-1); +//// result = result_next; +//// pc--; +//// } while(pc>2); +//// nValue = result; +//// return nValue; +////} +////double BasePolynomial::factorial(double nValue){ +//// if (nValue == 0) return (1); +//// else return (nValue * factorial(nValue - 1)); +////} +//double BasePolynomial::factorial(double nValue){ +// double value = 1; +// for(int i = 2; i <= nValue; i++){ +// value = value * i; +// } +// return value; +//} +// +//double BasePolynomial::binom(double nValue, double nValue2){ +// double result; +// if(nValue2 == 1) return nValue; +// result = (factorial(nValue)) / (factorial(nValue2)*factorial((nValue - nValue2))); +// nValue2 = result; +// return nValue2; +//} +// +/////Helper functions to calculate the D vector: +//std::vector BasePolynomial::fracIntCentralDvector(int m, double x, double xbase){ +// std::vector D; +// double tmp; +// if (m<1) throw ValueError(format("You have to provide coefficients, a vector length of %d is not a valid. ",m)); +// for (int j=0; j const& coefficients, double x, double xbase){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running fracIntCentral(std::vector, " << x << ", " << xbase << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// int m = coefficients.size(); +// std::vector D = fracIntCentralDvector(m, x, xbase); +// double result = 0; +// for(int j=0; jPOLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +// +///// Horner function generator implementations +///** Represent polynomials according to Horner's scheme. +// * This avoids unnecessary multiplication and thus +// * speeds up calculation. +// */ +//double BasePolynomial::baseHorner(std::vector const& coefficients, double x){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running baseHorner(std::vector, " << x << "): "; +// } +// double result = 0; +// for(int i=coefficients.size()-1; i>=0; i--) { +// result *= x; +// result += coefficients[i]; +// } +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +//double BasePolynomial::baseHorner(std::vector< std::vector > const& coefficients, double x, double y){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running baseHorner(std::vector, " << x << ", " << y << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// double result = 0; +// for(int i=coefficients.size()-1; i>=0; i--) { +// result *= x; +// result += baseHorner(coefficients[i], y); +// } +// this->POLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +/////Indefinite integral in x-direction +//double BasePolynomial::baseHornerInt(std::vector const& coefficients, double x){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running baseHornerInt(std::vector, " << x << "): "; +// } +// double result = 0; +// for(int i=coefficients.size()-1; i>=0; i--) { +// result *= x; +// result += coefficients[i]/(i+1.); +// } +// result = result * x; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +/////Indefinite integral in y-direction only +//double BasePolynomial::baseHornerInt(std::vector > const& coefficients, double x, double y){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running baseHornerInt(std::vector, " << x << ", " << y << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// double result = 0; +// for(int i=coefficients.size()-1; i>=0; i--) { +// result *= x; +// result += baseHornerInt(coefficients[i], y); +// } +// this->POLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +/////Indefinite integral in x-direction of a polynomial divided by its independent variable +//double BasePolynomial::baseHornerFracInt(std::vector const& coefficients, double x){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running baseHornerFra(std::vector, " << x << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// double result = 0; +// if (coefficients.size() > 1) { +// for(int i=coefficients.size()-1; i>=1; i--) { +// result *= x; +// result += coefficients[i]/(i); +// } +// result *= x; +// } +// result += coefficients[0] * log(x); +// this->POLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +/////Indefinite integral in y-direction of a polynomial divided by its 2nd independent variable +//double BasePolynomial::baseHornerFracInt(std::vector > const& coefficients, double x, double y){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running baseHornerFra(std::vector, " << x << ", " << y << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// +// double result = 0; +// for(int i=coefficients.size()-1; i>=0; i--) { +// result *= x; +// result += baseHornerFracInt(coefficients[i], y); +// } +// +// this->POLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +// +///** Alternatives +// * Simple functions that heavily rely on other parts of this file. +// * We still need to check which combinations yield the best +// * performance. +// */ +/////Derivative in x-direction +//double BasePolynomial::deriveIn2Steps(std::vector const& coefficients, double x){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running deriveIn2Steps(std::vector, " << x << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// double result = polyval(deriveCoeffs(coefficients),x); +// this->POLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +/////Derivative in terms of x(axis=true) or y(axis=false). +//double BasePolynomial::deriveIn2Steps(std::vector< std::vector > const& coefficients, double x, double y, bool axis){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running deriveIn2Steps(std::vector, " << x << ", " << y << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// double result = polyval(deriveCoeffs(coefficients,axis),x,y); +// this->POLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +/////Indefinite integral in x-direction +//double BasePolynomial::integrateIn2Steps(std::vector const& coefficients, double x){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running integrateIn2Steps(std::vector, " << x << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// double result = polyval(integrateCoeffs(coefficients),x); +// this->POLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +/////Indefinite integral in terms of x(axis=true) or y(axis=false). +//double BasePolynomial::integrateIn2Steps(std::vector< std::vector > const& coefficients, double x, double y, bool axis){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running integrateIn2Steps(std::vector, " << x << ", " << y << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// double result = polyval(integrateCoeffs(coefficients,axis),x,y); +// this->POLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +/////Indefinite integral in x-direction of a polynomial divided by its independent variable +//double BasePolynomial::fracIntIn2Steps(std::vector const& coefficients, double x){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running fracIntIn2Steps(std::vector, " << x << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// double result = coefficients[0] * log(x); +// if (coefficients.size() > 1) { +// std::vector newCoeffs(coefficients.begin() + 1, coefficients.end()); +// result += polyint(newCoeffs,x); +// } +// this->POLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +/////Indefinite integral in y-direction of a polynomial divided by its 2nd independent variable +//double BasePolynomial::fracIntIn2Steps(std::vector > const& coefficients, double x, double y){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running fracIntIn2Steps(std::vector, " << x << ", " << y << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// std::vector newCoeffs; +// for (unsigned int i=0; iPOLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +/////Indefinite integral in y-direction of a centred polynomial divided by its 2nd independent variable +//double BasePolynomial::fracIntCentral2Steps(std::vector > const& coefficients, double x, double y, double ybase){ +// if (this->POLYMATH_DEBUG) { +// std::cout << "Running fracIntCentral2Steps(std::vector, " << x << ", " << y << ", " << ybase << "): "; +// } +// bool db = this->POLYMATH_DEBUG; +// this->POLYMATH_DEBUG = false; +// std::vector newCoeffs; +// for (unsigned int i=0; iPOLYMATH_DEBUG = db; +// if (this->POLYMATH_DEBUG) { +// std::cout << result << std::endl; +// } +// return result; +//} +// +// +// +// +///** Implements the function wrapper interface and can be +// * used by the solvers. This is only an example and you should +// * use local redefinitions of the class. +// * TODO: Make multidimensional +// */ +//PolyResidual::PolyResidual(){ +// this->dim = -1; +//} +// +//PolyResidual::PolyResidual(const std::vector &coefficients, double y){ +// this->output = y; +// this->firstDim = 0; +// this->coefficients.clear(); +// this->coefficients.push_back(coefficients); +// this->dim = i1D; +//} +// +//PolyResidual::PolyResidual(const std::vector< std::vector > &coefficients, double x, double z){ +// this->output = z; +// this->firstDim = x; +// this->coefficients = coefficients; +// this->dim = i2D; +//} +// +//double PolyResidual::call(double x){ +// double polyRes = -1; +// switch (this->dim) { +// case i1D: +// polyRes = this->poly.polyval(this->coefficients[0], x); +// break; +// case i2D: +// polyRes = this->poly.polyval(this->coefficients, this->firstDim, x); +// break; +// default: +// throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); +// } +// return polyRes - this->output; +//} +// +//double PolyResidual::deriv(double x){ +// double polyRes = -1; +// switch (this->dim) { +// case i1D: +// polyRes = this->poly.polyder(this->coefficients[0], x); +// break; +// case i2D: +// polyRes = this->poly.polyder(this->coefficients, this->firstDim, x); +// break; +// default: +// throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); +// } +// return polyRes; +//} +// +//double PolyIntResidual::call(double x){ +// double polyRes = -1; +// switch (this->dim) { +// case i1D: +// polyRes = this->poly.polyint(this->coefficients[0], x); +// break; +// case i2D: +// polyRes = this->poly.polyint(this->coefficients, this->firstDim, x); +// break; +// default: +// throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); +// } +// return polyRes - this->output; +//} +// +//double PolyIntResidual::deriv(double x){ +// double polyRes = -1; +// switch (this->dim) { +// case i1D: +// polyRes = this->poly.polyval(this->coefficients[0], x); +// break; +// case i2D: +// polyRes = this->poly.polyval(this->coefficients, this->firstDim, x); +// break; +// default: +// throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); +// } +// return polyRes; +//} +// +//double PolyFracIntResidual::call(double x){ +// double polyRes = -1; +// switch (this->dim) { +// case i1D: +// polyRes = this->poly.polyfracint(this->coefficients[0], x); +// break; +// case i2D: +// polyRes = this->poly.polyfracint(this->coefficients, this->firstDim, x); +// break; +// default: +// throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); +// } +// return polyRes - this->output; +//} +// +//double PolyFracIntResidual::deriv(double x){ +// double polyRes = -1; +// switch (this->dim) { +// case i1D: +// polyRes = this->poly.polyfracval(this->coefficients[0], x); +// break; +// case i2D: +// polyRes = this->poly.polyfracval(this->coefficients, this->firstDim, x); +// break; +// default: +// throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); +// } +// return polyRes; +//} +// +//double PolyFracIntCentralResidual::call(double x){ +// double polyRes = -1; +// switch (this->dim) { +// case i1D: +// polyRes = this->poly.polyfracintcentral(this->coefficients[0], x, this->baseVal); +// break; +// case i2D: +// polyRes = this->poly.polyfracintcentral(this->coefficients, this->firstDim, x, this->baseVal); +// break; +// default: +// throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); +// } +// return polyRes - this->output; +//} +// +//double PolyFracIntCentralResidual::deriv(double x){ +// throw CoolProp::NotImplementedError("Derivative of a polynomial frac int is not defined."); +//} +// +//double PolyDerResidual::call(double x){ +// double polyRes = -1; +// switch (this->dim) { +// case i1D: +// polyRes = this->poly.polyder(this->coefficients[0], x); +// break; +// case i2D: +// polyRes = this->poly.polyder(this->coefficients, this->firstDim, x); +// break; +// default: +// throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); +// } +// return polyRes - this->output; +//} +// +//double PolyDerResidual::deriv(double x){ +// throw CoolProp::NotImplementedError("2nd derivative of a polynomial is not defined."); +//} +// +// +// +// +///** Implements the same public functions as the BasePolynomial +// * but solves the polynomial for the given value +// * instead of evaluating it. +// * TODO: This class does not check for bijective +// * polynomials and is therefore a little +// * fragile. +// */ +//PolynomialSolver::PolynomialSolver(){ +// this->POLYMATH_DEBUG = false; +// this->macheps = DBL_EPSILON; +// this->tol = DBL_EPSILON*1e3; +// this->maxiter = 50; +//} +// +///** Everything related to the normal polynomials goes in this +// * section, holds all the functions for solving polynomials. +// */ +///// Solves a one-dimensional polynomial for the given coefficients +///// @param coefficients vector containing the ordered coefficients +///// @param y double value that represents the current input +//double PolynomialSolver::polyval(const std::vector &coefficients, double y) { +// PolyResidual residual = PolyResidual(coefficients, y); +// return this->solve(residual); +//} +// +///// Solves a two-dimensional polynomial for the given coefficients +///// @param coefficients vector containing the ordered coefficients +///// @param x double value that represents the current input in the 1st dimension +///// @param z double value that represents the current output +//double PolynomialSolver::polyval(const std::vector< std::vector > &coefficients, double x, double z){ +// PolyResidual residual = PolyResidual(coefficients, x, z); +// return this->solve(residual); +//} +// +// +///** Everything related to the integrated polynomials goes in this +// * section, holds all the functions for solving polynomials. +// */ +///// Solves the indefinite integral of a one-dimensional polynomial +///// @param coefficients vector containing the ordered coefficients +///// @param y double value that represents the current output +//double PolynomialSolver::polyint(const std::vector &coefficients, double y){ +// PolyIntResidual residual = PolyIntResidual(coefficients, y); +// return this->solve(residual); +//} +// +///// Solves the indefinite integral of a two-dimensional polynomial along the 2nd axis (y) +///// @param coefficients vector containing the ordered coefficients +///// @param x double value that represents the current input in the 1st dimension +///// @param z double value that represents the current output +//double PolynomialSolver::polyint(const std::vector< std::vector > &coefficients, double x, double z){ +// PolyIntResidual residual = PolyIntResidual(coefficients, x, z); +// return this->solve(residual); +//} +// +// +///** Everything related to the derived polynomials goes in this +// * section, holds all the functions for solving polynomials. +// */ +///// Solves the derivative of a one-dimensional polynomial +///// @param coefficients vector containing the ordered coefficients +///// @param y double value that represents the current output +//double PolynomialSolver::polyder(const std::vector &coefficients, double y){ +// PolyDerResidual residual = PolyDerResidual(coefficients, y); +// return this->solve(residual); +//} +// +///// Solves the derivative of a two-dimensional polynomial along the 2nd axis (y) +///// @param coefficients vector containing the ordered coefficients +///// @param x double value that represents the current input in the 1st dimension +///// @param z double value that represents the current output +//double PolynomialSolver::polyder(const std::vector< std::vector > &coefficients, double x, double z){ +// PolyDerResidual residual = PolyDerResidual(coefficients, x, z); +// return this->solve(residual); +//} +// +// +///** Everything related to the polynomials divided by one variable goes in this +// * section, holds all the functions for solving polynomials. +// */ +///// Solves the indefinite integral of a one-dimensional polynomial divided by its independent variable +///// @param coefficients vector containing the ordered coefficients +///// @param y double value that represents the current output +//double PolynomialSolver::polyfracval(const std::vector &coefficients, double y){ +// throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function +//} +// +///// Solves the indefinite integral of a two-dimensional polynomial divided by its 2nd independent variable +///// @param coefficients vector containing the ordered coefficients +///// @param x double value that represents the current input in the 1st dimension +///// @param z double value that represents the current output +//double PolynomialSolver::polyfracval(const std::vector< std::vector > &coefficients, double x, double z){ +// throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function +//} +// +// +///** Everything related to the integrated polynomials divided by one variable goes in this +// * section, holds all the functions for solving polynomials. +// */ +///// Solves the indefinite integral of a one-dimensional polynomial divided by its independent variable +///// @param coefficients vector containing the ordered coefficients +///// @param y double value that represents the current output +//double PolynomialSolver::polyfracint(const std::vector &coefficients, double y){ +// PolyFracIntResidual residual = PolyFracIntResidual(coefficients, y); +// return this->solve(residual); +//} +// +///// Solves the indefinite integral of a two-dimensional polynomial divided by its 2nd independent variable +///// @param coefficients vector containing the ordered coefficients +///// @param x double value that represents the current input in the 1st dimension +///// @param z double value that represents the current output +//double PolynomialSolver::polyfracint(const std::vector< std::vector > &coefficients, double x, double z){ +// PolyFracIntResidual residual = PolyFracIntResidual(coefficients, x, z); +// return this->solve(residual); +//} +// +///// Solves the indefinite integral of a centred one-dimensional polynomial divided by its independent variable +///// @param coefficients vector containing the ordered coefficients +///// @param y double value that represents the current output +///// @param xbase central x-value for fitted function +//double PolynomialSolver::polyfracintcentral(const std::vector &coefficients, double y, double xbase){ +// PolyFracIntCentralResidual residual = PolyFracIntCentralResidual(coefficients, y, xbase); +// return this->solve(residual); +//} +// +///// Solves the indefinite integral of a centred two-dimensional polynomial divided by its 2nd independent variable +///// @param coefficients vector containing the ordered coefficients +///// @param x double value that represents the current input in the 1st dimension +///// @param z double value that represents the current output +///// @param ybase central y-value for fitted function +//double PolynomialSolver::polyfracintcentral(const std::vector< std::vector > &coefficients, double x, double z, double ybase){ +// PolyFracIntCentralResidual residual = PolyFracIntCentralResidual(coefficients, x, z, ybase); +// return this->solve(residual); +//} +// +// +///** Everything related to the derived polynomials divided by one variable goes in this +// * section, holds all the functions for solving polynomials. +// */ +///// Solves the derivative of a one-dimensional polynomial divided by its independent variable +///// @param coefficients vector containing the ordered coefficients +///// @param y double value that represents the current output +//double PolynomialSolver::polyfracder(const std::vector &coefficients, double y){ +// throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function +//} +// +///// Solves the derivative of a two-dimensional polynomial divided by its 2nd independent variable +///// @param coefficients vector containing the ordered coefficients +///// @param x double value that represents the current input in the 1st dimension +///// @param z double value that represents the current output +//double PolynomialSolver::polyfracder(const std::vector< std::vector > &coefficients, double x, double z){ +// throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function +//} +// +///// Solves the derivative of a centred one-dimensional polynomial divided by its independent variable +///// @param coefficients vector containing the ordered coefficients +///// @param y double value that represents the current output +///// @param xbase central x-value for fitted function +//double PolynomialSolver::polyfracdercentral(const std::vector &coefficients, double y, double xbase){ +// throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function +//} +// +///// Solves the derivative of a centred two-dimensional polynomial divided by its 2nd independent variable +///// @param coefficients vector containing the ordered coefficients +///// @param x double value that represents the current input in the 1st dimension +///// @param z double value that represents the current output +///// @param ybase central y-value for fitted function +//double PolynomialSolver::polyfracdercentral(const std::vector< std::vector > &coefficients, double x, double z, double ybase){ +// throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function +//} +// +// +///** Set the solvers and updates either the guess values or the +// * boundaries for the variable to solve for. +// */ +///// Sets the guess value for the Newton solver and enables it. +///// @param guess double value that represents the guess value +//void PolynomialSolver::setGuess(double guess){ +// this->uses = iNewton; +// this->guess = guess; +// this->min = -1; +// this->max = -1; +//} +///// Sets the limits for the Brent solver and enables it. +///// @param min double value that represents the lower boundary +///// @param max double value that represents the upper boundary +//void PolynomialSolver::setLimits(double min, double max){ +// this->uses = iBrent; +// this->guess = -1; +// this->min = min; +// this->max = max; +//} +// +///// Solves the equations based on previously defined parameters. +///// @param min double value that represents the lower boundary +///// @param max double value that represents the upper boundary +//double PolynomialSolver::solve(PolyResidual &res){ +// std::string errstring; +// double result = -1.0; +// switch (this->uses) { +// case iNewton: ///< Newton solver with derivative and guess value +// if (res.is2D()) { +// throw CoolProp::NotImplementedError("The Newton solver is not suitable for 2D polynomials, yet."); +// } +// result = Newton(res, this->guess, this->tol, this->maxiter, errstring); +// break; +// +// case iBrent: ///< Brent solver with bounds +// result = Brent(res, this->min, this->max, this->macheps, this->tol, this->maxiter, errstring); +// break; +// +// default: +// throw CoolProp::NotImplementedError("This solver has not been implemented or you forgot to select a solver..."); +// } +// return result; +//} +// +// +///** Here we define the functions that should be to evaluate exponential +// * functions. Not really polynomials, I know... +// */ +// +//BaseExponential::BaseExponential(){ +// this->POLYMATH_DEBUG = false; +//// this->poly = new BaseExponential(); +//} +//// +////BaseExponential::~BaseExponential(){ +//// delete this->poly; +////} +// +///// Evaluates an exponential function for the given coefficients +///// @param coefficients vector containing the ordered coefficients +///// @param x double value that represents the current input +///// @param n int value that determines the kind of exponential function +//double BaseExponential::expval(const std::vector &coefficients, double x, int n){ +// double result = 0.; +// if (n==1) { +// this->poly.checkCoefficients(coefficients,3); +// result = exp(coefficients[0]/(x+coefficients[1]) - coefficients[2]); +// } else if (n==2) { +// result = exp(this->poly.polyval(coefficients, x)); +// } else { +// throw ValueError(format("There is no function defined for this input (%d). ",n)); +// } +// return result; +//} +// +///// Evaluates an exponential function for the given coefficients +///// @param coefficients vector containing the ordered coefficients +///// @param x double value that represents the current input in the 1st dimension +///// @param y double value that represents the current input in the 2nd dimension +///// @param n int value that determines the kind of exponential function +//double BaseExponential::expval(const std::vector< std::vector > &coefficients, double x, double y, int n){ +// double result = 0.; +// if (n==2) { +// result = exp(this->poly.polyval(coefficients, x, y)); +// } else { +// throw ValueError(format("There is no function defined for this input (%d). ",n)); +// } +// return result; +//} +// +// +//} +// +// +//#ifdef ENABLE_CATCH +//#include +//#include "catch.hpp" +// +//class PolynomialConsistencyFixture { +//public: +// CoolProp::BasePolynomial poly; +// CoolProp::PolynomialSolver solver; +//// enum dims {i1D, i2D}; +//// double firstDim; +//// int dim; +//// std::vector< std::vector > coefficients; +//// +//// void setInputs(const std::vector &coefficients){ +//// this->firstDim = 0; +//// this->coefficients.clear(); +//// this->coefficients.push_back(coefficients); +//// this->dim = i1D; +//// } +//// +//// void setInputs(const std::vector< std::vector > &coefficients, double x){ +//// this->firstDim = x; +//// this->coefficients = coefficients; +//// this->dim = i2D; +//// } +//}; +// +// +//TEST_CASE("Internal consistency checks with PolyMath objects","[PolyMath]") //{ +// CoolProp::BasePolynomial poly; +// CoolProp::PolynomialSolver solver; +// // /// Test case for "SylthermXLT" by "Dow Chemicals" // std::vector cHeat; // cHeat.clear(); @@ -1250,29 +1056,31 @@ TEST_CASE("Internal consistency checks with PolyMath objects","[PolyMath]") // cHeat.push_back(+7.7175381057E-07); // cHeat.push_back(-3.7008444051E-20); // -// //setInputs(cHeat); // double deltaT = 0.1; +// double Tmin = 273.15- 50; +// double Tmax = 273.15+250; +// double Tinc = 15; +// // double val1,val2,val3,val4; // -// SECTION("DerFromVal1D") { -// for (double T = 273.15-50; T<273.15+250; T+=15) { -// val1 = this->poly.polyval(cHeat, T-deltaT); -// val2 = this->poly.polyval(cHeat, T+deltaT); +// SECTION("DerFromVal1D") { +// for (double T = Tmin; Tpoly.polyder(cHeat, T); +// val4 = poly.polyder(cHeat, T); // CAPTURE(T); // CAPTURE(val3); // CAPTURE(val4); // CHECK( (1.0-fabs(val4/val3)) < 1e-1); // } -// } -// -// SECTION("ValFromInt1D") { -// for (double T = 273.15-50; T<273.15+250; T+=15) { -// val1 = this->poly.polyint(cHeat, T-deltaT); -// val2 = this->poly.polyint(cHeat, T+deltaT); +// } +// SECTION("ValFromInt1D") { +// for (double T = Tmin; Tpoly.polyval(cHeat, T); +// val4 = poly.polyval(cHeat, T); // CAPTURE(T); // CAPTURE(val3); // CAPTURE(val4); @@ -1280,31 +1088,89 @@ TEST_CASE("Internal consistency checks with PolyMath objects","[PolyMath]") // } // } // -// SECTION("Solve1DNewton") { -// for (double T = 273.15-50; T<273.15+250; T+=15) { -// val1 = this->poly.polyval(cHeat, T); -// this->solver.setGuess(T+100); -// val2 = this->solver.polyval(cHeat, val1); +// SECTION("Solve1DNewton") { +// for (double T = Tmin; Tpoly.polyval(cHeat, T); -// this->solver.setLimits(T-300,T+300); -// val2 = this->solver.polyval(cHeat, val1); +// SECTION("Solve1DBrent") { +// for (double T = Tmin; T > cHeat2D; +// /// Test case for 2D +// double xDim1 = 0.3; +// std::vector< std::vector > cHeat2D; // cHeat2D.clear(); // cHeat2D.push_back(cHeat); // cHeat2D.push_back(cHeat); @@ -1313,11 +1179,11 @@ TEST_CASE("Internal consistency checks with PolyMath objects","[PolyMath]") // //setInputs(cHeat2D, 0.3); // // SECTION("DerFromVal2D") { -// for (double T = 273.15-50; T<273.15+250; T+=15) { -// val1 = this->poly.polyval(cHeat, T-deltaT); -// val2 = this->poly.polyval(cHeat, T+deltaT); +// for (double T = Tmin; Tpoly.polyder(cHeat, T); +// val4 = poly.polyder(cHeat2D, xDim1, T); // CAPTURE(T); // CAPTURE(val3); // CAPTURE(val4); @@ -1326,11 +1192,11 @@ TEST_CASE("Internal consistency checks with PolyMath objects","[PolyMath]") // } // // SECTION("ValFromInt2D") { -// for (double T = 273.15-50; T<273.15+250; T+=15) { -// val1 = this->poly.polyint(cHeat, T-deltaT); -// val2 = this->poly.polyint(cHeat, T+deltaT); +// for (double T = Tmin; Tpoly.polyval(cHeat, T); +// val4 = poly.polyval(cHeat2D, xDim1, T); // CAPTURE(T); // CAPTURE(val3); // CAPTURE(val4); @@ -1338,22 +1204,54 @@ TEST_CASE("Internal consistency checks with PolyMath objects","[PolyMath]") // } // } // -// SECTION("Solve2DNewton") { -// for (double T = 273.15-50; T<273.15+250; T+=15) { -// val1 = this->poly.polyval(cHeat, T); -// this->solver.setGuess(T+100); -// val2 = this->solver.polyval(cHeat, val1); +//// SECTION("Solve2DNewton") { +//// for (double T = Tmin; Tpoly.polyval(cHeat, T); -// this->solver.setLimits(T-300,T+300); -// val2 = this->solver.polyval(cHeat, val1); +// +// val1 = poly.polyint(cHeat2D, xDim1, T); +// solver.setLimits(T-300,T+300); +// val2 = solver.polyint(cHeat2D, xDim1, val1); +// CAPTURE(T); +// CAPTURE(val1); +// CAPTURE(val2); +// CHECK(fabs(T-val2) < 1e-1); +// +// val1 = poly.polyder(cHeat2D, xDim1, T); +// solver.setLimits(T-300,T+300); +// val2 = solver.polyder(cHeat2D, xDim1, val1); +// CAPTURE(T); +// CAPTURE(val1); +// CAPTURE(val2); +// CHECK(fabs(T-val2) < 1e-1); +// +// val1 = poly.polyfracint(cHeat2D, xDim1, T); +// solver.setLimits(T-100,T+100); +// val2 = solver.polyfracint(cHeat2D, xDim1, val1); +// CAPTURE(T); +// CAPTURE(val1); +// CAPTURE(val2); +// CHECK(fabs(T-val2) < 1e-1); +// +// val1 = poly.polyfracintcentral(cHeat2D, xDim1, T, 250); +// solver.setLimits(T-100,T+100); +// val2 = solver.polyfracintcentral(cHeat2D, xDim1, val1, 250); // CAPTURE(T); // CAPTURE(val1); // CAPTURE(val2); @@ -1363,57 +1261,180 @@ TEST_CASE("Internal consistency checks with PolyMath objects","[PolyMath]") // //} // -//TEST_CASE("Check against hard coded data","[PolyMath]") -//{ -// CHECK(fabs(HumidAir::f_factor(-60+273.15,101325)/(1.00708)-1) < 1e-3); -// CHECK(fabs(HumidAir::f_factor( 80+273.15,101325)/(1.00573)-1) < 1e-3); -// CHECK(fabs(HumidAir::f_factor(-60+273.15,10000e3)/(2.23918)-1) < 1e-3); -// CHECK(fabs(HumidAir::f_factor(300+273.15,10000e3)/(1.04804)-1) < 1e-3); -//} - - - -//int main() { +////TEST_CASE_METHOD(PolynomialConsistencyFixture,"Internal consistency checks","[PolyMath]") +////{ +//// /// Test case for "SylthermXLT" by "Dow Chemicals" +//// std::vector cHeat; +//// cHeat.clear(); +//// cHeat.push_back(+1.1562261074E+03); +//// cHeat.push_back(+2.0994549103E+00); +//// cHeat.push_back(+7.7175381057E-07); +//// cHeat.push_back(-3.7008444051E-20); +//// +//// //setInputs(cHeat); +//// double deltaT = 0.1; +//// double val1,val2,val3,val4; +//// +//// SECTION("DerFromVal1D") { +//// for (double T = 273.15-50; T<273.15+250; T+=15) { +//// val1 = this->poly.polyval(cHeat, T-deltaT); +//// val2 = this->poly.polyval(cHeat, T+deltaT); +//// val3 = (val2-val1)/2/deltaT; +//// val4 = this->poly.polyder(cHeat, T); +//// CAPTURE(T); +//// CAPTURE(val3); +//// CAPTURE(val4); +//// CHECK( (1.0-fabs(val4/val3)) < 1e-1); +//// } +//// } +//// +//// SECTION("ValFromInt1D") { +//// for (double T = 273.15-50; T<273.15+250; T+=15) { +//// val1 = this->poly.polyint(cHeat, T-deltaT); +//// val2 = this->poly.polyint(cHeat, T+deltaT); +//// val3 = (val2-val1)/2/deltaT; +//// val4 = this->poly.polyval(cHeat, T); +//// CAPTURE(T); +//// CAPTURE(val3); +//// CAPTURE(val4); +//// CHECK( (1.0-fabs(val4/val3)) < 1e-1); +//// } +//// } +//// +//// SECTION("Solve1DNewton") { +//// for (double T = 273.15-50; T<273.15+250; T+=15) { +//// val1 = this->poly.polyval(cHeat, T); +//// this->solver.setGuess(T+100); +//// val2 = this->solver.polyval(cHeat, val1); +//// CAPTURE(T); +//// CAPTURE(val1); +//// CAPTURE(val2); +//// CHECK(fabs(T-val2) < 1e-1); +//// } +//// } +//// SECTION("Solve1DBrent") { +//// for (double T = 273.15-50; T<273.15+250; T+=15) { +//// val1 = this->poly.polyval(cHeat, T); +//// this->solver.setLimits(T-300,T+300); +//// val2 = this->solver.polyval(cHeat, val1); +//// CAPTURE(T); +//// CAPTURE(val1); +//// CAPTURE(val2); +//// CHECK(fabs(T-val2) < 1e-1); +//// } +//// } +//// +//// /// Test case for 2D +//// std::vector< std::vector > cHeat2D; +//// cHeat2D.clear(); +//// cHeat2D.push_back(cHeat); +//// cHeat2D.push_back(cHeat); +//// cHeat2D.push_back(cHeat); +//// +//// //setInputs(cHeat2D, 0.3); +//// +//// SECTION("DerFromVal2D") { +//// for (double T = 273.15-50; T<273.15+250; T+=15) { +//// val1 = this->poly.polyval(cHeat, T-deltaT); +//// val2 = this->poly.polyval(cHeat, T+deltaT); +//// val3 = (val2-val1)/2/deltaT; +//// val4 = this->poly.polyder(cHeat, T); +//// CAPTURE(T); +//// CAPTURE(val3); +//// CAPTURE(val4); +//// CHECK( (1.0-fabs(val4/val3)) < 1e-1); +//// } +//// } +//// +//// SECTION("ValFromInt2D") { +//// for (double T = 273.15-50; T<273.15+250; T+=15) { +//// val1 = this->poly.polyint(cHeat, T-deltaT); +//// val2 = this->poly.polyint(cHeat, T+deltaT); +//// val3 = (val2-val1)/2/deltaT; +//// val4 = this->poly.polyval(cHeat, T); +//// CAPTURE(T); +//// CAPTURE(val3); +//// CAPTURE(val4); +//// CHECK( (1.0-fabs(val4/val3)) < 1e-1); +//// } +//// } +//// +//// SECTION("Solve2DNewton") { +//// for (double T = 273.15-50; T<273.15+250; T+=15) { +//// val1 = this->poly.polyval(cHeat, T); +//// this->solver.setGuess(T+100); +//// val2 = this->solver.polyval(cHeat, val1); +//// CAPTURE(T); +//// CAPTURE(val1); +//// CAPTURE(val2); +//// CHECK(fabs(T-val2) < 1e-1); +//// } +//// } +//// SECTION("Solve2DBrent") { +//// for (double T = 273.15-50; T<273.15+250; T+=15) { +//// val1 = this->poly.polyval(cHeat, T); +//// this->solver.setLimits(T-300,T+300); +//// val2 = this->solver.polyval(cHeat, val1); +//// CAPTURE(T); +//// CAPTURE(val1); +//// CAPTURE(val2); +//// CHECK(fabs(T-val2) < 1e-1); +//// } +//// } +//// +////} +//// +////TEST_CASE("Check against hard coded data","[PolyMath]") +////{ +//// CHECK(fabs(HumidAir::f_factor(-60+273.15,101325)/(1.00708)-1) < 1e-3); +//// CHECK(fabs(HumidAir::f_factor( 80+273.15,101325)/(1.00573)-1) < 1e-3); +//// CHECK(fabs(HumidAir::f_factor(-60+273.15,10000e3)/(2.23918)-1) < 1e-3); +//// CHECK(fabs(HumidAir::f_factor(300+273.15,10000e3)/(1.04804)-1) < 1e-3); +////} // -// Catch::ConfigData &config = session.configData(); -// config.testsOrTags.clear(); -// config.testsOrTags.push_back("[fast]"); -// session.useConfigData(config); -// return session.run(); // -//} - -#endif /* CATCH_ENABLED */ - - -//int main() { // -// std::vector cHeat; -// cHeat.clear(); -// cHeat.push_back(+1.1562261074E+03); -// cHeat.push_back(+2.0994549103E+00); -// cHeat.push_back(+7.7175381057E-07); -// cHeat.push_back(-3.7008444051E-20); +////int main() { +//// +//// Catch::ConfigData &config = session.configData(); +//// config.testsOrTags.clear(); +//// config.testsOrTags.push_back("[fast]"); +//// session.useConfigData(config); +//// return session.run(); +//// +////} // -// CoolProp::BasePolynomial base = CoolProp::BasePolynomial(); -// CoolProp::PolynomialSolver solve = CoolProp::PolynomialSolver(); +//#endif /* CATCH_ENABLED */ // -// double T = 273.15+50; // -// double c = base.polyval(cHeat,T); -// printf("Should be : c = %3.3f \t J/kg/K \n",1834.746); -// printf("From object: c = %3.3f \t J/kg/K \n",c); -// -// T = 0.0; -// solve.setGuess(75+273.15); -// T = solve.polyval(cHeat,c); -// printf("Should be : T = %3.3f \t K \n",273.15+50.0); -// printf("From object: T = %3.3f \t K \n",T); -// -// T = 0.0; -// solve.setLimits(273.15+10,273.15+100); -// T = solve.polyval(cHeat,c); -// printf("Should be : T = %3.3f \t K \n",273.15+50.0); -// printf("From object: T = %3.3f \t K \n",T); -// -//} +////int main() { +//// +//// std::vector cHeat; +//// cHeat.clear(); +//// cHeat.push_back(+1.1562261074E+03); +//// cHeat.push_back(+2.0994549103E+00); +//// cHeat.push_back(+7.7175381057E-07); +//// cHeat.push_back(-3.7008444051E-20); +//// +//// CoolProp::BasePolynomial base = CoolProp::BasePolynomial(); +//// CoolProp::PolynomialSolver solve = CoolProp::PolynomialSolver(); +//// +//// double T = 273.15+50; +//// +//// double c = base.polyval(cHeat,T); +//// printf("Should be : c = %3.3f \t J/kg/K \n",1834.746); +//// printf("From object: c = %3.3f \t J/kg/K \n",c); +//// +//// T = 0.0; +//// solve.setGuess(75+273.15); +//// T = solve.polyval(cHeat,c); +//// printf("Should be : T = %3.3f \t K \n",273.15+50.0); +//// printf("From object: T = %3.3f \t K \n",T); +//// +//// T = 0.0; +//// solve.setLimits(273.15+10,273.15+100); +//// T = solve.polyval(cHeat,c); +//// printf("Should be : T = %3.3f \t K \n",273.15+50.0); +//// printf("From object: T = %3.3f \t K \n",T); +//// +////} diff --git a/src/Tests/eigenTest.cpp b/src/Tests/eigenTest.cpp index 76a3d89c..66803074 100644 --- a/src/Tests/eigenTest.cpp +++ b/src/Tests/eigenTest.cpp @@ -41,4 +41,31 @@ std::cout << CoolProp::vec_to_string(vec0) << std::endl; std::cout << CoolProp::vec_to_string(vec1) << std::endl; std::cout << CoolProp::vec_to_string(vec2) << std::endl; +Eigen::Matrix mat; +mat.setConstant(2,2,0.25); +std::vector< std::vector > vec; + +CoolProp::convert(mat, vec); +std::cout << CoolProp::vec_to_string(vec) << std::endl; + +//Eigen::Matrix mat; +//mat.resize(6,2); + +Eigen::Matrix mat2; +CoolProp::convert(vec2, mat2); +CoolProp::convert(mat2, vec); +std::cout << CoolProp::vec_to_string(vec) << std::endl; + +//std::vector< std::vector > vec(vec2); +//CoolProp::convert(mat,vec); + +//std::cout << CoolProp::vec_to_string() << std::endl; + +//Eigen::Matrix2d mat2 = CoolProp::convert(vec2); + +//Eigen::MatrixXd mat2(10,10); +//CoolProp::convert(vec2, mat2); + +//std::cout << CoolProp::vec_to_string(CoolProp::convert(mat2)) << std::endl; + }