diff --git a/dev/incompressible_liquids/CPIncomp/BaseObjects.py b/dev/incompressible_liquids/CPIncomp/BaseObjects.py
index 6c0e7fc2..7507d285 100644
--- a/dev/incompressible_liquids/CPIncomp/BaseObjects.py
+++ b/dev/incompressible_liquids/CPIncomp/BaseObjects.py
@@ -27,7 +27,7 @@ class IncompressibleData(object):
         self.maxLog = np.log(np.finfo(np.float64).max-1)
         self.minLog = -self.maxLog
         
-        self.DEBUG = False
+        self.DEBUG = True
         
         
     ### Base functions that handle the custom data type, just a place holder to show the structure.
@@ -113,7 +113,7 @@ class IncompressibleData(object):
         return (r,c,np.reshape(array,(r,c)))
         
         
-    def fit(self, x, y=0.0, xbase=0.0, ybase=0.0):
+    def fit(self, x=0.0, y=0.0, xbase=0.0, ybase=0.0):
         """ 
         A function that selects the correct equations and
         fits coefficients. Some functions require a start
@@ -134,15 +134,23 @@ class IncompressibleData(object):
         
         cr,cc,_ = self.shapeArray(self.coeffs)
         if self.DEBUG: print("Coefficients: ({0},{1})".format(cr,cc))
+        if (xr==1 and xc==1 and cr>1):
+            if self.DEBUG: print("Discarding coefficient rows, {0} -> {1}".format(cr,xr))
+            self.coeffs = self.coeffs[0]
         if (yr==1 and yc==1 and cc>1):
-            if self.DEBUG: print("Discarding coefficient dimension, {0} -> {1}".format(cc,yc))
-            self.coeffs = self.coeffs.T[0]
+            if self.DEBUG: print("Discarding coefficient columns, {0} -> {1}".format(cc,yc))
+            self.coeffs = self.coeffs.T[0].T
         cr,cc,_ = self.shapeArray(self.coeffs)
-        
+                
         if self.DEBUG: print("Coefficients before fitting: \n{0}".format(self.coeffs))
         
         # Polynomial fitting works for both 1D and 2D functions
         if self.type==self.INCOMPRESSIBLE_POLYNOMIAL or self.type==self.INCOMPRESSIBLE_EXPPOLYNOMIAL:
+            if self.DEBUG: print("polynomial detected, fitting {0}".format(self.type))
+            if cr==1 and cc==1: 
+                if self.DEBUG: print("No coefficients left to fit, aborting procedure.")
+                self.coeffs = np.array([[0.0]])
+                return
             if (xr<cr): raise ValueError("Less data points than coefficients in first dimension ({0} < {1}), aborting.".format(xr,cr))
             if (yc<cc): raise ValueError("Less data points than coefficients in second dimension ({0} < {1}), aborting.".format(yc,cc))
             x_input = np.array(x.flat)-xbase
@@ -157,7 +165,7 @@ class IncompressibleData(object):
         
         # Select if 1D or 2D fitting
         if yc==1: # 1D fitting, only one input
-            if self.DEBUG: print("1D function detected, fitting {0}".format(self.type))
+            if self.DEBUG: print("1D function detected, fitting {0}".format(self.type))      
             x_input = x-xbase
             self.coeffs = self.getCoeffsIterative1D(x_input, coeffs_start=None)
             if self.DEBUG: print("Coefficients after fitting: \n{0}".format(self.coeffs))
diff --git a/dev/incompressible_liquids/CPIncomp/CoefficientObjects.py b/dev/incompressible_liquids/CPIncomp/CoefficientObjects.py
index 1c625808..3dc77588 100644
--- a/dev/incompressible_liquids/CPIncomp/CoefficientObjects.py
+++ b/dev/incompressible_liquids/CPIncomp/CoefficientObjects.py
@@ -238,12 +238,12 @@ class SecCoolExample(CoefficientData):
            4.482000E-09]))
 
         self.T_freeze.type = self.T_freeze.INCOMPRESSIBLE_POLYOFFSET
-        self.T_freeze.coeffs = np.array([[
+        self.T_freeze.coeffs = np.array([
            27.755555600/100.0,
           -22.973221700, 
           -1.1040507200*100.0, 
           -0.0120762281*100.0*100.0, 
-          -9.343458E-05*100.0*100.0*100.0]])  
+          -9.343458E-05*100.0*100.0*100.0])  
         
         
 class MelinderExample(CoefficientData):
diff --git a/dev/incompressible_liquids/CPIncomp/WriterObjects.py b/dev/incompressible_liquids/CPIncomp/WriterObjects.py
index 71440ae7..9a735cb5 100644
--- a/dev/incompressible_liquids/CPIncomp/WriterObjects.py
+++ b/dev/incompressible_liquids/CPIncomp/WriterObjects.py
@@ -69,9 +69,9 @@ class SolutionDataWriter(object):
             pass
 
         try:        
-            data.T_freeze.coeffs = np.copy(std_coeffs)[0]
+            data.T_freeze.coeffs = np.copy(std_coeffs)
             data.T_freeze.type   = data.T_freeze.INCOMPRESSIBLE_POLYNOMIAL
-            data.T_freeze.fit(x,0.0,data.xbase,0.0)
+            data.T_freeze.fit(0.0,x,0.0,data.xbase)
         except errList as ve:
             if self.verbose: print("Could not fit TFreeze coefficients:", ve)
             pass
diff --git a/dev/incompressible_liquids/ExampleMelinder.json b/dev/incompressible_liquids/ExampleMelinder.json
index 0a54cfc3..18ca7228 100644
--- a/dev/incompressible_liquids/ExampleMelinder.json
+++ b/dev/incompressible_liquids/ExampleMelinder.json
@@ -122,12 +122,8 @@
   "name": "ExampleMelinder", 
   "reference": "Melinder Book", 
   "saturation_pressure": {
-    "coeffs": [
-      -5.000000000000000e+03, 
-       3.000000000000000e+01, 
-      -1.000000000000000e+01
-    ], 
-    "type": "exponential"
+    "coeffs": "null", 
+    "type": "notdefined"
   }, 
   "specific_heat": {
     "coeffs": [
diff --git a/dev/incompressible_liquids/ExamplePure.json b/dev/incompressible_liquids/ExamplePure.json
index 9906151c..928423b9 100644
--- a/dev/incompressible_liquids/ExamplePure.json
+++ b/dev/incompressible_liquids/ExamplePure.json
@@ -1,11 +1,9 @@
 {
   "T_freeze": {
     "coeffs": [
-       0.000000000000000e+00, 
-       0.000000000000000e+00, 
-       0.000000000000000e+00, 
-       0.000000000000000e+00, 
-       0.000000000000000e+00
+      [
+         0.000000000000000e+00
+      ]
     ], 
     "type": "polynomial"
   }, 
@@ -73,10 +71,10 @@
   "saturation_pressure": {
     "coeffs": [
       [
-        -3.373625698981663e+01
+        -3.373625698981662e+01
       ], 
       [
-         1.484757141205828e-01
+         1.484757141205827e-01
       ], 
       [
         -1.429827445582900e-04
diff --git a/dev/incompressible_liquids/ExampleSecCool.json b/dev/incompressible_liquids/ExampleSecCool.json
index 4187080e..b8d7e828 100644
--- a/dev/incompressible_liquids/ExampleSecCool.json
+++ b/dev/incompressible_liquids/ExampleSecCool.json
@@ -1,13 +1,11 @@
 {
   "T_freeze": {
     "coeffs": [
-      [
-         2.775555560000000e-01, 
-        -2.297322170000000e+01, 
-        -1.104050720000000e+02, 
-        -1.207622810000000e+02, 
-        -9.343458000000001e+01
-      ]
+       2.775555560000000e-01, 
+      -2.297322170000000e+01, 
+      -1.104050720000000e+02, 
+      -1.207622810000000e+02, 
+      -9.343458000000001e+01
     ], 
     "type": "polyoffset"
   }, 
@@ -97,12 +95,8 @@
   "name": "ExampleSecCool", 
   "reference": "SecCool software", 
   "saturation_pressure": {
-    "coeffs": [
-      -5.000000000000000e+03, 
-       3.000000000000000e+01, 
-      -1.000000000000000e+01
-    ], 
-    "type": "exponential"
+    "coeffs": "null", 
+    "type": "notdefined"
   }, 
   "specific_heat": {
     "coeffs": [
diff --git a/dev/incompressible_liquids/ExampleSolution.json b/dev/incompressible_liquids/ExampleSolution.json
index 9db098fb..831a672e 100644
--- a/dev/incompressible_liquids/ExampleSolution.json
+++ b/dev/incompressible_liquids/ExampleSolution.json
@@ -1,11 +1,27 @@
 {
   "T_freeze": {
     "coeffs": [
-       0.000000000000000e+00, 
-       0.000000000000000e+00, 
-       0.000000000000000e+00, 
-       0.000000000000000e+00, 
-       0.000000000000000e+00
+      [
+         0.000000000000000e+00, 
+         0.000000000000000e+00, 
+         0.000000000000000e+00, 
+         0.000000000000000e+00, 
+         0.000000000000000e+00
+      ], 
+      [
+         0.000000000000000e+00, 
+         0.000000000000000e+00, 
+         0.000000000000000e+00, 
+         0.000000000000000e+00, 
+         0.000000000000000e+00
+      ], 
+      [
+         0.000000000000000e+00, 
+         0.000000000000000e+00, 
+         0.000000000000000e+00, 
+         0.000000000000000e+00, 
+         0.000000000000000e+00
+      ]
     ], 
     "type": "polynomial"
   }, 
@@ -42,23 +58,23 @@
   "density": {
     "coeffs": [
       [
-        -2.399430982040806e+02, 
-         1.250815901394219e+03, 
-         4.235073668358468e+01, 
-        -1.816072631001721e+02, 
-         3.787878787861905e+01
+        -2.399430982040804e+02, 
+         1.250815901394215e+03, 
+         4.235073668359494e+01, 
+        -1.816072631001761e+02, 
+         3.787878787862606e+01
       ], 
       [
          1.185772755113077e+01, 
-        -1.350134561596798e+01, 
-         7.888548752001322e-01, 
-         5.820105820111443e-01, 
+        -1.350134561596795e+01, 
+         7.888548752001263e-01, 
+         5.820105820110935e-01, 
          0.000000000000000e+00
       ], 
       [
-        -2.674489795934265e-02, 
-         3.140022675723779e-02, 
-        -3.287981859249434e-03, 
+        -2.674489795948098e-02, 
+         3.140022675718737e-02, 
+        -3.287981859254933e-03, 
          0.000000000000000e+00, 
          0.000000000000000e+00
       ]
diff --git a/include/IncompressibleFluid.h b/include/IncompressibleFluid.h
index ac0237eb..35d252a3 100644
--- a/include/IncompressibleFluid.h
+++ b/include/IncompressibleFluid.h
@@ -183,11 +183,11 @@ public:
 	/// 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__));}
+	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__));}
+	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__));}
+	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__));}
 
diff --git a/include/MatrixMath.h b/include/MatrixMath.h
index 1fd5b305..a9500edc 100644
--- a/include/MatrixMath.h
+++ b/include/MatrixMath.h
@@ -141,6 +141,42 @@ template <class T> Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> vec_to_eigen(c
 }
 
 
+/// Convert 1D matrix to vector
+/** Returns either a row- or a column-based
+ *  vector. By default, Eigen prefers column
+ *  major ordering, just like Fortran.
+ */
+template< class T> Eigen::Matrix<T,Eigen::Dynamic,1> makeVector(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &matrix) {
+	return makeColVector(matrix);
+}
+template< class T> Eigen::Matrix<T,Eigen::Dynamic,1> makeColVector(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &matrix){
+	std::size_t r = matrix.rows();
+	std::size_t c = matrix.cols();
+	Eigen::Matrix<T,Eigen::Dynamic,1> vector;
+	if (r==1&&c>=1) { // Check passed, matrix can be transformed
+		vector = matrix.transpose().block(0,0,c,r);
+	} else if ( r>=1&&c==1) { // Check passed, matrix can be transformed
+		vector = matrix.block(0,0,r,c);
+	} else { // Check failed, throw error
+		throw ValueError(format("Your matrix (%d,%d) cannot be converted into a vector (x,1).",r,c));
+	}
+	return vector;
+}
+template< class T> Eigen::Matrix<T,1,Eigen::Dynamic> makeRowVector(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &matrix){
+	std::size_t r = matrix.rows();
+	std::size_t c = matrix.cols();
+	Eigen::Matrix<T,1,Eigen::Dynamic> vector;
+	if (r==1&&c>=1) { // Check passed, matrix can be transformed
+		vector = matrix.block(0,0,r,c);
+	} else if ( r>=1&&c==1) { // Check passed, matrix can be transformed
+		vector = matrix.transpose().block(0,0,c,r);
+	} else { // Check failed, throw error
+		throw ValueError(format("Your matrix (%d,%d) cannot be converted into a vector (1,x).",r,c));
+	}
+	return vector;
+}
+
+
 /// Remove rows and columns from matrices
 /** A set of convenience functions inspired by http://stackoverflow.com/questions/13290395/how-to-remove-a-certain-row-or-column-while-using-eigen-library-c
  *  but altered to respect templates.
diff --git a/src/Backends/Incompressible/IncompressibleFluid.cpp b/src/Backends/Incompressible/IncompressibleFluid.cpp
index ec780fa3..4ea25566 100644
--- a/src/Backends/Incompressible/IncompressibleFluid.cpp
+++ b/src/Backends/Incompressible/IncompressibleFluid.cpp
@@ -40,14 +40,16 @@ bool IncompressibleFluid::is_pure() {
 
 /// Base functions that handle the custom data type, just a place holder to show the structure.
 double IncompressibleFluid::baseExponential(IncompressibleData data, double y, double ybase){
-	size_t r=data.coeffs.rows(),c=data.coeffs.cols();
+	Eigen::VectorXd coeffs = makeVector(data.coeffs);
+	size_t r=coeffs.rows(),c=coeffs.cols();
 	if (strict && (r!=3 || c!=1) ) throw ValueError(format("%s (%d): You have to provide a 3,1 matrix of coefficients, not  (%d,%d).",__FILE__,__LINE__,r,c));
-	return exp( (double) (data.coeffs(0,0) / ( (y-ybase)+data.coeffs(1,0) ) - data.coeffs(2,0) ) );
+	return exp( (double) (coeffs[0] / ( (y-ybase)+coeffs[1] ) - coeffs[2] ) );
 }
 double IncompressibleFluid::baseExponentialOffset(IncompressibleData data, double y){
-	size_t r=data.coeffs.rows(),c=data.coeffs.cols();
+	Eigen::VectorXd coeffs = makeVector(data.coeffs);
+	size_t r=coeffs.rows(),c=coeffs.cols();
 	if (strict && (r!=4 || c!=1) ) throw ValueError(format("%s (%d): You have to provide a 4,1 matrix of coefficients, not  (%d,%d).",__FILE__,__LINE__,r,c));
-	return exp( (double) (data.coeffs(1,0) / ( (y-data.coeffs(0,0))+data.coeffs(2,0) ) - data.coeffs(3,0) ) );
+	return exp( (double) (coeffs[1] / ( (y-coeffs[0])+coeffs[2] ) - coeffs[3] ) );
 }
 double IncompressibleFluid::basePolyOffset(IncompressibleData data, double y, double z){
 	size_t r=data.coeffs.rows(),c=data.coeffs.cols();
@@ -60,7 +62,7 @@ double IncompressibleFluid::basePolyOffset(IncompressibleData data, double y, do
 			coeffs = Eigen::MatrixXd(data.coeffs.block(0,1,r,c-1));
 			in = z;
 		} else if (r>1 && c==1) { // column vector -> function of y
-			coeffs = Eigen::MatrixXd(data.coeffs.block(1,0,r-1,c)).transpose();
+			coeffs = Eigen::MatrixXd(data.coeffs.block(1,0,r-1,c));
 			in = y;
 		} else {
 			throw ValueError(format("%s (%d): You have to provide a vector (1D matrix) of coefficients, not  (%d,%d).",__FILE__,__LINE__,r,c));
diff --git a/src/Backends/Incompressible/IncompressibleLibrary.cpp b/src/Backends/Incompressible/IncompressibleLibrary.cpp
index 1bb0bff6..1fa4ed50 100644
--- a/src/Backends/Incompressible/IncompressibleLibrary.cpp
+++ b/src/Backends/Incompressible/IncompressibleLibrary.cpp
@@ -35,7 +35,7 @@ IncompressibleData JSONIncompressibleLibrary::parse_coefficients(rapidjson::Valu
 				}
 				else if (!type.compare("polyoffset")){
 					fluidData.type = CoolProp::IncompressibleData::INCOMPRESSIBLE_POLYOFFSET;
-					fluidData.coeffs = vec_to_eigen(cpjson::get_double_array2D(obj[id.c_str()]["coeffs"]));
+					fluidData.coeffs = vec_to_eigen(cpjson::get_double_array(obj[id.c_str()]["coeffs"]));
 					return fluidData;
 				}
 				else if (vital){
@@ -95,10 +95,11 @@ void JSONIncompressibleLibrary::add_one(rapidjson::Value &fluid_json) {
 
 	// Create an instance of the fluid
 	IncompressibleFluid &fluid = fluid_map[index];
+	fluid.setName("unloaded");
     try
     {
-
 	    fluid.setName(cpjson::get_string(fluid_json, "name"));
+    	if (get_debug_level()>=20) std::cout << format("Incompressible library: Loading base values for %s ",fluid.getName().c_str()) << std::endl;
 	    fluid.setDescription(cpjson::get_string(fluid_json, "description"));
 	    fluid.setReference(cpjson::get_string(fluid_json, "reference"));
 	    fluid.setTmax(parse_value(fluid_json, "Tmax", true, 0.0));
@@ -111,6 +112,7 @@ void JSONIncompressibleLibrary::add_one(rapidjson::Value &fluid_json) {
 	    fluid.setxbase(parse_value(fluid_json, "xbase", false, 0.0));
 
 	    /// Setters for the coefficients
+	    if (get_debug_level()>=20) std::cout << format("Incompressible library: Loading coefficients for %s ",fluid.getName().c_str()) << std::endl;
 	    fluid.setDensity(parse_coefficients(fluid_json, "density", true));
 	    fluid.setSpecificHeat(parse_coefficients(fluid_json, "specific_heat", true));
 	    fluid.setViscosity(parse_coefficients(fluid_json, "viscosity", false));
@@ -120,6 +122,7 @@ void JSONIncompressibleLibrary::add_one(rapidjson::Value &fluid_json) {
 	    fluid.setVolToMass(parse_coefficients(fluid_json, "volume2mass", false));
 	    fluid.setMassToMole(parse_coefficients(fluid_json, "mass2mole", false));
 
+	    if (get_debug_level()>=20) std::cout << format("Incompressible library: Loading reference state for %s ",fluid.getName().c_str()) << std::endl;
 	    fluid.set_reference_state(
 			    parse_value(fluid_json, "Tref", false, 25+273.15) ,
 			    parse_value(fluid_json, "pref", false, 1.01325e5) ,
diff --git a/src/CoolProp.cpp b/src/CoolProp.cpp
index 080fa4ed..2bb7509d 100644
--- a/src/CoolProp.cpp
+++ b/src/CoolProp.cpp
@@ -295,7 +295,7 @@ bool has_solution_concentration(const std::string &fluid_string)
     // Start at the "::" if it is found to chop off the delimiter between backend and fluid
     std::size_t i = fluid_string.find("::");
     // If can find "-", expect mass fractions encoded as string
-    if (fluid_string.find('-',i+2))
+    if (fluid_string.find('-',i+2)!=std::string::npos)
     	return true;
     return false;
 }
@@ -425,6 +425,7 @@ double PropsSI(const std::string &Output, const std::string &Name1, double Prop1
         }
         else
         {
+        	if (get_debug_level()>10) std::cout << format("%s:%d: Filling fractions for %s with [1.0].",__FILE__,__LINE__,Ref.c_str()) << std::endl;
             // Here it must be a pure fluid since the fractions are not coded in the string.  If not, it will fail internally
             return _PropsSI(Output, Name1, Prop1, Name2, Prop2, Ref, std::vector<double>(1,1));
         }
diff --git a/src/PolyMath.cpp b/src/PolyMath.cpp
index 24f68051..257d6e37 100644
--- a/src/PolyMath.cpp
+++ b/src/PolyMath.cpp
@@ -157,10 +157,7 @@ Eigen::MatrixXd Polynomial2D::deriveCoeffs(const Eigen::MatrixXd &coefficients,
 /// @param coefficients vector containing the ordered coefficients
 /// @param x_in double value that represents the current input
 double Polynomial2D::evaluate(const Eigen::MatrixXd &coefficients, const double &x_in){
-	if (coefficients.rows() != 1) {
-		throw ValueError(format("%s (%d): You have a 2D coefficient matrix (%d,%d), please use the 2D functions. ",__FILE__,__LINE__,coefficients.rows(),coefficients.cols()));
-	}
-	double result = Eigen::poly_eval( Eigen::RowVectorXd(coefficients), x_in );
+	double result = Eigen::poly_eval( makeVector(coefficients), x_in );
 	if (this->do_debug()) std::cout << "Running      1D evaluate(" << mat_to_string(coefficients) << ", x_in:" << vec_to_string(x_in) << "): " << result << std::endl;
 	return result;
 }
@@ -243,7 +240,7 @@ Eigen::VectorXd Polynomial2D::solve(const Eigen::MatrixXd &coefficients, const d
 	case 1:
 		tmpCoefficients = Eigen::VectorXd::Zero(c);
 		for(size_t i=0; i<c; i++) {
-			tmpCoefficients(i,0) = evaluate(coefficients.col(i).transpose(), in);
+			tmpCoefficients(i,0) = evaluate(coefficients.col(i), in);
 		}
 		break;
 	default:
@@ -454,16 +451,23 @@ Eigen::MatrixXd Polynomial2DFrac::deriveCoeffs(const Eigen::MatrixXd &coefficien
 /// @param firstExponent integer value that represents the lowest exponent of the polynomial
 /// @param x_base double value that represents the base value for a centred fit in the 1st dimension
 double Polynomial2DFrac::evaluate(const Eigen::MatrixXd &coefficients, const double &x_in, const int &firstExponent, const double &x_base){
-	if (coefficients.rows() != 1) {
+	size_t r = coefficients.rows();
+	size_t c = coefficients.cols();
+
+	if ( (r!=1) && (c!=1) ) {
 		throw ValueError(format("%s (%d): You have a 2D coefficient matrix (%d,%d), please use the 2D functions. ",__FILE__,__LINE__,coefficients.rows(),coefficients.cols()));
 	}
 	if ( (firstExponent<0) && (fabs(x_in-x_base)<DBL_EPSILON)) {
 		throw ValueError(format("%s (%d): A fraction cannot be evaluated with zero as denominator, x_in-x_base=%f ",__FILE__,__LINE__,x_in-x_base));
 	}
+
 	Eigen::MatrixXd tmpCoeffs(coefficients);
+	if ( c==1 ) {
+		tmpCoeffs.transposeInPlace();
+		c = r;
+		r = 1;
+	}
 	Eigen::MatrixXd newCoeffs;
-	size_t r=1;
-	size_t c;
 	double negExp = 0;// First we treat the negative exponents
 	double posExp = 0;// then the positive exponents
 
@@ -1367,1437 +1371,3 @@ TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp","
 
 #endif /* ENABLE_CATCH */
 
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/// Constructors for the base class for all Polynomials
-//Polynomial1D::Polynomial1D();
-//bool Polynomial2D::setCoefficients(const Eigen::MatrixXd &coefficients){
-//	this.coefficients = coefficients;
-//	return this.coefficients == coefficients;
-//}
-//bool Polynomial2D::setCoefficients(const std::vector< std::vector<double> > &coefficients){
-//	return this->setCoefficients(convert(coefficients));
-//}
-
-
-
-
-
-
-
-//namespace CoolProp{
-//
-//BasePolynomial::BasePolynomial(){
-//  this->POLYMATH_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<double> &coefficients, unsigned int n){
-//	if (coefficients.size() == n){
-//		return true;
-//	} else {
-//		throw ValueError(format("The number of coefficients %d does not match with %d. ",coefficients.size(),n));
-//	}
-//	return false;
-//}
-//bool BasePolynomial::checkCoefficients(std::vector< std::vector<double> > const& coefficients, unsigned int rows, unsigned int columns){
-//	if (coefficients.size() == rows){
-//		bool result = true;
-//		for(unsigned int i=0; i<rows; i++) {
-//			result = result && checkCoefficients(coefficients[i],columns);
-//		}
-//		return result;
-//	} else {
-//		throw ValueError(format("The number of rows %d does not match with %d. ",coefficients.size(),rows));
-//	}
-//	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<double> BasePolynomial::integrateCoeffs(std::vector<double> const& coefficients){
-//	std::vector<double> newCoefficients;
-//	unsigned int sizeX = coefficients.size();
-//	if (sizeX<1) throw ValueError(format("You have to provide coefficients, a vector length of %d is not a valid. ",sizeX));
-//	// pushing a zero elevates the order by 1
-//	newCoefficients.push_back(0.0);
-//	for(unsigned int i=0; i<coefficients.size(); i++) {
-//		newCoefficients.push_back(coefficients[i]/(i+1.));
-//	}
-//	return newCoefficients;
-//}
-//std::vector< std::vector<double> > BasePolynomial::integrateCoeffs(std::vector< std::vector<double> > const& coefficients, bool axis){
-//	std::vector< std::vector<double> > 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<double> > tmpCoefficients;
-//		tmpCoefficients = transpose(coefficients);
-//		unsigned int sizeY = tmpCoefficients.size();
-//		for(unsigned int i=0; i<sizeY; i++) {
-//			newCoefficients.push_back(integrateCoeffs(tmpCoefficients[i]));
-//		}
-//		return transpose(newCoefficients);
-//	} else if (axis==false){
-//		for(unsigned int i=0; i<sizeX; i++) {
-//			newCoefficients.push_back(integrateCoeffs(coefficients[i]));
-//		}
-//		return newCoefficients;
-//	} else {
-//		throw ValueError(format("You can only use x-axis (0) and y-axis (1) for integration. %d is not a valid input. ",axis));
-//	}
-//	return newCoefficients;
-//}
-//
-//
-///** Deriving coefficients for polynomials is done by multiplying the
-// *  original coefficients with i and lowering the order by 1.
-// */
-//std::vector<double> BasePolynomial::deriveCoeffs(std::vector<double> const& coefficients){
-//	std::vector<double> newCoefficients;
-//	unsigned int sizeX = coefficients.size();
-//	if (sizeX<1) throw ValueError(format("You have to provide coefficients, a vector length of %d is not a valid. ",sizeX));
-//	// skipping the first element lowers the order
-//	for(unsigned int i=1; i<coefficients.size(); i++) {
-//		newCoefficients.push_back(coefficients[i]*i);
-//	}
-//	return newCoefficients;
-//}
-//std::vector< std::vector<double> > BasePolynomial::deriveCoeffs(const std::vector< std::vector<double> > &coefficients, unsigned int axis){
-//	std::vector< std::vector<double> > 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<double> > tmpCoefficients;
-//		tmpCoefficients = transpose(coefficients);
-//		unsigned int sizeY = tmpCoefficients.size();
-//		for(unsigned int i=0; i<sizeY; i++) {
-//			newCoefficients.push_back(deriveCoeffs(tmpCoefficients[i]));
-//		}
-//		return transpose(newCoefficients);
-//	} else if (axis==1){
-//		for(unsigned int i=0; i<sizeX; i++) {
-//			newCoefficients.push_back(deriveCoeffs(coefficients[i]));
-//		}
-//		return newCoefficients;
-//	} else {
-//		throw ValueError(format("You can only use x-axis (0) and y-axis (1) for derivation. %d is not a valid input. ",axis));
-//	}
-//	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
-// *  the implementations, but mark them as deprecated.
-// *  The old functions are good for debugging since the
-// *  structure is easier to read than the backward Horner-scheme
-// *  or the recursive Horner-scheme.
-// */
-//
-///// Simple polynomial function generator. <- Deprecated due to poor performance, use Horner-scheme instead
-///** Base function to produce n-th order polynomials
-// *  based on the length of the coefficient vector.
-// *  Starts with only the first coefficient at x^0. */
-//double BasePolynomial::simplePolynomial(std::vector<double> const& coefficients, double x){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running simplePolynomial(std::vector, " << x << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//	double result = 0.;
-//	for(unsigned int i=0; i<coefficients.size();i++) {
-//		result += coefficients[i] * pow(x,(int)i);
-//	}
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//double BasePolynomial::simplePolynomial(std::vector<std::vector<double> > const& coefficients, double x, double y){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running simplePolynomial(std::vector, " << x << ", " << y << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//	double result = 0;
-//	for(unsigned int i=0; i<coefficients.size();i++) {
-//		result += pow(x,(int)i) * simplePolynomial(coefficients[i], y);
-//	}
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-//
-///// Simple integrated polynomial function generator.
-///** Base function to produce integrals of n-th order
-// *  polynomials based on the length of the coefficient
-// *  vector.
-// *  Starts with only the first coefficient at x^0 */
-/////Indefinite integral in x-direction
-//double BasePolynomial::simplePolynomialInt(std::vector<double> const& coefficients, double x){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running simplePolynomialInt(std::vector, " << x << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//	double result = 0.;
-//	for(unsigned int i=0; i<coefficients.size();i++) {
-//		result += 1./(i+1.) * coefficients[i] * pow(x,(int)(i+1.));
-//	}
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-/////Indefinite integral in y-direction only
-//double BasePolynomial::simplePolynomialInt(std::vector<std::vector<double> > const& coefficients, double x, double y){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running simplePolynomialInt(std::vector, " << x << ", " << y << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//	double result = 0.;
-//	for(unsigned int i=0; i<coefficients.size();i++) {
-//		result += pow(x,(int)i) * simplePolynomialInt(coefficients[i], y);
-//	}
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-//
-///// Simple integrated polynomial function generator divided by independent variable.
-///** Base function to produce integrals of n-th order
-// *  polynomials based on the length of the coefficient
-// *  vector.
-// *  Starts with only the first coefficient at x^0 */
-//double BasePolynomial::simpleFracInt(std::vector<double> const& coefficients, double x){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running      simpleFracInt(std::vector, " << x << "): ";
-//	}
-//	double result = coefficients[0] * log(x);
-//	if (coefficients.size() > 1) {
-//		for (unsigned int i=1; i<coefficients.size(); i++){
-//			result += 1./(i) * coefficients[i] * pow(x,(int)(i));
-//		}
-//	}
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-//double BasePolynomial::simpleFracInt(std::vector< std::vector<double> > const& coefficients, double x, double y){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running      simpleFracInt(std::vector, " << x << ", " << y << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//	double result = 0;
-//	for (unsigned int i=0; i<coefficients.size(); i++){
-//		result += pow(x,(int)i) * polyfracint(coefficients[i],y);
-//	}
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-//
-///** Simple integrated centred(!) polynomial function generator divided by independent variable.
-// *  We need to rewrite some of the functions in order to
-// *  use central fit. Having a central temperature xbase
-// *  allows for a better fit, but requires a different
-// *  formulation of the fracInt function group. Other
-// *  functions are not affected.
-// *  Starts with only the first coefficient at x^0 */
-//
-/////Helper functions to calculate binomial coefficients: http://rosettacode.org/wiki/Evaluate_binomial_coefficients#C.2B.2B
-////double BasePolynomial::factorial(double nValue){
-////   double result = nValue;
-////   double result_next;
-////   double pc = nValue;
-////   do {
-////	   result_next = result*(pc-1);
-////	   result = result_next;
-////	   pc--;
-////   } while(pc>2);
-////   nValue = result;
-////   return nValue;
-////}
-////double BasePolynomial::factorial(double nValue){
-////	if (nValue == 0) return (1);
-////	else return (nValue * factorial(nValue - 1));
-////}
-//double BasePolynomial::factorial(double nValue){
-//    double value = 1;
-//    for(int i = 2; i <= nValue; i++){
-//        value = value * i;
-//    }
-//    return value;
-//}
-//
-//double BasePolynomial::binom(double nValue, double nValue2){
-//   double result;
-//   if(nValue2 == 1) return nValue;
-//   result = (factorial(nValue)) / (factorial(nValue2)*factorial((nValue - nValue2)));
-//   nValue2 = result;
-//   return nValue2;
-//}
-//
-/////Helper functions to calculate the D vector:
-//std::vector<double> BasePolynomial::fracIntCentralDvector(int m, double x, double xbase){
-//	std::vector<double> 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<m; j++){ // loop through row
-//		tmp = pow(-1.0,j) * log(x) * pow(xbase,(int)j);
-//		for(int k=0; k<j; k++) { // internal loop for every entry
-//			tmp += binom(j,k) * pow(-1.0,k) / (j-k) * pow(x,j-k) * pow(xbase,k);
-//		}
-//		D.push_back(tmp);
-//	}
-//	return D;
-//}
-//
-/////Indefinite integral of a centred polynomial divided by its independent variable
-//double BasePolynomial::fracIntCentral(std::vector<double> const& coefficients, double x, double xbase){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running    fracIntCentral(std::vector, " << x << ", " << xbase << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//	int m = coefficients.size();
-//	std::vector<double> D = fracIntCentralDvector(m, x, xbase);
-//	double result = 0;
-//	for(int j=0; j<m; j++) {
-//		result += coefficients[j] * D[j];
-//	}
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-//
-///// Horner function generator implementations
-///** Represent polynomials according to Horner's scheme.
-// *  This avoids unnecessary multiplication and thus
-// *  speeds up calculation.
-// */
-//double BasePolynomial::baseHorner(std::vector<double> const& coefficients, double x){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running       baseHorner(std::vector, " << x << "): ";
-//	}
-//	double result = 0;
-//	for(int i=coefficients.size()-1; i>=0; i--) {
-//		result *= x;
-//		result += coefficients[i];
-//	}
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-//double BasePolynomial::baseHorner(std::vector< std::vector<double> > const& coefficients, double x, double y){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running       baseHorner(std::vector, " << x << ", " << y << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//	double result = 0;
-//	for(int i=coefficients.size()-1; i>=0; i--) {
-//		result *= x;
-//		result += baseHorner(coefficients[i], y);
-//	}
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-/////Indefinite integral in x-direction
-//double BasePolynomial::baseHornerInt(std::vector<double> const& coefficients, double x){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running       baseHornerInt(std::vector, " << x << "): ";
-//	}
-//	double result = 0;
-//	for(int i=coefficients.size()-1; i>=0; i--) {
-//		result *= x;
-//		result += coefficients[i]/(i+1.);
-//	}
-//	result = result * x;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-/////Indefinite integral in y-direction only
-//double BasePolynomial::baseHornerInt(std::vector<std::vector<double> > const& coefficients, double x, double y){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running       baseHornerInt(std::vector, " << x << ", " << y << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//	double result = 0;
-//	for(int i=coefficients.size()-1; i>=0; i--) {
-//		result *= x;
-//		result += baseHornerInt(coefficients[i], y);
-//	}
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-/////Indefinite integral in x-direction of a polynomial divided by its independent variable
-//double BasePolynomial::baseHornerFracInt(std::vector<double> const& coefficients, double x){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running      baseHornerFra(std::vector, " << x << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//	double result = 0;
-//	if (coefficients.size() > 1) {
-//		for(int i=coefficients.size()-1; i>=1; i--) {
-//			result *= x;
-//			result += coefficients[i]/(i);
-//		}
-//		result *= x;
-//	}
-//	result += coefficients[0] * log(x);
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-/////Indefinite integral in y-direction of a polynomial divided by its 2nd independent variable
-//double BasePolynomial::baseHornerFracInt(std::vector<std::vector<double> > const& coefficients, double x, double y){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running      baseHornerFra(std::vector, " << x << ", " << y << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//
-//	double result = 0;
-//	for(int i=coefficients.size()-1; i>=0; i--) {
-//		result *= x;
-//		result += baseHornerFracInt(coefficients[i], y);
-//	}
-//
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-//
-///** Alternatives
-// *  Simple functions that heavily rely on other parts of this file.
-// *  We still need to check which combinations yield the best
-// *  performance.
-// */
-/////Derivative in x-direction
-//double BasePolynomial::deriveIn2Steps(std::vector<double> const& coefficients, double x){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running   deriveIn2Steps(std::vector, " << x << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//	double result =  polyval(deriveCoeffs(coefficients),x);
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-/////Derivative in terms of x(axis=true) or y(axis=false).
-//double BasePolynomial::deriveIn2Steps(std::vector< std::vector<double> > const& coefficients, double x, double y, bool axis){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running   deriveIn2Steps(std::vector, " << x << ", " << y << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//	double result = polyval(deriveCoeffs(coefficients,axis),x,y);
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-/////Indefinite integral in x-direction
-//double BasePolynomial::integrateIn2Steps(std::vector<double> const& coefficients, double x){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running   integrateIn2Steps(std::vector, " << x << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//	double result =  polyval(integrateCoeffs(coefficients),x);
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-/////Indefinite integral in terms of x(axis=true) or y(axis=false).
-//double BasePolynomial::integrateIn2Steps(std::vector< std::vector<double> > const& coefficients, double x, double y, bool axis){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running   integrateIn2Steps(std::vector, " << x << ", " << y << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//	double result = polyval(integrateCoeffs(coefficients,axis),x,y);
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-/////Indefinite integral in x-direction of a polynomial divided by its independent variable
-//double BasePolynomial::fracIntIn2Steps(std::vector<double> const& coefficients, double x){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running    fracIntIn2Steps(std::vector, " << x << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//	double result = coefficients[0] * log(x);
-//	if (coefficients.size() > 1) {
-//		std::vector<double> newCoeffs(coefficients.begin() + 1, coefficients.end());
-//		result += polyint(newCoeffs,x);
-//	}
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-/////Indefinite integral in y-direction of a polynomial divided by its 2nd independent variable
-//double BasePolynomial::fracIntIn2Steps(std::vector<std::vector<double> > const& coefficients, double x, double y){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running    fracIntIn2Steps(std::vector, " << x << ", " << y << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//	std::vector<double> newCoeffs;
-//	for (unsigned int i=0; i<coefficients.size(); i++){
-//		newCoeffs.push_back(polyfracint(coefficients[i],y));
-//	}
-//	double result = polyval(newCoeffs,x);
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-/////Indefinite integral in y-direction of a centred polynomial divided by its 2nd independent variable
-//double BasePolynomial::fracIntCentral2Steps(std::vector<std::vector<double> > const& coefficients, double x, double y, double ybase){
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << "Running    fracIntCentral2Steps(std::vector, " << x << ", " << y << ", " << ybase << "): ";
-//	}
-//	bool db = this->POLYMATH_DEBUG;
-//	this->POLYMATH_DEBUG = false;
-//	std::vector<double> newCoeffs;
-//	for (unsigned int i=0; i<coefficients.size(); i++){
-//		newCoeffs.push_back(fracIntCentral(coefficients[i], y, ybase));
-//	}
-//	double result = polyval(newCoeffs,x);
-//	this->POLYMATH_DEBUG = db;
-//	if (this->POLYMATH_DEBUG) {
-//		std::cout << result << std::endl;
-//	}
-//	return result;
-//}
-//
-//
-//
-//
-///** Implements the function wrapper interface and can be
-// *  used by the solvers. This is only an example and you should
-// *  use local redefinitions of the class.
-// *  TODO: Make multidimensional
-// */
-//PolyResidual::PolyResidual(){
-//	this->dim = -1;
-//}
-//
-//PolyResidual::PolyResidual(const std::vector<double> &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<double> > &coefficients, double x, double z){
-//	this->output = z;
-//	this->firstDim = x;
-//	this->coefficients = coefficients;
-//	this->dim = i2D;
-//}
-//
-//double PolyResidual::call(double x){
-//	double polyRes = -1;
-//	switch (this->dim) {
-//	case i1D:
-//		polyRes = this->poly.polyval(this->coefficients[0], x);
-//		break;
-//	case i2D:
-//		polyRes = this->poly.polyval(this->coefficients, this->firstDim, x);
-//		break;
-//	default:
-//		throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D.");
-//	}
-//	return polyRes - this->output;
-//}
-//
-//double PolyResidual::deriv(double x){
-//	double polyRes = -1;
-//	switch (this->dim) {
-//	case i1D:
-//		polyRes = this->poly.polyder(this->coefficients[0], x);
-//		break;
-//	case i2D:
-//		polyRes = this->poly.polyder(this->coefficients, this->firstDim, x);
-//		break;
-//	default:
-//		throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D.");
-//	}
-//	return polyRes;
-//}
-//
-//double PolyIntResidual::call(double x){
-//	double polyRes = -1;
-//	switch (this->dim) {
-//	case i1D:
-//		polyRes = this->poly.polyint(this->coefficients[0], x);
-//		break;
-//	case i2D:
-//		polyRes = this->poly.polyint(this->coefficients, this->firstDim, x);
-//		break;
-//	default:
-//		throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D.");
-//	}
-//	return polyRes - this->output;
-//}
-//
-//double PolyIntResidual::deriv(double x){
-//	double polyRes = -1;
-//	switch (this->dim) {
-//	case i1D:
-//		polyRes = this->poly.polyval(this->coefficients[0], x);
-//		break;
-//	case i2D:
-//		polyRes = this->poly.polyval(this->coefficients, this->firstDim, x);
-//		break;
-//	default:
-//		throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D.");
-//	}
-//	return polyRes;
-//}
-//
-//double PolyFracIntResidual::call(double x){
-//	double polyRes = -1;
-//	switch (this->dim) {
-//	case i1D:
-//		polyRes = this->poly.polyfracint(this->coefficients[0], x);
-//		break;
-//	case i2D:
-//		polyRes = this->poly.polyfracint(this->coefficients, this->firstDim, x);
-//		break;
-//	default:
-//		throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D.");
-//	}
-//	return polyRes - this->output;
-//}
-//
-//double PolyFracIntResidual::deriv(double x){
-//	double polyRes = -1;
-//	switch (this->dim) {
-//	case i1D:
-//		polyRes = this->poly.polyfracval(this->coefficients[0], x);
-//		break;
-//	case i2D:
-//		polyRes = this->poly.polyfracval(this->coefficients, this->firstDim, x);
-//		break;
-//	default:
-//		throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D.");
-//	}
-//	return polyRes;
-//}
-//
-//double PolyFracIntCentralResidual::call(double x){
-//	double polyRes = -1;
-//	switch (this->dim) {
-//	case i1D:
-//		polyRes = this->poly.polyfracintcentral(this->coefficients[0], x, this->baseVal);
-//		break;
-//	case i2D:
-//		polyRes = this->poly.polyfracintcentral(this->coefficients, this->firstDim, x, this->baseVal);
-//		break;
-//	default:
-//		throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D.");
-//	}
-//	return polyRes - this->output;
-//}
-//
-//double PolyFracIntCentralResidual::deriv(double x){
-//	throw CoolProp::NotImplementedError("Derivative of a polynomial frac int is not defined.");
-//}
-//
-//double PolyDerResidual::call(double x){
-//	double polyRes = -1;
-//	switch (this->dim) {
-//	case i1D:
-//		polyRes = this->poly.polyder(this->coefficients[0], x);
-//		break;
-//	case i2D:
-//		polyRes = this->poly.polyder(this->coefficients, this->firstDim, x);
-//		break;
-//	default:
-//		throw CoolProp::NotImplementedError("There are only 1D and 2D, a polynomial's live is not 3D.");
-//	}
-//	return polyRes - this->output;
-//}
-//
-//double PolyDerResidual::deriv(double x){
-//	throw CoolProp::NotImplementedError("2nd derivative of a polynomial is not defined.");
-//}
-//
-//
-//
-//
-///** Implements the same public functions as the BasePolynomial
-// *  but solves the polynomial for the given value
-// *  instead of evaluating it.
-// *  TODO: This class does not check for bijective
-// *        polynomials and is therefore a little
-// *        fragile.
-// */
-//PolynomialSolver::PolynomialSolver(){
-//  this->POLYMATH_DEBUG   = false;
-//  this->macheps = DBL_EPSILON;
-//  this->tol     = DBL_EPSILON*1e3;
-//  this->maxiter = 50;
-//}
-//
-///** Everything related to the normal polynomials goes in this
-// *  section, holds all the functions for solving polynomials.
-// */
-///// Solves a one-dimensional polynomial for the given coefficients
-///// @param coefficients vector containing the ordered coefficients
-///// @param y double value that represents the current input
-//double PolynomialSolver::polyval(const std::vector<double> &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<double> > &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<double> &coefficients, double y){
-//	PolyIntResidual residual = PolyIntResidual(coefficients, y);
-//	return this->solve(residual);
-//}
-//
-///// Solves the indefinite integral of a two-dimensional polynomial along the 2nd axis (y)
-///// @param coefficients vector containing the ordered coefficients
-///// @param x double value that represents the current input in the 1st dimension
-///// @param z double value that represents the current output
-//double PolynomialSolver::polyint(const std::vector< std::vector<double> > &coefficients, double x, double z){
-//	PolyIntResidual residual = PolyIntResidual(coefficients, x, z);
-//	return this->solve(residual);
-//}
-//
-//
-///** Everything related to the derived polynomials goes in this
-// *  section, holds all the functions for solving polynomials.
-// */
-///// Solves the derivative of a one-dimensional polynomial
-///// @param coefficients vector containing the ordered coefficients
-///// @param y double value that represents the current output
-//double PolynomialSolver::polyder(const std::vector<double> &coefficients, double y){
-//	PolyDerResidual residual = PolyDerResidual(coefficients, y);
-//	return this->solve(residual);
-//}
-//
-///// Solves the derivative of a two-dimensional polynomial along the 2nd axis (y)
-///// @param coefficients vector containing the ordered coefficients
-///// @param x double value that represents the current input in the 1st dimension
-///// @param z double value that represents the current output
-//double PolynomialSolver::polyder(const std::vector< std::vector<double> > &coefficients, double x, double z){
-//	PolyDerResidual residual = PolyDerResidual(coefficients, x, z);
-//	return this->solve(residual);
-//}
-//
-//
-///** Everything related to the polynomials divided by one variable goes in this
-// *  section, holds all the functions for solving polynomials.
-// */
-///// Solves the indefinite integral of a one-dimensional polynomial divided by its independent variable
-///// @param coefficients vector containing the ordered coefficients
-///// @param y double value that represents the current output
-//double PolynomialSolver::polyfracval(const std::vector<double> &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<double> > &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<double> &coefficients, double y){
-//	PolyFracIntResidual residual = PolyFracIntResidual(coefficients, y);
-//	return this->solve(residual);
-//}
-//
-///// Solves the indefinite integral of a two-dimensional polynomial divided by its 2nd independent variable
-///// @param coefficients vector containing the ordered coefficients
-///// @param x double value that represents the current input in the 1st dimension
-///// @param z double value that represents the current output
-//double PolynomialSolver::polyfracint(const std::vector< std::vector<double> > &coefficients, double x, double z){
-//	PolyFracIntResidual residual = PolyFracIntResidual(coefficients, x, z);
-//	return this->solve(residual);
-//}
-//
-///// Solves the indefinite integral of a centred one-dimensional polynomial divided by its independent variable
-///// @param coefficients vector containing the ordered coefficients
-///// @param y double value that represents the current output
-///// @param xbase central x-value for fitted function
-//double PolynomialSolver::polyfracintcentral(const std::vector<double> &coefficients, double y, double xbase){
-//	PolyFracIntCentralResidual residual = PolyFracIntCentralResidual(coefficients, y, xbase);
-//	return this->solve(residual);
-//}
-//
-///// Solves the indefinite integral of a centred two-dimensional polynomial divided by its 2nd independent variable
-///// @param coefficients vector containing the ordered coefficients
-///// @param x double value that represents the current input in the 1st dimension
-///// @param z double value that represents the current output
-///// @param ybase central y-value for fitted function
-//double PolynomialSolver::polyfracintcentral(const std::vector< std::vector<double> > &coefficients, double x, double z, double ybase){
-//	PolyFracIntCentralResidual residual = PolyFracIntCentralResidual(coefficients, x, z, ybase);
-//	return this->solve(residual);
-//}
-//
-//
-///** Everything related to the derived polynomials divided by one variable goes in this
-// *  section, holds all the functions for solving polynomials.
-// */
-///// Solves the derivative of a one-dimensional polynomial divided by its independent variable
-///// @param coefficients vector containing the ordered coefficients
-///// @param y double value that represents the current output
-//double PolynomialSolver::polyfracder(const std::vector<double> &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<double> > &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<double> &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<double> > &coefficients, double x, double z, double ybase){
-//	throw CoolProp::NotImplementedError("This solver has not been implemented, yet."); // TODO: Implement function
-//}
-//
-//
-///** Set the solvers and updates either the guess values or the
-// *  boundaries for the variable to solve for.
-// */
-///// Sets the guess value for the Newton solver and enables it.
-///// @param guess double value that represents the guess value
-//void PolynomialSolver::setGuess(double guess){
-//	this->uses  = iNewton;
-//	this->guess = guess;
-//	this->min   = -1;
-//	this->max   = -1;
-//}
-///// Sets the limits for the Brent solver and enables it.
-///// @param min double value that represents the lower boundary
-///// @param max double value that represents the upper boundary
-//void PolynomialSolver::setLimits(double min, double max){
-//	this->uses  = iBrent;
-//	this->guess = -1;
-//	this->min   = min;
-//	this->max   = max;
-//}
-//
-///// Solves the equations based on previously defined parameters.
-///// @param min double value that represents the lower boundary
-///// @param max double value that represents the upper boundary
-//double PolynomialSolver::solve(PolyResidual &res){
-//	std::string errstring;
-//	double result = -1.0;
-//	switch (this->uses) {
-//	case iNewton: ///< Newton solver with derivative and guess value
-//		if (res.is2D()) {
-//			throw CoolProp::NotImplementedError("The Newton solver is not suitable for 2D polynomials, yet.");
-//		}
-//		result = Newton(res, this->guess, this->tol, this->maxiter, errstring);
-//		break;
-//
-//	case iBrent: ///< Brent solver with bounds
-//		result = Brent(res, this->min, this->max, this->macheps, this->tol, this->maxiter, errstring);
-//		break;
-//
-//	default:
-//		throw CoolProp::NotImplementedError("This solver has not been implemented or you forgot to select a solver...");
-//	}
-//	return result;
-//}
-//
-//
-///** Here we define the functions that should be to evaluate exponential
-// *  functions. Not really polynomials, I know...
-// */
-//
-//BaseExponential::BaseExponential(){
-//  this->POLYMATH_DEBUG = false;
-////  this->poly = new BaseExponential();
-//}
-////
-////BaseExponential::~BaseExponential(){
-////  delete this->poly;
-////}
-//
-///// Evaluates an exponential function for the given coefficients
-///// @param coefficients vector containing the ordered coefficients
-///// @param x double value that represents the current input
-///// @param n int value that determines the kind of exponential function
-//double BaseExponential::expval(const std::vector<double> &coefficients, double x, int n){
-//	double result = 0.;
-//	if (n==1) {
-//		this->poly.checkCoefficients(coefficients,3);
-//		result = exp(coefficients[0]/(x+coefficients[1]) - coefficients[2]);
-//	} else if (n==2) {
-//		result = exp(this->poly.polyval(coefficients, x));
-//	} else {
-//		throw ValueError(format("There is no function defined for this input (%d). ",n));
-//	}
-//	return result;
-//}
-//
-///// Evaluates an exponential function for the given coefficients
-///// @param coefficients vector containing the ordered coefficients
-///// @param x double value that represents the current input in the 1st dimension
-///// @param y double value that represents the current input in the 2nd dimension
-///// @param n int value that determines the kind of exponential function
-//double BaseExponential::expval(const std::vector< std::vector<double> > &coefficients, double x, double y, int n){
-//	double result = 0.;
-//	if (n==2) {
-//		result = exp(this->poly.polyval(coefficients, x, y));
-//	} else {
-//		throw ValueError(format("There is no function defined for this input (%d). ",n));
-//	}
-//	return result;
-//}
-//
-//
-//}
-//
-//
-//#ifdef ENABLE_CATCH
-//#include <math.h>
-//#include "catch.hpp"
-//
-//class PolynomialConsistencyFixture {
-//public:
-//	CoolProp::BasePolynomial poly;
-//	CoolProp::PolynomialSolver solver;
-////	enum dims {i1D, i2D};
-////	double firstDim;
-////	int dim;
-////	std::vector< std::vector<double> > coefficients;
-////
-////    void setInputs(const std::vector<double> &coefficients){
-////    	this->firstDim = 0;
-////    	this->coefficients.clear();
-////    	this->coefficients.push_back(coefficients);
-////    	this->dim = i1D;
-////    }
-////
-////    void setInputs(const std::vector< std::vector<double> > &coefficients, double x){
-////    	this->firstDim = x;
-////    	this->coefficients = coefficients;
-////    	this->dim = i2D;
-////    }
-//};
-//
-//
-//TEST_CASE("Internal consistency checks with PolyMath objects","[PolyMath]")
-//{
-//	CoolProp::BasePolynomial poly;
-//	CoolProp::PolynomialSolver solver;
-//
-//	/// Test case for "SylthermXLT" by "Dow Chemicals"
-//	std::vector<double> cHeat;
-//	cHeat.clear();
-//	cHeat.push_back(+1.1562261074E+03);
-//	cHeat.push_back(+2.0994549103E+00);
-//	cHeat.push_back(+7.7175381057E-07);
-//	cHeat.push_back(-3.7008444051E-20);
-//
-//	double deltaT = 0.1;
-//	double Tmin   = 273.15- 50;
-//	double Tmax   = 273.15+250;
-//	double Tinc   = 15;
-//
-//	double val1,val2,val3,val4;
-//
-//	SECTION("DerFromVal1D") {
-//		for (double T = Tmin; T<Tmax; T+=Tinc) {
-//			val1 = poly.polyval(cHeat, T-deltaT);
-//			val2 = poly.polyval(cHeat, T+deltaT);
-//			val3 = (val2-val1)/2/deltaT;
-//			val4 = poly.polyder(cHeat, T);
-//			CAPTURE(T);
-//			CAPTURE(val3);
-//			CAPTURE(val4);
-//			CHECK( (1.0-fabs(val4/val3)) < 1e-1);
-//		}
-//	}
-//	SECTION("ValFromInt1D") {
-//		for (double T = Tmin; T<Tmax; T+=Tinc) {
-//			val1 = poly.polyint(cHeat, T-deltaT);
-//			val2 = poly.polyint(cHeat, T+deltaT);
-//			val3 = (val2-val1)/2/deltaT;
-//			val4 = poly.polyval(cHeat, T);
-//			CAPTURE(T);
-//			CAPTURE(val3);
-//			CAPTURE(val4);
-//			CHECK( (1.0-fabs(val4/val3)) < 1e-1);
-//		}
-//	}
-//
-//	SECTION("Solve1DNewton") {
-//		for (double T = Tmin; T<Tmax; T+=Tinc) {
-//			val1 = poly.polyval(cHeat, T);
-//			solver.setGuess(T+100);
-//			val2 = solver.polyval(cHeat, val1);
-//			CAPTURE(T);
-//			CAPTURE(val1);
-//			CAPTURE(val2);
-//			CHECK(fabs(T-val2) < 1e-1);
-//
-//			val1 = poly.polyint(cHeat, T);
-//			solver.setGuess(T+100);
-//			val2 = solver.polyint(cHeat, val1);
-//			CAPTURE(T);
-//			CAPTURE(val1);
-//			CAPTURE(val2);
-//			CHECK(fabs(T-val2) < 1e-1);
-//
-////			val1 = poly.polyder(cHeat, T);
-////			solver.setGuess(T+100);
-////			val2 = solver.polyder(cHeat, val1);
-////			CAPTURE(T);
-////			CAPTURE(val1);
-////			CAPTURE(val2);
-////			CHECK(fabs(T-val2) < 1e-1);
-////
-////			val1 = poly.polyfracint(cHeat, T);
-////			solver.setGuess(T+100);
-////			val2 = solver.polyfracint(cHeat, val1);
-////			CAPTURE(T);
-////			CAPTURE(val1);
-////			CAPTURE(val2);
-////			CHECK(fabs(T-val2) < 1e-1);
-//		}
-//	}
-//	SECTION("Solve1DBrent") {
-//		for (double T = Tmin; T<Tmax; T+=Tinc) {
-//			val1 = poly.polyval(cHeat, T);
-//			solver.setLimits(T-300,T+300);
-//			val2 = solver.polyval(cHeat, val1);
-//			CAPTURE(T);
-//			CAPTURE(val1);
-//			CAPTURE(val2);
-//			CHECK(fabs(T-val2) < 1e-1);
-//
-//			val1 = poly.polyint(cHeat, T);
-//			solver.setLimits(T-300,T+300);
-//			val2 = solver.polyint(cHeat, val1);
-//			CAPTURE(T);
-//			CAPTURE(val1);
-//			CAPTURE(val2);
-//			CHECK(fabs(T-val2) < 1e-1);
-//
-//			val1 = poly.polyder(cHeat, T);
-//			solver.setLimits(T-300,T+300);
-//			val2 = solver.polyder(cHeat, val1);
-//			CAPTURE(T);
-//			CAPTURE(val1);
-//			CAPTURE(val2);
-//			CHECK(fabs(T-val2) < 1e-1);
-//
-//			val1 = poly.polyfracint(cHeat, T);
-//			solver.setLimits(T-100,T+100);
-//			val2 = solver.polyfracint(cHeat, val1);
-//			CAPTURE(T);
-//			CAPTURE(val1);
-//			CAPTURE(val2);
-//			CHECK(fabs(T-val2) < 1e-1);
-//
-//			val1 = poly.polyfracintcentral(cHeat, T, 250.0);
-//			solver.setLimits(T-100,T+100);
-//			val2 = solver.polyfracintcentral(cHeat, val1, 250.0);
-//			CAPTURE(T);
-//			CAPTURE(val1);
-//			CAPTURE(val2);
-//			CHECK(fabs(T-val2) < 1e-1);
-//
-//		}
-//	}
-//
-//	/// Test case for 2D
-//	double xDim1 = 0.3;
-//	std::vector< std::vector<double> > cHeat2D;
-//	cHeat2D.clear();
-//	cHeat2D.push_back(cHeat);
-//	cHeat2D.push_back(cHeat);
-//	cHeat2D.push_back(cHeat);
-//
-//	//setInputs(cHeat2D, 0.3);
-//
-//	SECTION("DerFromVal2D") {
-//		for (double T = Tmin; T<Tmax; T+=Tinc) {
-//			val1 = poly.polyval(cHeat2D, xDim1, T-deltaT);
-//			val2 = poly.polyval(cHeat2D, xDim1, T+deltaT);
-//			val3 = (val2-val1)/2/deltaT;
-//			val4 = poly.polyder(cHeat2D, xDim1, T);
-//			CAPTURE(T);
-//			CAPTURE(val3);
-//			CAPTURE(val4);
-//			CHECK( (1.0-fabs(val4/val3)) < 1e-1);
-//		}
-//	}
-//
-//	SECTION("ValFromInt2D") {
-//		for (double T = Tmin; T<Tmax; T+=Tinc) {
-//			val1 = poly.polyint(cHeat2D, xDim1, T-deltaT);
-//			val2 = poly.polyint(cHeat2D, xDim1, T+deltaT);
-//			val3 = (val2-val1)/2/deltaT;
-//			val4 = poly.polyval(cHeat2D, xDim1, T);
-//			CAPTURE(T);
-//			CAPTURE(val3);
-//			CAPTURE(val4);
-//			CHECK( (1.0-fabs(val4/val3)) < 1e-1);
-//		}
-//	}
-//
-////	SECTION("Solve2DNewton") {
-////		for (double T = Tmin; T<Tmax; T+=Tinc) {
-////			val1 = poly.polyval(cHeat2D, xDim1, T);
-////			solver.setGuess(T+100);
-////			val2 = solver.polyval(cHeat2D, xDim1, val1);
-////			CAPTURE(T);
-////			CAPTURE(val1);
-////			CAPTURE(val2);
-////			CHECK(fabs(T-val2) < 1e-1);
-////		}
-////	}
-//	SECTION("Solve2DBrent") {
-//		for (double T = Tmin; T<Tmax; T+=Tinc) {
-//			val1 = poly.polyval(cHeat2D, xDim1, T);
-//			solver.setLimits(T-300,T+300);
-//			val2 = solver.polyval(cHeat2D, xDim1, val1);
-//			CAPTURE(T);
-//			CAPTURE(val1);
-//			CAPTURE(val2);
-//			CHECK(fabs(T-val2) < 1e-1);
-//
-//			val1 = poly.polyint(cHeat2D, xDim1, T);
-//			solver.setLimits(T-300,T+300);
-//			val2 = solver.polyint(cHeat2D, xDim1, val1);
-//			CAPTURE(T);
-//			CAPTURE(val1);
-//			CAPTURE(val2);
-//			CHECK(fabs(T-val2) < 1e-1);
-//
-//			val1 = poly.polyder(cHeat2D, xDim1, T);
-//			solver.setLimits(T-300,T+300);
-//			val2 = solver.polyder(cHeat2D, xDim1, val1);
-//			CAPTURE(T);
-//			CAPTURE(val1);
-//			CAPTURE(val2);
-//			CHECK(fabs(T-val2) < 1e-1);
-//
-//			val1 = poly.polyfracint(cHeat2D, xDim1, T);
-//			solver.setLimits(T-100,T+100);
-//			val2 = solver.polyfracint(cHeat2D, xDim1, val1);
-//			CAPTURE(T);
-//			CAPTURE(val1);
-//			CAPTURE(val2);
-//			CHECK(fabs(T-val2) < 1e-1);
-//
-//			val1 = poly.polyfracintcentral(cHeat2D, xDim1, T, 250);
-//			solver.setLimits(T-100,T+100);
-//			val2 = solver.polyfracintcentral(cHeat2D, xDim1, val1, 250);
-//			CAPTURE(T);
-//			CAPTURE(val1);
-//			CAPTURE(val2);
-//			CHECK(fabs(T-val2) < 1e-1);
-//		}
-//	}
-//
-//}
-//
-////TEST_CASE_METHOD(PolynomialConsistencyFixture,"Internal consistency checks","[PolyMath]")
-////{
-////	/// Test case for "SylthermXLT" by "Dow Chemicals"
-////	std::vector<double> cHeat;
-////	cHeat.clear();
-////	cHeat.push_back(+1.1562261074E+03);
-////	cHeat.push_back(+2.0994549103E+00);
-////	cHeat.push_back(+7.7175381057E-07);
-////	cHeat.push_back(-3.7008444051E-20);
-////
-////	//setInputs(cHeat);
-////	double deltaT = 0.1;
-////	double val1,val2,val3,val4;
-////
-////    SECTION("DerFromVal1D") {
-////		for (double T = 273.15-50; T<273.15+250; T+=15) {
-////			val1 = this->poly.polyval(cHeat, T-deltaT);
-////			val2 = this->poly.polyval(cHeat, T+deltaT);
-////			val3 = (val2-val1)/2/deltaT;
-////			val4 = this->poly.polyder(cHeat, T);
-////			CAPTURE(T);
-////			CAPTURE(val3);
-////			CAPTURE(val4);
-////			CHECK( (1.0-fabs(val4/val3)) < 1e-1);
-////		}
-////    }
-////
-////    SECTION("ValFromInt1D") {
-////		for (double T = 273.15-50; T<273.15+250; T+=15) {
-////			val1 = this->poly.polyint(cHeat, T-deltaT);
-////			val2 = this->poly.polyint(cHeat, T+deltaT);
-////			val3 = (val2-val1)/2/deltaT;
-////			val4 = this->poly.polyval(cHeat, T);
-////			CAPTURE(T);
-////			CAPTURE(val3);
-////			CAPTURE(val4);
-////			CHECK( (1.0-fabs(val4/val3)) < 1e-1);
-////		}
-////	}
-////
-////    SECTION("Solve1DNewton") {
-////		for (double T = 273.15-50; T<273.15+250; T+=15) {
-////			val1 = this->poly.polyval(cHeat, T);
-////			this->solver.setGuess(T+100);
-////			val2 = this->solver.polyval(cHeat, val1);
-////			CAPTURE(T);
-////			CAPTURE(val1);
-////			CAPTURE(val2);
-////			CHECK(fabs(T-val2) < 1e-1);
-////		}
-////	}
-////    SECTION("Solve1DBrent") {
-////		for (double T = 273.15-50; T<273.15+250; T+=15) {
-////			val1 = this->poly.polyval(cHeat, T);
-////			this->solver.setLimits(T-300,T+300);
-////			val2 = this->solver.polyval(cHeat, val1);
-////			CAPTURE(T);
-////			CAPTURE(val1);
-////			CAPTURE(val2);
-////			CHECK(fabs(T-val2) < 1e-1);
-////		}
-////	}
-////
-////    /// Test case for 2D
-////    std::vector< std::vector<double> > cHeat2D;
-////	cHeat2D.clear();
-////	cHeat2D.push_back(cHeat);
-////	cHeat2D.push_back(cHeat);
-////	cHeat2D.push_back(cHeat);
-////
-////	//setInputs(cHeat2D, 0.3);
-////
-////	SECTION("DerFromVal2D") {
-////		for (double T = 273.15-50; T<273.15+250; T+=15) {
-////			val1 = this->poly.polyval(cHeat, T-deltaT);
-////			val2 = this->poly.polyval(cHeat, T+deltaT);
-////			val3 = (val2-val1)/2/deltaT;
-////			val4 = this->poly.polyder(cHeat, T);
-////			CAPTURE(T);
-////			CAPTURE(val3);
-////			CAPTURE(val4);
-////			CHECK( (1.0-fabs(val4/val3)) < 1e-1);
-////		}
-////	}
-////
-////	SECTION("ValFromInt2D") {
-////		for (double T = 273.15-50; T<273.15+250; T+=15) {
-////			val1 = this->poly.polyint(cHeat, T-deltaT);
-////			val2 = this->poly.polyint(cHeat, T+deltaT);
-////			val3 = (val2-val1)/2/deltaT;
-////			val4 = this->poly.polyval(cHeat, T);
-////			CAPTURE(T);
-////			CAPTURE(val3);
-////			CAPTURE(val4);
-////			CHECK( (1.0-fabs(val4/val3)) < 1e-1);
-////		}
-////	}
-////
-////    SECTION("Solve2DNewton") {
-////		for (double T = 273.15-50; T<273.15+250; T+=15) {
-////			val1 = this->poly.polyval(cHeat, T);
-////			this->solver.setGuess(T+100);
-////			val2 = this->solver.polyval(cHeat, val1);
-////			CAPTURE(T);
-////			CAPTURE(val1);
-////			CAPTURE(val2);
-////			CHECK(fabs(T-val2) < 1e-1);
-////		}
-////	}
-////    SECTION("Solve2DBrent") {
-////		for (double T = 273.15-50; T<273.15+250; T+=15) {
-////			val1 = this->poly.polyval(cHeat, T);
-////			this->solver.setLimits(T-300,T+300);
-////			val2 = this->solver.polyval(cHeat, val1);
-////			CAPTURE(T);
-////			CAPTURE(val1);
-////			CAPTURE(val2);
-////			CHECK(fabs(T-val2) < 1e-1);
-////		}
-////	}
-////
-////}
-////
-////TEST_CASE("Check against hard coded data","[PolyMath]")
-////{
-////    CHECK(fabs(HumidAir::f_factor(-60+273.15,101325)/(1.00708)-1) < 1e-3);
-////    CHECK(fabs(HumidAir::f_factor( 80+273.15,101325)/(1.00573)-1) < 1e-3);
-////    CHECK(fabs(HumidAir::f_factor(-60+273.15,10000e3)/(2.23918)-1) < 1e-3);
-////    CHECK(fabs(HumidAir::f_factor(300+273.15,10000e3)/(1.04804)-1) < 1e-3);
-////}
-//
-//
-//
-////int main() {
-////
-////	Catch::ConfigData &config = session.configData();
-////	config.testsOrTags.clear();
-////	config.testsOrTags.push_back("[fast]");
-////	session.useConfigData(config);
-////	return session.run();
-////
-////}
-//
-//#endif /* CATCH_ENABLED */
-//
-//
-////int main() {
-////
-////	std::vector<double> cHeat;
-////	cHeat.clear();
-////	cHeat.push_back(+1.1562261074E+03);
-////	cHeat.push_back(+2.0994549103E+00);
-////	cHeat.push_back(+7.7175381057E-07);
-////	cHeat.push_back(-3.7008444051E-20);
-////
-////	CoolProp::BasePolynomial base = CoolProp::BasePolynomial();
-////	CoolProp::PolynomialSolver solve = CoolProp::PolynomialSolver();
-////
-////	double T = 273.15+50;
-////
-////	double c = base.polyval(cHeat,T);
-////	printf("Should be  :      c = %3.3f \t J/kg/K    \n",1834.746);
-////	printf("From object:      c = %3.3f \t J/kg/K    \n",c);
-////
-////	T = 0.0;
-////	solve.setGuess(75+273.15);
-////	T = solve.polyval(cHeat,c);
-////	printf("Should be  :      T = %3.3f \t K    \n",273.15+50.0);
-////	printf("From object:      T = %3.3f \t K    \n",T);
-////
-////	T = 0.0;
-////	solve.setLimits(273.15+10,273.15+100);
-////	T = solve.polyval(cHeat,c);
-////	printf("Should be  :      T = %3.3f \t K    \n",273.15+50.0);
-////	printf("From object:      T = %3.3f \t K    \n",T);
-////
-////}
diff --git a/src/Tests/TestObjects.cpp b/src/Tests/TestObjects.cpp
index 050561e7..09f96d5f 100644
--- a/src/Tests/TestObjects.cpp
+++ b/src/Tests/TestObjects.cpp
@@ -173,7 +173,7 @@ CoolProp::IncompressibleFluid CoolPropTesting::incompressibleFluidObject(){
 
 	// After preparing the coefficients, we have to create the objects
 	CoolProp::IncompressibleFluid CH3OH;
-	CH3OH.setName("CH3OH");
+	CH3OH.setName("CH3OH-testing");
 	CH3OH.setDescription("Methanol solution");
 	CH3OH.setReference("SecCool software");
 	CH3OH.setTmax( 20 + 273.15);