diff --git a/dev/scripts/buildFiles.sh b/dev/scripts/buildFiles.sh new file mode 100644 index 00000000..388c438e --- /dev/null +++ b/dev/scripts/buildFiles.sh @@ -0,0 +1,16 @@ +#!/bin/bash +# +set -x +# Specify the files +SOURCES=( "CoolPropTools" "MatrixMath" "Solvers" "PolyMath" ) +ALLSRCS="" +# +LENGTH=${#SOURCES[@]} +# +# Loop through the sources and compile them +for (( i=0; i<${LENGTH}; i++ )); do + g++ -I include -c -o ${SOURCES[$i]}.o src/${SOURCES[$i]}.cpp + ALLSRCS="$ALLSRCS ${SOURCES[$i]}.o" +done +g++ $ALLSRCS -o test +exit 0 diff --git a/include/PolyMath.h b/include/PolyMath.h index 4fd1457a..94636e2e 100644 --- a/include/PolyMath.h +++ b/include/PolyMath.h @@ -25,12 +25,12 @@ public: // Destructor. No implementation virtual ~BasePolynomial(){}; -protected: +public: /// Basic checks for coefficient vectors. /** Starts with only the first coefficient dimension * and checks the vector length against parameter n. */ - 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); + 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 @@ -40,8 +40,8 @@ protected: * 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); + 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. @@ -49,8 +49,8 @@ protected: * 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); + 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 @@ -65,77 +65,77 @@ private: /// 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 T^0. */ - DEPRECATED(long double simplePolynomial(const std::vector &coefficients, long double T)); - DEPRECATED(long double simplePolynomial(const std::vector > &coefficients, long double x, long double T)); + * 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 T^0 */ - ///Indefinite integral in T-direction - long double simplePolynomialInt(const std::vector &coefficients, long double T); - ///Indefinite integral in T-direction only - long double simplePolynomialInt(const std::vector > &coefficients, long double x, long double T); + * 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 T^0 */ + * vector. Starts with only the first coefficient at x^0 */ ///Indefinite integral of a polynomial divided by its independent variable - long double simpleFracInt(const std::vector &coefficients, long double T); + double simpleFracInt(const std::vector &coefficients, double x); ///Indefinite integral of a polynomial divided by its 2nd independent variable - long double simpleFracInt(const std::vector > &coefficients, long double x, long double T); + 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 Tbase + * 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 T^0 */ + * Starts with only the first coefficient at x^0 */ ///Helper function to calculate the D vector: - long double factorial(long double nValue); - long double binom(long double nValue, long double nValue2); - std::vector fracIntCentralDvector(int m, long double T, long double Tbase); + 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 - long double fracIntCentral(const std::vector &coefficients, long double T, long double Tbase); + 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. */ - long double baseHorner(const std::vector &coefficients, long double T); - long double baseHorner(const std::vector< std::vector > &coefficients, long double x, long double T); - ///Indefinite integral in T-direction - long double baseHornerInt(const std::vector &coefficients, long double T); - ///Indefinite integral in T-direction only - long double baseHornerInt(const std::vector > &coefficients, long double x, long double T); + 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 - long double baseHornerFracInt(const std::vector &coefficients, long double T); + double baseHornerFracInt(const std::vector &coefficients, double x); ///Indefinite integral of a polynomial divided by its 2nd independent variable - long double baseHornerFracInt(const std::vector > &coefficients, long double x, long double T); + 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 T-direction - long double deriveIn2Steps(const std::vector &coefficients, long double T); - ///Derivative in terms of x(axis=true) or T(axis=false). - long double deriveIn2Steps(const std::vector< std::vector > &coefficients, long double x, long double T, bool axis); - ///Indefinite integral in T-direction - long double integrateIn2Steps(const std::vector &coefficients, long double T); - ///Indefinite integral in terms of x(axis=true) or T(axis=false). - long double integrateIn2Steps(const std::vector< std::vector > &coefficients, long double x, long double T, bool axis); - ///Indefinite integral in T-direction of a polynomial divided by its independent variable - long double fracIntIn2Steps(const std::vector &coefficients, long double T); - ///Indefinite integral in T-direction of a polynomial divided by its 2nd independent variable - long double fracIntIn2Steps(const std::vector > &coefficients, long double x, long double T); + ///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 - long double fracIntCentral2Steps(const std::vector > &coefficients, long double x, long double T, long double Tbase); + double fracIntCentral2Steps(const std::vector > &coefficients, double x, double y, double ybase); public: /** Here we define the functions that should be used by the @@ -146,491 +146,381 @@ public: */ /** Everything related to the normal polynomials goes in this - * section, holds functions for both evaluation and solving - * of polynomials. + * 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 long double value that represents the current input - inline long double polyval(const std::vector &coefficients, long double x){ + /// @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 long double value that represents the current input in the 1st dimension - /// @param y long double value that represents the current input in the 2nd dimension - inline long double polyval(const std::vector< std::vector > &coefficients, long double x, long double y){ + /// @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 functions for both evaluation and solving - * of polynomials. + * 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 T long double value that represents the current input - inline long double polyint(const std::vector &coefficients, long double T){ - return baseHornerInt(coefficients,T); + /// @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 (T) + /// Evaluates the indefinite integral of a two-dimensional polynomial along the 2nd axis (y) /// @param coefficients vector containing the ordered coefficients - /// @param x long double value that represents the current input in the 1st dimension - /// @param T long double value that represents the current input in the 2nd dimension - inline long double polyint(const std::vector< std::vector > &coefficients, long double x, long double T){ - return baseHornerInt(coefficients,x,T); + /// @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 functions for both evaluation and solving - * of polynomials. + * section, holds all the functions for evaluating polynomials. */ /// Evaluates the derivative of a one-dimensional polynomial /// @param coefficients vector containing the ordered coefficients - /// @param T long double value that represents the current input - inline long double polyder(const std::vector &coefficients, long double T){ - return deriveIn2Steps(coefficients,T); + /// @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 (T) + /// Evaluates the derivative of a two-dimensional polynomial along the 2nd axis (y) /// @param coefficients vector containing the ordered coefficients - /// @param x long double value that represents the current input in the 1st dimension - /// @param T long double value that represents the current input in the 2nd dimension - inline long double polyder(const std::vector< std::vector > &coefficients, long double x, long double T){ - return deriveIn2Steps(coefficients,x,T,false); + /// @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 functions for both evaluation and solving - * of polynomials. + * 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 T long double value that represents the current position - inline long double polyfracval(const std::vector &coefficients, long double T){ - return baseHornerFracInt(coefficients,T); + /// @param x double value that represents the current position + virtual inline double polyfracval(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 long double value that represents the current input in the 1st dimension - /// @param T long double value that represents the current input in the 2nd dimension - inline long double polyfracval(const std::vector< std::vector > &coefficients, long double x, long double T){ - return baseHornerFracInt(coefficients,x,T); + /// @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 baseHornerFracInt(coefficients,x,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 T long double value that represents the current position - inline long double polyfracint(const std::vector &coefficients, long double T){ - return baseHornerFracInt(coefficients,T); + /// @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 long double value that represents the current input in the 1st dimension - /// @param T long double value that represents the current input in the 2nd dimension - inline long double polyfracint(const std::vector< std::vector > &coefficients, long double x, long double T){ - return baseHornerFracInt(coefficients,x,T); + /// @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 T long double value that represents the current position - /// @param Tbase central temperature for fitted function - inline long double polyfracintcentral(const std::vector &coefficients, long double T, long double Tbase){ - return fracIntCentral(coefficients,T,Tbase); + /// @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 long double value that represents the current input in the 1st dimension - /// @param T long double value that represents the current input in the 2nd dimension - /// @param Tbase central temperature for fitted function - inline long double polyfracintcentral(const std::vector< std::vector > &coefficients, long double x, long double T, long double Tbase){ - return fracIntCentral2Steps(coefficients,x,T,Tbase); + /// @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(){}; + virtual double call(double x); + virtual double deriv(double x); +}; +class PolyIntResidual : public PolyResidual { +public: + virtual double call(double x); + virtual double deriv(double x); +}; +class PolyFracIntResidual : public PolyResidual { +public: + virtual double call(double x); + virtual double deriv(double x); +}; +class PolyDerResidual : public PolyResidual { +public: + 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 DEBUG; + +public: + BaseExponential(); + virtual ~BaseExponential(){}; + +public: /// Evaluates an exponential function for the given coefficients /// @param coefficients vector containing the ordered coefficients - /// @param T long double value that represents the current input + /// @param x double value that represents the current input /// @param n int value that determines the kind of exponential function - long double expval(const std::vector &coefficients, long double T, int n); + 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 long double value that represents the current input in the 1st dimension - /// @param T long double value that represents the current input in the 2nd dimension + /// @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 - long double expval(const std::vector< std::vector > &coefficients, long double x, long double T, int n); -}; - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -/// The classes for Polynomials -class PolynomialImpl1D : public BasePolynomial{ -protected: - std::vector coefficients; - /// A nested class that is used by the solvers to calculate - /// residuals and derivatives during the solution process. - class Residual : public FuncWrapper1D { - private: - PolynomialImpl1D *poly; - long double y; - Residual(); - public: - Residual(PolynomialImpl1D *poly, long double y); - virtual double call(double x); - virtual double deriv(double x); - }; - class ResidualInt : public Residual { - public: - virtual double call(double x); - virtual double deriv(double x); - }; - class ResidualDer : public Residual { - public: - virtual double call(double x); - virtual double deriv(double x); - }; -private: - PolynomialImpl1D(); -public: - PolynomialImpl1D(const std::vector &coefficients); - virtual ~PolynomialImpl1D(){}; - /// Evaluates a one-dimensional polynomial for the given coefficients - /// @param x long double value that represents the current input - virtual long double eval(long double x); - /// Evaluates the indefinite integral of a one-dimensional polynomial - /// @param x long double value that represents the current input - virtual long double integ(long double x); - /// Evaluates the derivative of a one-dimensional polynomial - /// @param x long double value that represents the current input - virtual long double deriv(long double x); - /// Solves a one-dimensional polynomial for the given coefficients - /// @param y long double value that represents the current function output - /// @param x0 long double value that represents the first guess for x - virtual long double solve(long double y, long double x0); - /// Solves a one-dimensional polynomial for the given coefficients - /// @param y long double value that represents the current function output - /// @param xmin long double value that represents the lower limit for x - /// @param xmax long double value that represents the upper limit for x - virtual long double solve(long double y, long double xmin, long double xmax); - /// Solves an integrated one-dimensional polynomial for the given coefficients - /// @param y long double value that represents the current function output - /// @param x0 long double value that represents the first guess for x - virtual long double solveInt(long double y, long double x0); - /// Solves an integrated one-dimensional polynomial for the given coefficients - /// @param y long double value that represents the current function output - /// @param xmin long double value that represents the lower limit for x - /// @param xmax long double value that represents the upper limit for x - virtual long double solveInt(long double y, long double xmin, long double xmax); - /// Solves the derivative of a one-dimensional polynomial for the given coefficients - /// @param y long double value that represents the current function output - /// @param x0 long double value that represents the first guess for x - virtual long double solveDer(long double y, long double x0); - /// Solves the derivative of a one-dimensional polynomial for the given coefficients - /// @param y long double value that represents the current function output - /// @param xmin long double value that represents the lower limit for x - /// @param xmax long double value that represents the upper limit for x - virtual long double solveDer(long double y, long double xmin, long double xmax); -}; - -class PolynomialImpl2D : public BasePolynomial{ -protected: - std::vector< std::vector > coefficients; - /// A nested class that is used by the solvers to calculate - /// residuals and derivatives during the solution process. - class Residual : public FuncWrapper1D { - private: - PolynomialImpl2D *poly; - long double x, y; - Residual(); - public: - Residual(PolynomialImpl2D *poly, long double y, long double x); - virtual double call(double z); - virtual double deriv(double z); - }; - class ResidualInt : public Residual { - public: - virtual double call(double z); - virtual double deriv(double z); - }; - class ResidualDer : public Residual { - public: - virtual double call(double z); - virtual double deriv(double z); - }; -private: - PolynomialImpl2D(); -public: - PolynomialImpl2D(const std::vector< std::vector > &coefficients); - virtual ~PolynomialImpl2D(){}; - /// Evaluates a two-dimensional polynomial for the given coefficients - /// @param x long double value that represents the current input in the 1st dimension - /// @param z long double value that represents the current input in the 2nd dimension - virtual long double eval(long double x, long double z); - /// Evaluates the indefinite integral of a two-dimensional polynomial along the 2nd axis (z) - /// @param x long double value that represents the current input in the 1st dimension - /// @param z long double value that represents the current input in the 2nd dimension - virtual long double integ(long double x, long double z); - /// Evaluates the derivative of a two-dimensional polynomial along the 2nd axis (z) - /// @param coefficients vector containing the ordered coefficients - /// @param x long double value that represents the current input in the 1st dimension - /// @param z long double value that represents the current input in the 2nd dimension - virtual long double deriv(long double x, long double z); - /// Solves a two-dimensional polynomial for the 2nd input (z) - /// @param y long double value that represents the current function output - /// @param x long double value that represents the current input in the 1st dimension - /// @param z0 long double value that represents the first guess for z - virtual long double solve(long double y, long double x, long double z0); - /// Solves a two-dimensional polynomial for the 2nd input (z) - /// @param y long double value that represents the current function output - /// @param x long double value that represents the current input in the 1st dimension - /// @param zmin long double value that represents the lower limit for z - /// @param zmax long double value that represents the upper limit for z - virtual long double solve(long double y, long double x, long double zmin, long double zmax); - /// Solves an integrated two-dimensional polynomial for the 2nd input (z) - /// @param y long double value that represents the current function output - /// @param x long double value that represents the current input in the 1st dimension - /// @param z0 long double value that represents the first guess for z - virtual long double solveInt(long double y, long double x, long double z0); - /// Solves an integrated two-dimensional polynomial for the 2nd input (z) - /// @param y long double value that represents the current function output - /// @param x long double value that represents the current input in the 1st dimension - /// @param zmin long double value that represents the lower limit for z - /// @param zmax long double value that represents the upper limit for z - virtual long double solveInt(long double y, long double x, long double zmin, long double zmax); - /// Solves the derivative of a two-dimensional polynomial for the 2nd input (z) - /// @param y long double value that represents the current function output - /// @param x long double value that represents the current input in the 1st dimension - /// @param z0 long double value that represents the first guess for z - virtual long double solveDer(long double y, long double x, long double z0); - /// Solves the derivative of a two-dimensional polynomial for the 2nd input (z) - /// @param y long double value that represents the current function output - /// @param x long double value that represents the current input in the 1st dimension - /// @param zmin long double value that represents the lower limit for z - /// @param zmax long double value that represents the upper limit for z - virtual long double solveDer(long double y, long double x, long double zmin, long double zmax); -}; - -class PolynomialFrac1D : public PolynomialImpl1D{ -private: - PolynomialFrac1D(); -public: - PolynomialFrac1D(const std::vector &coefficients); - virtual ~PolynomialFrac1D(){}; - /// Evaluates a one-dimensional polynomial for the given coefficients - /// @param x long double value that represents the current input - virtual long double eval(long double x); - /// Evaluates the indefinite integral of a one-dimensional polynomial - /// @param x long double value that represents the current input - virtual long double integ(long double x); - /// Evaluates the derivative of a one-dimensional polynomial - /// @param x long double value that represents the current input - virtual long double deriv(long double x); - /// Solves a one-dimensional polynomial for the given coefficients - /// @param y long double value that represents the current function output - /// @param x0 long double value that represents the first guess for x - virtual long double solve(long double y, long double x0); - /// Solves a one-dimensional polynomial for the given coefficients - /// @param y long double value that represents the current function output - /// @param xmin long double value that represents the lower limit for x - /// @param xmax long double value that represents the upper limit for x - virtual long double solve(long double y, long double xmin, long double xmax); - /// Solves an integrated one-dimensional polynomial for the given coefficients - /// @param y long double value that represents the current function output - /// @param x0 long double value that represents the first guess for x - virtual long double solveInt(long double y, long double x0); - /// Solves an integrated one-dimensional polynomial for the given coefficients - /// @param y long double value that represents the current function output - /// @param xmin long double value that represents the lower limit for x - /// @param xmax long double value that represents the upper limit for x - virtual long double solveInt(long double y, long double xmin, long double xmax); - // virtual long double solveIntCentral(long double y, long double xBase, long double x0); // TODO: implement solveIntCentral with x0 - /// Solves an integrated one-dimensional polynomial - /// @param y long double value that represents the current function output - /// @param xBase long double value that represents the central value for x - /// @param xmin long double value that represents the lower limit for x - /// @param xmax long double value that represents the upper limit for x - virtual long double solveIntCentral(long double y, long double xBase, long double xmin, long double xmax); - /// Solves the derivative of a one-dimensional polynomial for the given coefficients - /// @param y long double value that represents the current function output - /// @param x0 long double value that represents the first guess for x - virtual long double solveDer(long double y, long double x0); - /// Solves the derivative of a one-dimensional polynomial for the given coefficients - /// @param y long double value that represents the current function output - /// @param xmin long double value that represents the lower limit for x - /// @param xmax long double value that represents the upper limit for x - virtual long double solveDer(long double y, long double xmin, long double xmax); - // virtual long double solveDerCentral(long double y, long double xBase, long double x0); // TODO: implement solveDerCentral with x0 - // virtual long double solveDerCentral(long double y, long double xBase, long double xmin, long double xmax); // TODO: implement solveDerCentral with xmin, xmax - - - - - - - - - /// Solves a one-dimensional polynomial for the given coefficients - /// @param coefficients vector containing the ordered coefficients - /// @param y long double value that represents the current function output - /// @param x0 long double value that represents the first guess for x - virtual long double solve(const std::vector &coefficients, long double y, long double x0); - /// Solves a two-dimensional polynomial for the 2nd input (z) - /// @param coefficients vector containing the ordered coefficients - /// @param y long double value that represents the current function output - /// @param x long double value that represents the current input in the 1st dimension - /// @param z0 long double value that represents the first guess for z - virtual long double solve(const std::vector< std::vector > &coefficients, long double y, long double x, long double z0); - /// Solves a one-dimensional polynomial for the given coefficients - /// @param coefficients vector containing the ordered coefficients - /// @param y long double value that represents the current function output - /// @param xmin long double value that represents the lower limit for x - /// @param xmax long double value that represents the upper limit for x - virtual long double solve(const std::vector &coefficients, long double y, long double xmin, long double xmax); - /// Solves a two-dimensional polynomial for the 2nd input (z) - /// @param coefficients vector containing the ordered coefficients - /// @param y long double value that represents the current function output - /// @param x long double value that represents the current input in the 1st dimension - /// @param zmin long double value that represents the lower limit for z - /// @param zmax long double value that represents the upper limit for z - virtual long double solve(const std::vector< std::vector > &coefficients, long double y, long double x, long double zmin, long double zmax); - /// Solves an integrated one-dimensional polynomial for the given coefficients - /// @param coefficients vector containing the ordered coefficients - /// @param y long double value that represents the current function output - /// @param x0 long double value that represents the first guess for x - virtual long double solveInt(const std::vector &coefficients, long double y, long double x0); - /// Solves an integrated two-dimensional polynomial for the 2nd input (z) - /// @param coefficients vector containing the ordered coefficients - /// @param y long double value that represents the current function output - /// @param x long double value that represents the current input in the 1st dimension - /// @param z0 long double value that represents the first guess for z - virtual long double solveInt(const std::vector< std::vector > &coefficients, long double y, long double x, long double z0); - /// Solves an integrated one-dimensional polynomial for the given coefficients - /// @param coefficients vector containing the ordered coefficients - /// @param y long double value that represents the current function output - /// @param xmin long double value that represents the lower limit for x - /// @param xmax long double value that represents the upper limit for x - virtual long double solveInt(const std::vector &coefficients, long double y, long double xmin, long double xmax); - /// Solves an integrated two-dimensional polynomial for the 2nd input (z) - /// @param coefficients vector containing the ordered coefficients - /// @param y long double value that represents the current function output - /// @param x long double value that represents the current input in the 1st dimension - /// @param zmin long double value that represents the lower limit for z - /// @param zmax long double value that represents the upper limit for z - virtual long double solveInt(const std::vector< std::vector > &coefficients, long double y, long double x, long double zmin, long double zmax); - /// Solves the derivative of a one-dimensional polynomial for the given coefficients - /// @param coefficients vector containing the ordered coefficients - /// @param y long double value that represents the current function output - /// @param x0 long double value that represents the first guess for x - virtual long double solveDer(const std::vector &coefficients, long double y, long double x0); - /// Solves the derivative of a one-dimensional polynomial for the given coefficients - /// @param coefficients vector containing the ordered coefficients - /// @param y long double value that represents the current function output - /// @param xmin long double value that represents the lower limit for x - /// @param xmax long double value that represents the upper limit for x - virtual long double solveDer(const std::vector &coefficients, long double y, long double xmin, long double xmax); - - - - - - - - - - - /// Solves the derivative of a two-dimensional polynomial for the 2nd input (z) - /// @param coefficients vector containing the ordered coefficients - /// @param y long double value that represents the current function output - /// @param x long double value that represents the current input in the 1st dimension - /// @param z0 long double value that represents the first guess for z - virtual long double solveDer(const std::vector< std::vector > &coefficients, long double y, long double x, long double z0); - /// Solves the derivative of a two-dimensional polynomial for the 2nd input (z) - /// @param coefficients vector containing the ordered coefficients - /// @param y long double value that represents the current function output - /// @param x long double value that represents the current input in the 1st dimension - /// @param zmin long double value that represents the lower limit for z - /// @param zmax long double value that represents the upper limit for z - virtual long double solveDer(const std::vector< std::vector > &coefficients, long double y, long double x, long double zmin, long double zmax); - // virtual long double solveDerCentral(long double y, long double x, long double zBase, long double z0); // TODO: implement solveDerCentral with z0 - // virtual long double solveDerCentral(long double y, long double x, long double zBase, long double zmin, long double zmax); // TODO: implement solveDerCentral with zmin, zmax - - - - - - - - - - - -}; - -class PolynomialFrac2D : public PolynomialImpl2D{ -private: - PolynomialFrac2D(); -public: - PolynomialFrac2D(const std::vector< std::vector > &coefficients); - virtual ~PolynomialFrac2D(){}; - /// Evaluates a two-dimensional polynomial for the given coefficients - /// @param x long double value that represents the current input in the 1st dimension - /// @param z long double value that represents the current input in the 2nd dimension - virtual long double eval(long double x, long double z); - /// Evaluates the indefinite integral of a two-dimensional polynomial along the 2nd axis (z) - /// @param x long double value that represents the current input in the 1st dimension - /// @param z long double value that represents the current input in the 2nd dimension - virtual long double integ(long double x, long double z); - /// Evaluates the indefinite integral of a two-dimensional polynomial along the 2nd axis (z) - /// @param x long double value that represents the current input in the 1st dimension - /// @param z long double value that represents the current input in the 2nd dimension - /// @param zBase long double value that represents the central value for z - virtual long double integCentral(long double x, long double z, long double zBase); - /// Evaluates the derivative of a two-dimensional polynomial along the 2nd axis (z) - /// @param coefficients vector containing the ordered coefficients - /// @param x long double value that represents the current input in the 1st dimension - /// @param z long double value that represents the current input in the 2nd dimension - virtual long double deriv(long double x, long double z); - - + double expval(const std::vector< std::vector > &coefficients, double x, double y, int n); }; diff --git a/src/PolyMath.cpp b/src/PolyMath.cpp index 5ab75e7c..675f9d20 100644 --- a/src/PolyMath.cpp +++ b/src/PolyMath.cpp @@ -3,6 +3,7 @@ #include "CoolPropTools.h" #include "Exceptions.h" +#include "MatrixMath.h" #include #include @@ -18,10 +19,11 @@ BasePolynomial::BasePolynomial(){ this->DEBUG = false; } + /// 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){ +bool BasePolynomial::checkCoefficients(const std::vector &coefficients, unsigned int n){ if (coefficients.size() == n){ return true; } else { @@ -29,8 +31,7 @@ bool BasePolynomial::checkCoefficients(const std::vector &coefficie } return false; } - -bool BasePolynomial::checkCoefficients(std::vector< std::vector > const& coefficients, unsigned int rows, unsigned int columns){ +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 > c return false; } + /** Integrating coefficients for polynomials is done by dividing the * original coefficients by (i+1) and elevating the order by 1. * Some reslicing needs to be applied to integrate along the x-axis. */ -std::vector BasePolynomial::integrateCoeffs(std::vector const& coefficients){ - std::vector newCoefficients; +std::vector 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 @@ -58,15 +60,13 @@ std::vector BasePolynomial::integrateCoeffs(std::vector > BasePolynomial::integrateCoeffs(std::vector< std::vector > const& coefficients, bool axis){ - std::vector< std::vector > newCoefficients; +std::vector< std::vector > 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; + std::vector< std::vector > tmpCoefficients; tmpCoefficients = transpose(coefficients); unsigned int sizeY = tmpCoefficients.size(); for(unsigned int i=0; i > BasePolynomial::integrateCoeffs(std::vec return newCoefficients; } + /** Deriving coefficients for polynomials is done by multiplying the * original coefficients with i and lowering the order by 1. */ -std::vector BasePolynomial::deriveCoeffs(std::vector const& coefficients){ - std::vector newCoefficients; +std::vector 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 @@ -97,14 +98,13 @@ std::vector BasePolynomial::deriveCoeffs(std::vector c } return newCoefficients; } - -std::vector< std::vector > BasePolynomial::deriveCoeffs(const std::vector< std::vector > &coefficients, unsigned int axis){ - std::vector< std::vector > newCoefficients; +std::vector< std::vector > 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; + std::vector< std::vector > tmpCoefficients; tmpCoefficients = transpose(coefficients); unsigned int sizeY = tmpCoefficients.size(); for(unsigned int i=0; i > BasePolynomial::deriveCoeffs(const std:: return newCoefficients; } + /** 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 @@ -134,16 +135,16 @@ std::vector< std::vector > BasePolynomial::deriveCoeffs(const std:: /// 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 T^0. */ -long double BasePolynomial::simplePolynomial(std::vector const& coefficients, long double T){ + * Starts with only the first coefficient at x^0. */ +double BasePolynomial::simplePolynomial(std::vector const& coefficients, double x){ if (this->DEBUG) { - std::cout << "Running simplePolynomial(std::vector, " << T << "): "; + std::cout << "Running simplePolynomial(std::vector, " << x << "): "; } bool db = this->DEBUG; this->DEBUG = false; - long double result = 0.; + double result = 0.; for(unsigned int i=0; iDEBUG = db; if (this->DEBUG) { @@ -151,15 +152,15 @@ long double BasePolynomial::simplePolynomial(std::vector const& coe } return result; } -long double BasePolynomial::simplePolynomial(std::vector > const& coefficients, long double x, long double T){ +double BasePolynomial::simplePolynomial(std::vector > const& coefficients, double x, double y){ if (this->DEBUG) { - std::cout << "Running simplePolynomial(std::vector, " << x << ", " << T << "): "; + std::cout << "Running simplePolynomial(std::vector, " << x << ", " << y << "): "; } bool db = this->DEBUG; this->DEBUG = false; - long double result = 0; + double result = 0; for(unsigned int i=0; iDEBUG = db; if (this->DEBUG) { @@ -173,17 +174,17 @@ long double BasePolynomial::simplePolynomial(std::vector const& coefficients, long double T){ + * 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->DEBUG) { - std::cout << "Running simplePolynomialInt(std::vector, " << T << "): "; + std::cout << "Running simplePolynomialInt(std::vector, " << x << "): "; } bool db = this->DEBUG; this->DEBUG = false; - long double result = 0.; + double result = 0.; for(unsigned int i=0; iDEBUG = db; if (this->DEBUG) { @@ -193,15 +194,15 @@ long double BasePolynomial::simplePolynomialInt(std::vector const& } ///Indefinite integral in y-direction only -long double BasePolynomial::simplePolynomialInt(std::vector > const& coefficients, long double x, long double T){ +double BasePolynomial::simplePolynomialInt(std::vector > const& coefficients, double x, double y){ if (this->DEBUG) { - std::cout << "Running simplePolynomialInt(std::vector, " << x << ", " << T << "): "; + std::cout << "Running simplePolynomialInt(std::vector, " << x << ", " << y << "): "; } bool db = this->DEBUG; this->DEBUG = false; - long double result = 0.; + double result = 0.; for(unsigned int i=0; iDEBUG = db; if (this->DEBUG) { @@ -215,15 +216,15 @@ long double BasePolynomial::simplePolynomialInt(std::vector const& coefficients, long double T){ + * Starts with only the first coefficient at x^0 */ +double BasePolynomial::simpleFracInt(std::vector const& coefficients, double x){ if (this->DEBUG) { - std::cout << "Running simpleFracInt(std::vector, " << T << "): "; + std::cout << "Running simpleFracInt(std::vector, " << x << "): "; } - long double result = coefficients[0] * log(T); + double result = coefficients[0] * log(x); if (coefficients.size() > 1) { for (unsigned int i=1; iDEBUG) { @@ -232,15 +233,15 @@ long double BasePolynomial::simpleFracInt(std::vector const& coeffi return result; } -long double BasePolynomial::simpleFracInt(std::vector< std::vector > const& coefficients, long double x, long double T){ +double BasePolynomial::simpleFracInt(std::vector< std::vector > const& coefficients, double x, double y){ if (this->DEBUG) { - std::cout << "Running simpleFracInt(std::vector, " << x << ", " << T << "): "; + std::cout << "Running simpleFracInt(std::vector, " << x << ", " << y << "): "; } bool db = this->DEBUG; this->DEBUG = false; - long double result = 0; + double result = 0; for (unsigned int i=0; iDEBUG = db; if (this->DEBUG) { @@ -252,17 +253,17 @@ long double BasePolynomial::simpleFracInt(std::vector< std::vector /** 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 Tbase + * 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 T^0 */ + * 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 -//long double BasePolynomial::factorial(long double nValue){ -// long double result = nValue; -// long double result_next; -// long double pc = nValue; +//double BasePolynomial::factorial(double nValue){ +// double result = nValue; +// double result_next; +// double pc = nValue; // do { // result_next = result*(pc-1); // result = result_next; @@ -271,20 +272,20 @@ long double BasePolynomial::simpleFracInt(std::vector< std::vector // nValue = result; // return nValue; //} -//long double BasePolynomial::factorial(long double nValue){ +//double BasePolynomial::factorial(double nValue){ // if (nValue == 0) return (1); // else return (nValue * factorial(nValue - 1)); //} -long double BasePolynomial::factorial(long double nValue){ - long double value = 1; +double BasePolynomial::factorial(double nValue){ + double value = 1; for(int i = 2; i <= nValue; i++){ value = value * i; } return value; } -long double BasePolynomial::binom(long double nValue, long double nValue2){ - long double result; +double BasePolynomial::binom(double nValue, double nValue2){ + double result; if(nValue2 == 1) return nValue; result = (factorial(nValue)) / (factorial(nValue2)*factorial((nValue - nValue2))); nValue2 = result; @@ -292,14 +293,14 @@ long double BasePolynomial::binom(long double nValue, long double nValue2){ } ///Helper functions to calculate the D vector: -std::vector BasePolynomial::fracIntCentralDvector(int m, long double T, long double Tbase){ - std::vector D; - long double tmp; +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 BasePolynomial::fracIntCentralDvector(int m, long doubl } ///Indefinite integral of a centred polynomial divided by its independent variable -long double BasePolynomial::fracIntCentral(std::vector const& coefficients, long double T, long double Tbase){ +double BasePolynomial::fracIntCentral(std::vector const& coefficients, double x, double xbase){ if (this->DEBUG) { - std::cout << "Running fracIntCentral(std::vector, " << T << ", " << Tbase << "): "; + std::cout << "Running fracIntCentral(std::vector, " << x << ", " << xbase << "): "; } bool db = this->DEBUG; this->DEBUG = false; int m = coefficients.size(); - std::vector D = fracIntCentralDvector(m, T, Tbase); - long double result = 0; + std::vector D = fracIntCentralDvector(m, x, xbase); + double result = 0; for(int j=0; j const& coeff * This avoids unnecessary multiplication and thus * speeds up calculation. */ -long double BasePolynomial::baseHorner(std::vector const& coefficients, long double T){ +double BasePolynomial::baseHorner(std::vector const& coefficients, double x){ if (this->DEBUG) { - std::cout << "Running baseHorner(std::vector, " << T << "): "; + std::cout << "Running baseHorner(std::vector, " << x << "): "; } - long double result = 0; + double result = 0; for(int i=coefficients.size()-1; i>=0; i--) { - result *= T; + result *= x; result += coefficients[i]; } if (this->DEBUG) { @@ -347,16 +348,16 @@ long double BasePolynomial::baseHorner(std::vector const& coefficie return result; } -long double BasePolynomial::baseHorner(std::vector< std::vector > const& coefficients, long double x, long double T){ +double BasePolynomial::baseHorner(std::vector< std::vector > const& coefficients, double x, double y){ if (this->DEBUG) { - std::cout << "Running baseHorner(std::vector, " << x << ", " << T << "): "; + std::cout << "Running baseHorner(std::vector, " << x << ", " << y << "): "; } bool db = this->DEBUG; this->DEBUG = false; - long double result = 0; + double result = 0; for(int i=coefficients.size()-1; i>=0; i--) { result *= x; - result += baseHorner(coefficients[i], T); + result += baseHorner(coefficients[i], y); } this->DEBUG = db; if (this->DEBUG) { @@ -365,34 +366,34 @@ long double BasePolynomial::baseHorner(std::vector< std::vector > c return result; } -///Indefinite integral in T-direction -long double BasePolynomial::baseHornerInt(std::vector const& coefficients, long double T){ +///Indefinite integral in x-direction +double BasePolynomial::baseHornerInt(std::vector const& coefficients, double x){ if (this->DEBUG) { - std::cout << "Running baseHornerInt(std::vector, " << T << "): "; + std::cout << "Running baseHornerInt(std::vector, " << x << "): "; } - long double result = 0; + double result = 0; for(int i=coefficients.size()-1; i>=0; i--) { - result *= T; + result *= x; result += coefficients[i]/(i+1.); } - result = result * T; + result = result * x; if (this->DEBUG) { std::cout << result << std::endl; } return result; } -///Indefinite integral in T-direction only -long double BasePolynomial::baseHornerInt(std::vector > const& coefficients, long double x, long double T){ +///Indefinite integral in y-direction only +double BasePolynomial::baseHornerInt(std::vector > const& coefficients, double x, double y){ if (this->DEBUG) { - std::cout << "Running baseHornerInt(std::vector, " << x << ", " << T << "): "; + std::cout << "Running baseHornerInt(std::vector, " << x << ", " << y << "): "; } bool db = this->DEBUG; this->DEBUG = false; - long double result = 0; + double result = 0; for(int i=coefficients.size()-1; i>=0; i--) { result *= x; - result += baseHornerInt(coefficients[i], T); + result += baseHornerInt(coefficients[i], y); } this->DEBUG = db; if (this->DEBUG) { @@ -401,22 +402,22 @@ long double BasePolynomial::baseHornerInt(std::vector > return result; } -///Indefinite integral in T-direction of a polynomial divided by its independent variable -long double BasePolynomial::baseHornerFracInt(std::vector const& coefficients, long double T){ +///Indefinite integral in x-direction of a polynomial divided by its independent variable +double BasePolynomial::baseHornerFracInt(std::vector const& coefficients, double x){ if (this->DEBUG) { - std::cout << "Running baseHornerFra(std::vector, " << T << "): "; + std::cout << "Running baseHornerFra(std::vector, " << x << "): "; } bool db = this->DEBUG; this->DEBUG = false; - long double result = 0; + double result = 0; if (coefficients.size() > 1) { for(int i=coefficients.size()-1; i>=1; i--) { - result *= T; + result *= x; result += coefficients[i]/(i); } - result *= T; + result *= x; } - result += coefficients[0] * log(T); + result += coefficients[0] * log(x); this->DEBUG = db; if (this->DEBUG) { std::cout << result << std::endl; @@ -424,18 +425,18 @@ long double BasePolynomial::baseHornerFracInt(std::vector const& co return result; } -///Indefinite integral in T-direction of a polynomial divided by its 2nd independent variable -long double BasePolynomial::baseHornerFracInt(std::vector > const& coefficients, long double x, long double T){ +///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->DEBUG) { - std::cout << "Running baseHornerFra(std::vector, " << x << ", " << T << "): "; + std::cout << "Running baseHornerFra(std::vector, " << x << ", " << y << "): "; } bool db = this->DEBUG; this->DEBUG = false; - long double result = 0; + double result = 0; for(int i=coefficients.size()-1; i>=0; i--) { result *= x; - result += baseHornerFracInt(coefficients[i], T); + result += baseHornerFracInt(coefficients[i], y); } this->DEBUG = db; @@ -451,14 +452,14 @@ long double BasePolynomial::baseHornerFracInt(std::vector const& coefficients, long double T){ +///Derivative in x-direction +double BasePolynomial::deriveIn2Steps(std::vector const& coefficients, double x){ if (this->DEBUG) { - std::cout << "Running deriveIn2Steps(std::vector, " << T << "): "; + std::cout << "Running deriveIn2Steps(std::vector, " << x << "): "; } bool db = this->DEBUG; this->DEBUG = false; - long double result = polyval(deriveCoeffs(coefficients),T); + double result = polyval(deriveCoeffs(coefficients),x); this->DEBUG = db; if (this->DEBUG) { std::cout << result << std::endl; @@ -466,14 +467,14 @@ long double BasePolynomial::deriveIn2Steps(std::vector const& coeff return result; } -///Derivative in terms of x(axis=true) or T(axis=false). -long double BasePolynomial::deriveIn2Steps(std::vector< std::vector > const& coefficients, long double x, long double T, bool axis){ +///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->DEBUG) { - std::cout << "Running deriveIn2Steps(std::vector, " << x << ", " << T << "): "; + std::cout << "Running deriveIn2Steps(std::vector, " << x << ", " << y << "): "; } bool db = this->DEBUG; this->DEBUG = false; - long double result = polyval(deriveCoeffs(coefficients,axis),x,T); + double result = polyval(deriveCoeffs(coefficients,axis),x,y); this->DEBUG = db; if (this->DEBUG) { std::cout << result << std::endl; @@ -481,14 +482,14 @@ long double BasePolynomial::deriveIn2Steps(std::vector< std::vector return result; } -///Indefinite integral in T-direction -long double BasePolynomial::integrateIn2Steps(std::vector const& coefficients, long double T){ +///Indefinite integral in x-direction +double BasePolynomial::integrateIn2Steps(std::vector const& coefficients, double x){ if (this->DEBUG) { - std::cout << "Running integrateIn2Steps(std::vector, " << T << "): "; + std::cout << "Running integrateIn2Steps(std::vector, " << x << "): "; } bool db = this->DEBUG; this->DEBUG = false; - long double result = polyval(integrateCoeffs(coefficients),T); + double result = polyval(integrateCoeffs(coefficients),x); this->DEBUG = db; if (this->DEBUG) { std::cout << result << std::endl; @@ -496,14 +497,14 @@ long double BasePolynomial::integrateIn2Steps(std::vector const& co return result; } -///Indefinite integral in terms of x(axis=true) or T(axis=false). -long double BasePolynomial::integrateIn2Steps(std::vector< std::vector > const& coefficients, long double x, long double T, bool axis){ +///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->DEBUG) { - std::cout << "Running integrateIn2Steps(std::vector, " << x << ", " << T << "): "; + std::cout << "Running integrateIn2Steps(std::vector, " << x << ", " << y << "): "; } bool db = this->DEBUG; this->DEBUG = false; - long double result = polyval(integrateCoeffs(coefficients,axis),x,T); + double result = polyval(integrateCoeffs(coefficients,axis),x,y); this->DEBUG = db; if (this->DEBUG) { std::cout << result << std::endl; @@ -511,17 +512,17 @@ long double BasePolynomial::integrateIn2Steps(std::vector< std::vector const& coefficients, long double T){ +///Indefinite integral in x-direction of a polynomial divided by its independent variable +double BasePolynomial::fracIntIn2Steps(std::vector const& coefficients, double x){ if (this->DEBUG) { - std::cout << "Running fracIntIn2Steps(std::vector, " << T << "): "; + std::cout << "Running fracIntIn2Steps(std::vector, " << x << "): "; } bool db = this->DEBUG; this->DEBUG = false; - long double result = coefficients[0] * log(T); + double result = coefficients[0] * log(x); if (coefficients.size() > 1) { - std::vector newCoeffs(coefficients.begin() + 1, coefficients.end()); - result += polyint(newCoeffs,T); + std::vector newCoeffs(coefficients.begin() + 1, coefficients.end()); + result += polyint(newCoeffs,x); } this->DEBUG = db; if (this->DEBUG) { @@ -530,18 +531,18 @@ long double BasePolynomial::fracIntIn2Steps(std::vector const& coef return result; } -///Indefinite integral in T-direction of a polynomial divided by its 2nd independent variable -long double BasePolynomial::fracIntIn2Steps(std::vector > const& coefficients, long double x, long double T){ +///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->DEBUG) { - std::cout << "Running fracIntIn2Steps(std::vector, " << x << ", " << T << "): "; + std::cout << "Running fracIntIn2Steps(std::vector, " << x << ", " << y << "): "; } bool db = this->DEBUG; this->DEBUG = false; - std::vector newCoeffs; + std::vector newCoeffs; for (unsigned int i=0; iDEBUG = db; if (this->DEBUG) { std::cout << result << std::endl; @@ -549,18 +550,18 @@ long double BasePolynomial::fracIntIn2Steps(std::vector return result; } -///Indefinite integral in T-direction of a centred polynomial divided by its 2nd independent variable -long double BasePolynomial::fracIntCentral2Steps(std::vector > const& coefficients, long double x, long double T, long double Tbase){ +///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->DEBUG) { - std::cout << "Running fracIntCentral2Steps(std::vector, " << x << ", " << T << ", " << Tbase << "): "; + std::cout << "Running fracIntCentral2Steps(std::vector, " << x << ", " << y << ", " << ybase << "): "; } bool db = this->DEBUG; this->DEBUG = false; - std::vector newCoeffs; + std::vector newCoeffs; for (unsigned int i=0; iDEBUG = db; if (this->DEBUG) { std::cout << result << std::endl; @@ -569,22 +570,335 @@ long double BasePolynomial::fracIntCentral2Steps(std::vectordim = -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){ + //throw CoolProp::NotImplementedError("Please redefine your classes locally."); + double polyRes = -1; + if (this->dim==i1D) { + polyRes = this->poly.polyval(this->coefficients[0], x); + } else if (this->dim==i2D) { + polyRes = this->poly.polyval(this->coefficients, this->firstDim, x); + } else { + 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){ +// throw CoolProp::NotImplementedError("Please redefine your classes locally."); + double polyRes = -1; + if (this->dim==i1D) { + polyRes = this->poly.polyder(this->coefficients[0], x); + } else if (this->dim==i2D) { + polyRes = this->poly.polyder(this->coefficients, this->firstDim, x); + } else { + throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); + } + return polyRes; +} + +double PolyIntResidual::call(double x){ + //throw CoolProp::NotImplementedError("Please redefine your classes locally."); + double polyRes = -1; + if (this->dim==i1D) { + polyRes = this->poly.polyint(this->coefficients[0], x); + } else if (this->dim==i2D) { + polyRes = this->poly.polyint(this->coefficients, this->firstDim, x); + } else { + 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){ +// throw CoolProp::NotImplementedError("Please redefine your classes locally."); + double polyRes = -1; + if (this->dim==i1D) { + polyRes = this->poly.polyval(this->coefficients[0], x); + } else if (this->dim==i2D) { + polyRes = this->poly.polyval(this->coefficients, this->firstDim, x); + } else { + throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D."); + } + return polyRes; +} + +double PolyDerResidual::call(double x){ + //throw CoolProp::NotImplementedError("Please redefine your classes locally."); + double polyRes = -1; + if (this->dim==i1D) { + polyRes = this->poly.polyder(this->coefficients[0], x); + } else if (this->dim==i2D) { + polyRes = this->poly.polyder(this->coefficients, this->firstDim, x); + } else { + 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 + * 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->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){ + throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function +} + +/// 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){ + throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function +} + + +/** 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){ + throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function +} + +/// 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){ + throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function +} + + +/** 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){ + 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::polyfracint(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 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){ + throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function +} + +/// 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){ + throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function +} + + +/** 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 + 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."); + } + return result; +} + + +/** Here we define the functions that should be to evaluate exponential + * functions. Not really polynomials, I know... + */ + +BaseExponential::BaseExponential(){ + this->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 T long double value that represents the current input +/// @param x double value that represents the current input /// @param n int value that determines the kind of exponential function -long double BasePolynomial::expval(std::vector const& coefficients, long double T, int n){ - long double result = 0.; +double BaseExponential::expval(const std::vector &coefficients, double x, int n){ + double result = 0.; if (n==1) { - checkCoefficients(coefficients,3); - result = exp(coefficients[0]/(T+coefficients[1]) - coefficients[2]); + this->poly.checkCoefficients(coefficients,3); + result = exp(coefficients[0]/(x+coefficients[1]) - coefficients[2]); } else if (n==2) { - result = exp(polyval(coefficients, T)); + result = exp(this->poly.polyval(coefficients, x)); } else { throw ValueError(format("There is no function defined for this input (%d). ",n)); } @@ -593,13 +907,13 @@ long double BasePolynomial::expval(std::vector const& coefficients, /// Evaluates an exponential function for the given coefficients /// @param coefficients vector containing the ordered coefficients -/// @param x long double value that represents the current input in the 1st dimension -/// @param T long double value that represents the current input in the 2nd dimension +/// @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 -long double BasePolynomial::expval(std::vector< std::vector > const& coefficients, long double x, long double T, int n){ - long double result = 0.; +double BaseExponential::expval(const std::vector< std::vector > &coefficients, double x, double y, int n){ + double result = 0.; if (n==2) { - result = exp(polyval(coefficients, x, T)); + result = exp(this->poly.polyval(coefficients, x, y)); } else { throw ValueError(format("There is no function defined for this input (%d). ",n)); } @@ -607,212 +921,38 @@ long double BasePolynomial::expval(std::vector< std::vector > const } - - - - -/** Here are the real implementations, use these classes if possible. - * It is not the most flexible way, but the classes below include - * the polynomial solvers and should be easy to use. - */ -PolynomialImpl1D::Residual::Residual(PolynomialImpl1D *poly, long double y) { - this->poly=poly; - this->y=y; -} - -virtual double PolynomialImpl1D::Residual::call(double x) { - return poly->eval(x) - y; -} - -virtual double PolynomialImpl1D::Residual::deriv(double x) { - return poly->deriv(x); -} - -virtual double PolynomialImpl1D::ResidualInt::call(double x) { - return poly->integ(x) - y; -} - -virtual double PolynomialImpl1D::ResidualInt::deriv(double x) { - return poly->eval(x); -} - -virtual double PolynomialImpl1D::ResidualDer::call(double x) { - return poly->deriv(x) - y; -} - -virtual double PolynomialImpl1D::ResidualDer::deriv(double x) { - throw CoolProp::NotImplementedError("Second derivatives have not been implemented."); - // TODO: implement call of derivative function - return 0.0; -} - -PolynomialImpl1D::PolynomialImpl1D(const std::vector &coefficients){ - this->coefficients=coefficients; -} - -virtual long double PolynomialImpl1D::eval(long double x){ - return this->polyval(this->coefficients, x); -} - -virtual long double PolynomialImpl1D::integ(long double x){ - return this->polyint(this->coefficients, x); -} - -virtual long double PolynomialImpl1D::deriv(long double x){ - return this->polyder(this->coefficients, x); -} - -virtual long double PolynomialImpl1D::solve(long double y, long double x0){ - Residual resid(this, y); - std::string errstring; - return Newton(resid, x0, DBL_EPSILON*1e3, 100, errstring); -} - -virtual long double PolynomialImpl1D::solve(long double y, long double xmin, long double xmax){ - Residual resid(this, y); - std::string errstring; - return Brent(resid,xmin,xmax,DBL_EPSILON, DBL_EPSILON*1e3,100,errstring); -} - -virtual long double PolynomialImpl1D::solveInt(long double y, long double x0){ - ResidualInt resid(this, y); - std::string errstring; - return Newton(resid, x0, DBL_EPSILON*1e3, 100, errstring); -} - -virtual long double PolynomialImpl1D::solveInt(long double y, long double xmin, long double xmax){ - ResidualInt resid(this, y); - std::string errstring; - return Brent(resid,xmin,xmax,DBL_EPSILON, DBL_EPSILON*1e3,100,errstring); -} - -virtual long double PolynomialImpl1D::solveDer(long double y, long double x0){ - ResidualDer resid(this, y); - std::string errstring; - return Newton(resid, x0, DBL_EPSILON*1e3, 100, errstring); -} - -virtual long double PolynomialImpl1D::solveDer(long double y, long double xmin, long double xmax){ - ResidualDer resid(this, y); - std::string errstring; - return Brent(resid,xmin,xmax,DBL_EPSILON, DBL_EPSILON*1e3,100,errstring); } -PolynomialImpl2D::Residual::Residual(PolynomialImpl2D *poly, long double y, long double x) { - this->poly=poly; - this->y=y; - this->x=x; -} - -virtual double PolynomialImpl2D::Residual::call(double z) { - return poly->eval(x,z) - y; -} - -virtual double PolynomialImpl2D::Residual::deriv(double z) { - return poly->deriv(x,z); -} - -virtual double PolynomialImpl2D::ResidualInt::call(double z) { - return poly->integ(x,z) - y; -} - -virtual double PolynomialImpl2D::ResidualInt::deriv(double z) { - return poly->eval(x,z); -} - -virtual double PolynomialImpl2D::ResidualDer::call(double z) { - return poly->deriv(x,z) - y; -} - -virtual double PolynomialImpl2D::ResidualDer::deriv(double z) { - throw CoolProp::NotImplementedError("Second derivatives have not been implemented."); - // TODO: implement call of derivative function - return 0.0; -} - -PolynomialImpl2D::PolynomialImpl2D(const std::vector< std::vector > &coefficients){ - this->coefficients=coefficients; -} - -virtual long double PolynomialImpl2D::eval(long double x, long double z){ - return this->polyval(this->coefficients, x, z); -} - -virtual long double PolynomialImpl2D::integ(long double x, long double z){ - return this->polyint(this->coefficients, x, z); -} - -virtual long double PolynomialImpl2D::deriv(long double x, long double z){ - return this->polyder(this->coefficients, x, z); -} - -virtual long double PolynomialImpl2D::solve(long double y, long double x, long double z0){ - Residual resid(this, y, x); - std::string errstring; - return Newton(resid, z0, DBL_EPSILON*1e3, 100, errstring); -} - -virtual long double PolynomialImpl2D::solve(long double y, long double x, long double zmin, long double zmax){ - Residual resid(this, y, x); - std::string errstring; - return Brent(resid,zmin,zmax,DBL_EPSILON, DBL_EPSILON*1e3,100,errstring); -} - -virtual long double PolynomialImpl2D::solveInt(long double y, long double x, long double z0){ - ResidualInt resid(this, y, x); - std::string errstring; - return Newton(resid, z0, DBL_EPSILON*1e3, 100, errstring); -} - -virtual long double PolynomialImpl2D::solveInt(long double y, long double x, long double zmin, long double zmax){ - ResidualInt resid(this, y, x); - std::string errstring; - return Brent(resid,zmin,zmax,DBL_EPSILON, DBL_EPSILON*1e3,100,errstring); -} - -virtual long double PolynomialImpl2D::solveDer(long double y, long double x, long double z0){ - ResidualDer resid(this, y, x); - std::string errstring; - return Newton(resid, z0, DBL_EPSILON*1e3, 100, errstring); -} - -virtual long double PolynomialImpl2D::solveDer(long double y, long double x, long double zmin, long double zmax){ - ResidualDer resid(this, y, x); - std::string errstring; - return Brent(resid,zmin,zmax,DBL_EPSILON, DBL_EPSILON*1e3,100,errstring); -} - - - - - - - - - - - - -} - +// //int main() { // -// SimpleIncompressible* liquid = new DowthermQClass(); -// long double AT = 150.0 + 273.15; -// long double Ap = 3e5; -// liquid->testInputs(AT,Ap); +// 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(); // -// SecCoolSolution* obj = new MethanolSolution(); -// long double x = 0.25; -// long double T = 5.0 + 273.15; -// long double p = 3e5; +// double T = 273.15+50; // -// obj->testInputs(T+00,p,x); -// obj->testInputs(T+05,p,x); -// obj->testInputs(T+10,p,x); -// obj->testInputs(T+15,p,x); +// 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); // //}