From ed35139d8a25c915f34a0104328b8dec968ddb07 Mon Sep 17 00:00:00 2001 From: jowr Date: Thu, 10 Jul 2014 12:08:09 +0200 Subject: [PATCH] Implement new, more general solvers for the polynomials. Added solvers for rho, c, u and s to the incompressibles. T as function of h is still missing, should got to the backend... --- include/IncompressibleFluid.h | 24 ++++ include/PolyMath.h | 96 +++++++++++++- .../Incompressible/IncompressibleBackend.cpp | 52 +++++++- .../Incompressible/IncompressibleFluid.cpp | 83 +++++++++++- src/PolyMath.cpp | 120 +++++++++++++----- 5 files changed, 330 insertions(+), 45 deletions(-) diff --git a/include/IncompressibleFluid.h b/include/IncompressibleFluid.h index eaf7a8a2..732ca953 100644 --- a/include/IncompressibleFluid.h +++ b/include/IncompressibleFluid.h @@ -165,6 +165,30 @@ public: /// Conversion from mass-based to mole-based composition. double M2M (double T, double x); + /* Some functions can be inverted directly, those are listed + * here. It is also possible to solve for other quantities, but + * that involves some more sophisticated processing and is not + * done here, but in the backend, T(h,p) for example. + */ + /// Temperature as a function of density, pressure and composition. + double T_rho (double Dmass, double p, double x=0.0); + /// Temperature as a function of heat capacities as a function of temperature, pressure and composition. + double T_c (double Cmass, double p, double x=0.0); + /// Temperature as a function of entropy as a function of temperature, pressure and composition. + double T_s (double Smass, double p, double x=0.0); + /// Temperature as a function of internal energy as a function of temperature, pressure and composition. + double T_u (double Umass, double p, double x=0.0); + /// Temperature as a function of enthalpy, pressure and composition. + double T_h (double Hmass, double p, double x=0.0){throw NotImplementedError(format("%s (%d): T from enthalpy is not implemented in the fluid, use the backend.",__FILE__,__LINE__));} + /// Viscosity as a function of temperature, pressure and composition. + double T_visc(double visc, double p, double x=0.0){throw NotImplementedError(format("%s (%d): T from viscosity is not implemented.",__FILE__,__LINE__));} + /// Thermal conductivity as a function of temperature, pressure and composition. + double T_cond(double cond, double p, double x=0.0){throw NotImplementedError(format("%s (%d): T from conductivity is not implemented.",__FILE__,__LINE__));} + /// Saturation pressure as a function of temperature and composition. + double T_psat(double psat, double x=0.0){throw NotImplementedError(format("%s (%d): T from psat is not implemented.",__FILE__,__LINE__));} + /// Composition as a function of freezing temperature and pressure. + double x_Tfreeze( double Tfreeze, double p){throw NotImplementedError(format("%s (%d): x from T_freeze is not implemented.",__FILE__,__LINE__));} + protected: diff --git a/include/PolyMath.h b/include/PolyMath.h index f407c691..aa248d8e 100644 --- a/include/PolyMath.h +++ b/include/PolyMath.h @@ -13,6 +13,9 @@ namespace CoolProp{ +// Just a forward declaration +class Poly2DResidual; + /// The base class for all Polynomials /** A clear and straight forward implementation of polynomial operations. Still * very basic, but serves its purpose. @@ -99,6 +102,19 @@ public: /// @param axis unsigned integer value that represents the axis to integrate for (0=x, 1=y) double integral(const Eigen::MatrixXd &coefficients, const double &x_in, const double &y_in, const int &axis); +protected: + /// Uses the Brent solver to find the roots of p(x_in,y_in)-z_in + /// @param res Poly2DResidual object to calculate residuals and derivatives + /// @param min double value that represents the minimum value + /// @param max double value that represents the maximum value + double solve_limits(Poly2DResidual res, const double &min, const double &max); + + /// Uses the Newton solver to find the roots of p(x_in,y_in)-z_in + /// @param res Poly2DResidual object to calculate residuals and derivatives + /// @param guess double value that represents the start value + double solve_guess(Poly2DResidual res, const double &guess); + +public: /// Returns a vector with ALL the real roots of p(x_in,y_in)-z_in /// @param coefficients vector containing the ordered coefficients /// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension) @@ -163,6 +179,12 @@ protected: Poly2DResidual(); public: + /// Residual of a polynomial + /// @param Polynomial2D &poly polynomial object used to evaluate the calls + /// @param const Eigen::MatrixXd &coefficients, + /// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension) + /// @param z_in double value that represents the current output in the 3rd dimension + /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) Poly2DResidual(Polynomial2D &poly, const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis); virtual ~Poly2DResidual(){}; @@ -274,7 +296,7 @@ public: /// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension /// @param x_base double value that represents the base value for a centred fit in the 1st dimension /// @param y_base double value that represents the base value for a centred fit in the 2nd dimension - Eigen::VectorXd solve(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base); + Eigen::VectorXd solve(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis, const int &x_exp, const int &y_exp, const double &x_base=0.0, const double &y_base=0.0); /// Uses the Brent solver to find the roots of p(x_in,y_in)-z_in /// @param coefficients vector containing the ordered coefficients @@ -287,7 +309,7 @@ public: /// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension /// @param x_base double value that represents the base value for a centred fit in the 1st dimension /// @param y_base double value that represents the base value for a centred fit in the 2nd dimension - double solve_limits(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &min, const double &max, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base); + double solve_limits(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &min, const double &max, const int &axis, const int &x_exp, const int &y_exp, const double &x_base=0.0, const double &y_base=0.0); /// Uses the Newton solver to find the roots of p(x_in,y_in)-z_in /// @param coefficients vector containing the ordered coefficients @@ -299,7 +321,34 @@ public: /// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension /// @param x_base double value that represents the base value for a centred fit in the 1st dimension /// @param y_base double value that represents the base value for a centred fit in the 2nd dimension - double solve_guess(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &guess, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base); + double solve_guess(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &guess, const int &axis, const int &x_exp, const int &y_exp, const double &x_base=0.0, const double &y_base=0.0); + + /// Uses the Brent solver to find the roots of Int(p(x_in,y_in))-z_in + /// @param coefficients vector containing the ordered coefficients + /// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension) + /// @param z_in double value that represents the current output in the 3rd dimension + /// @param min double value that represents the minimum value + /// @param max double value that represents the maximum value + /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) + /// @param x_exp integer value that represents the lowest exponent of the polynomial in the 1st dimension + /// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension + /// @param x_base double value that represents the base value for a centred fit in the 1st dimension + /// @param y_base double value that represents the base value for a centred fit in the 2nd dimension + /// @param int_axis axis for the integration (0=x, 1=y) + double solve_limitsInt(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &min, const double &max, const int &axis, const int &x_exp, const int &y_exp, const double &x_base=0.0, const double &y_base=0.0, const int &int_axis=0); + + /// Uses the Newton solver to find the roots of Int(p(x_in,y_in))-z_in + /// @param coefficients vector containing the ordered coefficients + /// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension) + /// @param z_in double value that represents the current output in the 3rd dimension + /// @param guess double value that represents the start value + /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) + /// @param x_exp integer value that represents the lowest exponent of the polynomial in the 1st dimension + /// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension + /// @param x_base double value that represents the base value for a centred fit in the 1st dimension + /// @param y_base double value that represents the base value for a centred fit in the 2nd dimension + /// @param int_axis axis for the integration (0=x, 1=y) + double solve_guessInt(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &guess, const int &axis, const int &x_exp, const int &y_exp, const double &x_base=0.0, const double &y_base=0.0, const int &int_axis=0); protected: /// @param nValue integer value that represents the order of the factorial @@ -334,12 +383,51 @@ protected: Poly2DFracResidual(); public: - Poly2DFracResidual(Polynomial2DFrac &poly, const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base); + /// Residual of a polynomial divided by the independent variable + /// @param Polynomial2DFrac &poly polynomial object used to evaluate the calls + /// @param const Eigen::MatrixXd &coefficients, + /// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension) + /// @param z_in double value that represents the current output in the 3rd dimension + /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) + /// @param x_exp first exponent in x-direction + /// @param y_exp first exponent in y-direction + /// @param x_base base value for x (x = x_in - x_base) + /// @param y_base base value for y (y = y_in - y_base) + Poly2DFracResidual(Polynomial2DFrac &poly, const Eigen::MatrixXd &coefficients, + const double &in, const double &z_in, const int &axis, + const int &x_exp, const int &y_exp, const double &x_base, const double &y_base); virtual ~Poly2DFracResidual(){}; double call(double target); double deriv(double target); }; +class Poly2DFracIntResidual : public Poly2DFracResidual { + +protected: + int int_axis; + Poly2DFracIntResidual(); + +public: + /// Residual of an integrated polynomial divided by the independent variable + /// @param Polynomial2DFrac &poly polynomial object used to evaluate the calls + /// @param const Eigen::MatrixXd &coefficients, + /// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension) + /// @param z_in double value that represents the current output in the 3rd dimension + /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) + /// @param x_exp first exponent in x-direction + /// @param y_exp first exponent in y-direction + /// @param x_base base value for x (x = x_in - x_base) + /// @param y_base base value for y (y = y_in - y_base) + /// @param int_axis axis for the integration (0=x, 1=y) + Poly2DFracIntResidual(Polynomial2DFrac &poly, const Eigen::MatrixXd &coefficients, + const double &in, const double &z_in, const int &axis, + const int &x_exp, const int &y_exp, const double &x_base, const double &y_base, + const int &int_axis); + virtual ~Poly2DFracIntResidual(){}; + double call(double target); + double deriv(double target); +}; + // diff --git a/src/Backends/Incompressible/IncompressibleBackend.cpp b/src/Backends/Incompressible/IncompressibleBackend.cpp index 0ac1afbe..29402e7e 100644 --- a/src/Backends/Incompressible/IncompressibleBackend.cpp +++ b/src/Backends/Incompressible/IncompressibleBackend.cpp @@ -27,9 +27,9 @@ IncompressibleBackend::IncompressibleBackend(const std::vector &com } void IncompressibleBackend::update(long input_pair, double value1, double value2) { - if (mass_fractions.empty()){ - throw ValueError("mass fractions have not been set"); - } + //if (mass_fractions.empty()){ + // throw ValueError("mass fractions have not been set"); + //} switch (input_pair) { case PT_INPUTS: { @@ -65,7 +65,7 @@ void IncompressibleBackend::update(long input_pair, double value1, double value2 @param mole_fractions The vector of mole fractions of the components */ void IncompressibleBackend::set_mole_fractions(const std::vector &mole_fractions) { - throw CoolProp::NotImplementedError("Cannot set mole fractions for incompressible fluid"); + throw NotImplementedError("Cannot set mole fractions for incompressible fluid"); } /// Set the mass fractions @@ -73,12 +73,54 @@ void IncompressibleBackend::set_mole_fractions(const std::vector &m @param mass_fractions The vector of mass fractions of the components */ void IncompressibleBackend::set_mass_fractions(const std::vector &mass_fractions) { + if (mass_fractions.size()!=1) throw ValueError(format("The incompressible backend only supports one entry in the mass fraction vector and not %d.",mass_fractions.size())); this->mass_fractions = mass_fractions; } /// Check if the mole fractions have been set, etc. void IncompressibleBackend::check_status() { - throw CoolProp::NotImplementedError("Cannot check status for incompressible fluid"); + throw NotImplementedError("Cannot check status for incompressible fluid"); } +///// Calculate T given pressure and density +///** +//@param rhomass The mass density in kg/m^3 +//@param p The pressure in Pa +//@returns T The temperature in K +//*/ +//long double IncompressibleBackend::DmassP_flash(long double rhomass, long double p){ +// +//} +///// Calculate T given pressure and enthalpy +///** +//@param hmass The mass enthalpy in J/kg +//@param p The pressure in Pa +//@returns T The temperature in K +//*/ +//long double IncompressibleBackend::HmassP_flash(long double hmass, long double p); +///// Calculate T given pressure and entropy +///** +//@param smass The mass entropy in J/kg/K +//@param p The pressure in Pa +//@returns T The temperature in K +//*/ +//long double IncompressibleBackend::PSmass_flash(long double p, long double smass); +// +///// Calculate T given pressure and internal energy +///** +//@param umass The mass internal energy in J/kg +//@param p The pressure in Pa +//@returns T The temperature in K +//*/ +//long double IncompressibleBackend::PUmass_flash(long double p, long double umass); +// + + + + + + + + + } diff --git a/src/Backends/Incompressible/IncompressibleFluid.cpp b/src/Backends/Incompressible/IncompressibleFluid.cpp index 7fcec3bb..0d697cfb 100644 --- a/src/Backends/Incompressible/IncompressibleFluid.cpp +++ b/src/Backends/Incompressible/IncompressibleFluid.cpp @@ -105,7 +105,7 @@ double IncompressibleFluid::c (double T, double p, double x){ throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type)); break; default: - throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for entropy.",__FILE__,__LINE__,specific_heat.type)); + throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for specific heat.",__FILE__,__LINE__,specific_heat.type)); break; } return _HUGE; @@ -263,6 +263,87 @@ double IncompressibleFluid::V2M (double T, double y){throw NotImplemen /// Conversion from mass-based to mole-based composition. double IncompressibleFluid::M2M (double T, double x){throw NotImplementedError("TODO");} +/* Some functions can be inverted directly, those are listed + * here. It is also possible to solve for other quantities, but + * that involves some more sophisticated processing and is not + * done here, but in the backend, T(h,p) for example. + */ +/// Temperature as a function of density, pressure and composition. +double IncompressibleFluid::T_rho (double Dmass, double p, double x){ + double d_raw = Dmass; // No changes needed, no reference values... + switch (density.type) { + case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: + return poly.solve_limits(density.coeffs, x, d_raw, Tmin, Tmax, 0, 0, 0, Tbase, xbase); + break; + case IncompressibleData::INCOMPRESSIBLE_NOT_SET: + throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type)); + break; + default: + throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for inverse density.",__FILE__,__LINE__,specific_heat.type)); + break; + } + return _HUGE; +} +/// Temperature as a function of heat capacities as a function of temperature, pressure and composition. +double IncompressibleFluid::T_c (double Cmass, double p, double x){ + double c_raw = Cmass; // No changes needed, no reference values... + switch (specific_heat.type) { + case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: + return poly.solve_limits(specific_heat.coeffs, x, c_raw, Tmin, Tmax, 0, 0, 0, Tbase, xbase); + break; + case IncompressibleData::INCOMPRESSIBLE_NOT_SET: + throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type)); + break; + default: + throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for inverse specific heat.",__FILE__,__LINE__,specific_heat.type)); + break; + } + return _HUGE; +} +/// Temperature as a function of entropy as a function of temperature, pressure and composition. +double IncompressibleFluid::T_s (double Smass, double p, double x){ + double s_raw = Smass + sref; + switch (specific_heat.type) { + case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: + return poly.solve_limitsInt(specific_heat.coeffs, x, s_raw, Tmin, Tmax, 0, -1, 0, Tbase, xbase, 0); + break; + case IncompressibleData::INCOMPRESSIBLE_NOT_SET: + throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type)); + break; + default: + throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for inverse entropy.",__FILE__,__LINE__,specific_heat.type)); + break; + } + return _HUGE; +} +/// Temperature as a function of internal energy as a function of temperature, pressure and composition. +double IncompressibleFluid::T_u (double Umass, double p, double x){ + double u_raw = Umass + uref; + switch (specific_heat.type) { + case IncompressibleData::INCOMPRESSIBLE_POLYNOMIAL: + return poly.solve_limitsInt(specific_heat.coeffs, x, u_raw, Tmin, Tmax, 0, 0, 0, Tbase, xbase, 0); + break; + case IncompressibleData::INCOMPRESSIBLE_NOT_SET: + throw ValueError(format("%s (%d): The function type is not specified (\"[%d]\"), are you sure the coefficients have been set?",__FILE__,__LINE__,specific_heat.type)); + break; + default: + throw ValueError(format("%s (%d): There is no predefined way to use this function type \"[%d]\" for inverse entropy.",__FILE__,__LINE__,specific_heat.type)); + break; + } + return _HUGE; +} +///// Temperature as a function of enthalpy, pressure and composition. +//double IncompressibleFluid::T_h (double Hmass, double p, double x){throw NotImplementedError(format("%s (%d): T from enthalpy is not implemented in the fluid, use the backend.",__FILE__,__LINE__));} +///// Viscosity as a function of temperature, pressure and composition. +//double IncompressibleFluid::T_visc(double visc, double p, double x){throw NotImplementedError(format("%s (%d): T from viscosity is not implemented.",__FILE__,__LINE__));} +///// Thermal conductivity as a function of temperature, pressure and composition. +//double IncompressibleFluid::T_cond(double cond, double p, double x){throw NotImplementedError(format("%s (%d): T from conductivity is not implemented.",__FILE__,__LINE__));} +///// Saturation pressure as a function of temperature and composition. +//double IncompressibleFluid::T_psat(double psat, double x){throw NotImplementedError(format("%s (%d): T from psat is not implemented.",__FILE__,__LINE__));} +///// Composition as a function of freezing temperature and pressure. +//double IncompressibleFluid::x_Tfreeze( double Tfreeze, double p){throw NotImplementedError(format("%s (%d): x from T_freeze is not implemented.",__FILE__,__LINE__));} + + /* * Some more functions to provide a single implementation * of important routines. diff --git a/src/PolyMath.cpp b/src/PolyMath.cpp index 902b75e9..3c6e1fa2 100644 --- a/src/PolyMath.cpp +++ b/src/PolyMath.cpp @@ -196,6 +196,33 @@ double Polynomial2D::integral(const Eigen::MatrixXd &coefficients, const double return this->evaluate(this->integrateCoeffs(coefficients, axis, 1), x_in,y_in); } +/// Uses the Brent solver to find the roots of p(x_in,y_in)-z_in +/// @param res Poly2DResidual object to calculate residuals and derivatives +/// @param min double value that represents the minimum value +/// @param max double value that represents the maximum value +double Polynomial2D::solve_limits(Poly2DResidual res, const double &min, const double &max){ + std::string errstring; + double macheps = DBL_EPSILON; + double tol = DBL_EPSILON*1e3; + int maxiter = 10; + double result = Brent(res, min, max, macheps, tol, maxiter, errstring); + if (this->do_debug()) std::cout << "Brent solver message: " << errstring << std::endl; + return result; +} + +/// Uses the Newton solver to find the roots of p(x_in,y_in)-z_in +/// @param res Poly2DResidual object to calculate residuals and derivatives +/// @param guess double value that represents the start value +double Polynomial2D::solve_guess(Poly2DResidual res, const double &guess){ + std::string errstring; + //set_debug_level(1000); + double tol = DBL_EPSILON*1e3; + int maxiter = 10; + double result = Newton(res, guess, tol, maxiter, errstring); + if (this->do_debug()) std::cout << "Newton solver message: " << errstring << std::endl; + return result; +} + /// @param coefficients vector containing the ordered coefficients /// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension) @@ -236,13 +263,7 @@ Eigen::VectorXd Polynomial2D::solve(const Eigen::MatrixXd &coefficients, const d /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) double Polynomial2D::solve_limits(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &min, const double &max, const int &axis){ Poly2DResidual res = Poly2DResidual(*this, coefficients, in, z_in, axis); - std::string errstring; - double macheps = DBL_EPSILON; - double tol = DBL_EPSILON*1e3; - int maxiter = 10; - double result = Brent(res, min, max, macheps, tol, maxiter, errstring); - if (this->do_debug()) std::cout << "Brent solver message: " << errstring << std::endl; - return result; + return solve_limits(res, min, max); } /// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension) @@ -251,13 +272,7 @@ double Polynomial2D::solve_limits(const Eigen::MatrixXd &coefficients, const dou /// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) double Polynomial2D::solve_guess(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &guess, const int &axis){ Poly2DResidual res = Poly2DResidual(*this, coefficients, in, z_in, axis); - std::string errstring; - //set_debug_level(1000); - double tol = DBL_EPSILON*1e3; - int maxiter = 10; - double result = Newton(res, guess, tol, maxiter, errstring); - if (this->do_debug()) std::cout << "Newton solver message: " << errstring << std::endl; - return result; + return solve_guess(res, guess); } @@ -650,7 +665,7 @@ double Polynomial2DFrac::integral(const Eigen::MatrixXd &coefficients, const dou /// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension /// @param x_base double value that represents the base value for a centred fit in the 1st dimension /// @param y_base double value that represents the base value for a centred fit in the 2nd dimension -Eigen::VectorXd Polynomial2DFrac::solve(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis, const int &x_exp, const int &y_exp, const double &x_base = 0.0, const double &y_base = 0.0){ +Eigen::VectorXd Polynomial2DFrac::solve(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base){ Eigen::MatrixXd newCoefficients; Eigen::VectorXd tmpCoefficients; @@ -715,16 +730,10 @@ Eigen::VectorXd Polynomial2DFrac::solve(const Eigen::MatrixXd &coefficients, con /// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension /// @param x_base double value that represents the base value for a centred fit in the 1st dimension /// @param y_base double value that represents the base value for a centred fit in the 2nd dimension -double Polynomial2DFrac::solve_limits(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &min, const double &max, const int &axis, const int &x_exp, const int &y_exp, const double &x_base = 0.0, const double &y_base = 0.0){ +double Polynomial2DFrac::solve_limits(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &min, const double &max, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base){ Poly2DFracResidual res = Poly2DFracResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base); - std::string errstring; - double macheps = DBL_EPSILON; - double tol = DBL_EPSILON*1e3; - int maxiter = 10; - double result = Brent(res, min, max, macheps, tol, maxiter, errstring); - if (this->do_debug()) std::cout << "Brent solver message: " << errstring << std::endl; - return result; -} + return Polynomial2D::solve_limits(res, min, max); +} //TODO: Implement tests for this solver /// Uses the Newton solver to find the roots of p(x_in,y_in)-z_in /// @param coefficients vector containing the ordered coefficients @@ -736,17 +745,43 @@ double Polynomial2DFrac::solve_limits(const Eigen::MatrixXd &coefficients, const /// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension /// @param x_base double value that represents the base value for a centred fit in the 1st dimension /// @param y_base double value that represents the base value for a centred fit in the 2nd dimension -double Polynomial2DFrac::solve_guess(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &guess, const int &axis, const int &x_exp, const int &y_exp, const double &x_base = 0.0, const double &y_base = 0.0){ +double Polynomial2DFrac::solve_guess(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &guess, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base){ Poly2DFracResidual res = Poly2DFracResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base); - std::string errstring; - //set_debug_level(1000); - double tol = DBL_EPSILON*1e3; - int maxiter = 10; - double result = Newton(res, guess, tol, maxiter, errstring); - if (this->do_debug()) std::cout << "Newton solver message: " << errstring << std::endl; - return result; -} + return Polynomial2D::solve_guess(res, guess); +} //TODO: Implement tests for this solver +/// Uses the Brent solver to find the roots of Int(p(x_in,y_in))-z_in +/// @param coefficients vector containing the ordered coefficients +/// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension) +/// @param z_in double value that represents the current output in the 3rd dimension +/// @param min double value that represents the minimum value +/// @param max double value that represents the maximum value +/// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) +/// @param x_exp integer value that represents the lowest exponent of the polynomial in the 1st dimension +/// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension +/// @param x_base double value that represents the base value for a centred fit in the 1st dimension +/// @param y_base double value that represents the base value for a centred fit in the 2nd dimension +/// @param int_axis axis for the integration (0=x, 1=y) +double Polynomial2DFrac::solve_limitsInt(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &min, const double &max, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base, const int &int_axis){ + Poly2DFracIntResidual res = Poly2DFracIntResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base, int_axis); + return Polynomial2D::solve_limits(res, min, max); +} //TODO: Implement tests for this solver + +/// Uses the Newton solver to find the roots of Int(p(x_in,y_in))-z_in +/// @param coefficients vector containing the ordered coefficients +/// @param in double value that represents the current input in x (1st dimension) or y (2nd dimension) +/// @param z_in double value that represents the current output in the 3rd dimension +/// @param guess double value that represents the start value +/// @param axis unsigned integer value that represents the axis to solve for (0=x, 1=y) +/// @param x_exp integer value that represents the lowest exponent of the polynomial in the 1st dimension +/// @param y_exp integer value that represents the lowest exponent of the polynomial in the 2nd dimension +/// @param x_base double value that represents the base value for a centred fit in the 1st dimension +/// @param y_base double value that represents the base value for a centred fit in the 2nd dimension +/// @param int_axis axis for the integration (0=x, 1=y) +double Polynomial2DFrac::solve_guessInt(const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const double &guess, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base, const int &int_axis){ + Poly2DFracIntResidual res = Poly2DFracIntResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base, int_axis); + return Polynomial2D::solve_guess(res, guess); +} //TODO: Implement tests for this solver @@ -822,16 +857,31 @@ Poly2DFracResidual::Poly2DFracResidual(Polynomial2DFrac &poly, const Eigen::Matr double Poly2DFracResidual::call(double target){ if (axis==iX) return poly.evaluate(coefficients,target,in,x_exp,y_exp,x_base,y_base)-z_in; if (axis==iY) return poly.evaluate(coefficients,in,target,x_exp,y_exp,x_base,y_base)-z_in; - return -_HUGE; + return _HUGE; } double Poly2DFracResidual::deriv(double target){ if (axis==iX) return poly.derivative(coefficients,target,in,axis,x_exp,y_exp,x_base,y_base); if (axis==iY) return poly.derivative(coefficients,in,target,axis,x_exp,y_exp,x_base,y_base); - return -_HUGE; + return _HUGE; } +Poly2DFracIntResidual::Poly2DFracIntResidual(Polynomial2DFrac &poly, const Eigen::MatrixXd &coefficients, const double &in, const double &z_in, const int &axis, const int &x_exp, const int &y_exp, const double &x_base, const double &y_base, const int &int_axis) + : Poly2DFracResidual(poly, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base){ + this->int_axis = int_axis; +} +double Poly2DFracIntResidual::call(double target){ + if (axis==iX) return poly.integral(coefficients,target,in,int_axis,x_exp,y_exp,x_base,y_base)-z_in; + if (axis==iY) return poly.integral(coefficients,in,target,int_axis,x_exp,y_exp,x_base,y_base)-z_in; + return _HUGE; +} + +double Poly2DFracIntResidual::deriv(double target){ + if (axis==iX) return poly.evaluate(coefficients,target,in,x_exp,y_exp,x_base,y_base); + if (axis==iY) return poly.evaluate(coefficients,in,target,x_exp,y_exp,x_base,y_base); + return _HUGE; +}