CoolProp  6.6.0
An open-source fluid property and humid air property database
PolyMath.cpp
Go to the documentation of this file.
1 
2 #include "PolyMath.h"
3 
4 #include "Exceptions.h"
5 #include "MatrixMath.h"
6 
7 #include <vector>
8 #include <string>
9 //#include <sstream>
10 //#include <numeric>
11 #include <math.h>
12 #include <iostream>
13 
14 #include "Solvers.h"
15 
16 #include "unsupported/Eigen/Polynomials"
17 
18 namespace CoolProp {
19 
20 constexpr double CPPOLY_EPSILON = DBL_EPSILON * 100.0;
21 constexpr double CPPOLY_DELTA = CPPOLY_EPSILON * 10.0;
22 
24 
30 bool Polynomial2D::checkCoefficients(const Eigen::MatrixXd& coefficients, const unsigned int rows, const unsigned int columns) {
31  if (static_cast<size_t>(coefficients.rows()) == rows) {
32  if (static_cast<size_t>(coefficients.cols()) == columns) {
33  return true;
34  } else {
35  throw ValueError(format("%s (%d): The number of columns %d does not match with %d. ", __FILE__, __LINE__, coefficients.cols(), columns));
36  }
37  } else {
38  throw ValueError(format("%s (%d): The number of rows %d does not match with %d. ", __FILE__, __LINE__, coefficients.rows(), rows));
39  }
40  return false;
41 }
42 
44 
55 Eigen::MatrixXd Polynomial2D::integrateCoeffs(const Eigen::MatrixXd& coefficients, const int& axis = -1, const int& times = 1) {
56  if (times < 0)
57  throw ValueError(format("%s (%d): You have to provide a positive order for integration, %d is not valid. ", __FILE__, __LINE__, times));
58  if (times == 0) return Eigen::MatrixXd(coefficients);
59  Eigen::MatrixXd oldCoefficients;
60  Eigen::MatrixXd newCoefficients(coefficients);
61 
62  switch (axis) {
63  case 0:
64  newCoefficients = Eigen::MatrixXd(coefficients);
65  break;
66  case 1:
67  newCoefficients = Eigen::MatrixXd(coefficients.transpose());
68  break;
69  default:
70  throw ValueError(
71  format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ", __FILE__, __LINE__, axis));
72  break;
73  }
74 
75  std::size_t r, c;
76  for (int k = 0; k < times; k++) {
77  oldCoefficients = Eigen::MatrixXd(newCoefficients);
78  r = oldCoefficients.rows(), c = oldCoefficients.cols();
79  newCoefficients = Eigen::MatrixXd::Zero(r + 1, c);
80  newCoefficients.block(1, 0, r, c) = oldCoefficients.block(0, 0, r, c);
81  for (size_t i = 0; i < r; ++i) {
82  for (size_t j = 0; j < c; ++j)
83  newCoefficients(i + 1, j) /= (i + 1.);
84  }
85  }
86 
87  switch (axis) {
88  case 0:
89  break;
90  case 1:
91  newCoefficients.transposeInPlace();
92  break;
93  default:
94  throw ValueError(
95  format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ", __FILE__, __LINE__, axis));
96  break;
97  }
98 
99  return newCoefficients;
100 }
101 
103 
109 Eigen::MatrixXd Polynomial2D::deriveCoeffs(const Eigen::MatrixXd& coefficients, const int& axis, const int& times) {
110  if (times < 0)
111  throw ValueError(format("%s (%d): You have to provide a positive order for derivation, %d is not valid. ", __FILE__, __LINE__, times));
112  if (times == 0) return Eigen::MatrixXd(coefficients);
113  // Recursion is also possible, but not recommended
114  //Eigen::MatrixXd newCoefficients;
115  //if (times > 1) newCoefficients = deriveCoeffs(coefficients, axis, times-1);
116  //else newCoefficients = Eigen::MatrixXd(coefficients);
117  Eigen::MatrixXd newCoefficients;
118 
119  switch (axis) {
120  case 0:
121  newCoefficients = Eigen::MatrixXd(coefficients);
122  break;
123  case 1:
124  newCoefficients = Eigen::MatrixXd(coefficients.transpose());
125  break;
126  default:
127  throw ValueError(
128  format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ", __FILE__, __LINE__, axis));
129  break;
130  }
131 
132  std::size_t r, c, i, j;
133  for (int k = 0; k < times; k++) {
134  r = newCoefficients.rows(), c = newCoefficients.cols();
135  for (i = 1; i < r; ++i) {
136  for (j = 0; j < c; ++j)
137  newCoefficients(i, j) *= i;
138  }
139  removeRow(newCoefficients, 0);
140  }
141 
142  switch (axis) {
143  case 0:
144  break;
145  case 1:
146  newCoefficients.transposeInPlace();
147  break;
148  default:
149  throw ValueError(
150  format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ", __FILE__, __LINE__, axis));
151  break;
152  }
153 
154  return newCoefficients;
155 }
156 
158 
167 double Polynomial2D::evaluate(const Eigen::MatrixXd& coefficients, const double& x_in) {
168  double result = Eigen::poly_eval(makeVector(coefficients), x_in);
169  if (this->do_debug())
170  std::cout << "Running 1D evaluate(" << mat_to_string(coefficients) << ", x_in:" << vec_to_string(x_in) << "): " << result << std::endl;
171  return result;
172 }
176 double Polynomial2D::evaluate(const Eigen::MatrixXd& coefficients, const double& x_in, const double& y_in) {
177  size_t r = coefficients.rows();
178  double result = evaluate(coefficients.row(r - 1), y_in);
179  for (int i = static_cast<int>(r) - 2; i >= 0; i--) {
180  result *= x_in;
181  result += evaluate(coefficients.row(i), y_in);
182  }
183  if (this->do_debug())
184  std::cout << "Running 2D evaluate(" << mat_to_string(coefficients) << ", x_in:" << vec_to_string(x_in)
185  << ", y_in:" << vec_to_string(y_in) << "): " << result << std::endl;
186  return result;
187 }
188 
193 double Polynomial2D::derivative(const Eigen::MatrixXd& coefficients, const double& x_in, const double& y_in, const int& axis = -1) {
194  return this->evaluate(this->deriveCoeffs(coefficients, axis, 1), x_in, y_in);
195 }
196 
201 double Polynomial2D::integral(const Eigen::MatrixXd& coefficients, const double& x_in, const double& y_in, const int& axis = -1) {
202  return this->evaluate(this->integrateCoeffs(coefficients, axis, 1), x_in, y_in);
203 }
204 
209 double Polynomial2D::solve_limits(Poly2DResidual* res, const double& min, const double& max) {
210  if (do_debug()) std::cout << format("Called solve_limits with: min=%f and max=%f", min, max) << std::endl;
211  double macheps = DBL_EPSILON;
212  double tol = DBL_EPSILON * 1e3;
213  int maxiter = 10;
214  double result = Brent(res, min, max, macheps, tol, maxiter);
215  if (this->do_debug()) std::cout << "Brent solver message: " << res->errstring << std::endl;
216  return result;
217 }
218 
222 double Polynomial2D::solve_guess(Poly2DResidual* res, const double& guess) {
223  if (do_debug()) std::cout << format("Called solve_guess with: guess=%f ", guess) << std::endl;
224  //set_debug_level(1000);
225  double tol = DBL_EPSILON * 1e3;
226  int maxiter = 10;
227  double result = Newton(res, guess, tol, maxiter);
228  if (this->do_debug()) std::cout << "Newton solver message: " << res->errstring << std::endl;
229  return result;
230 }
231 
236 Eigen::VectorXd Polynomial2D::solve(const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const int& axis = -1) {
237  std::size_t r = coefficients.rows(), c = coefficients.cols();
238  Eigen::VectorXd tmpCoefficients;
239  switch (axis) {
240  case 0:
241  tmpCoefficients = Eigen::VectorXd::Zero(r);
242  for (size_t i = 0; i < r; i++) {
243  tmpCoefficients(i, 0) = evaluate(coefficients.row(i), in);
244  }
245  break;
246  case 1:
247  tmpCoefficients = Eigen::VectorXd::Zero(c);
248  for (size_t i = 0; i < c; i++) {
249  tmpCoefficients(i, 0) = evaluate(coefficients.col(i), in);
250  }
251  break;
252  default:
253  throw ValueError(format("%s (%d): You have to provide a dimension, 0 or 1, for the solver, %d is not valid. ", __FILE__, __LINE__, axis));
254  break;
255  }
256  tmpCoefficients(0, 0) -= z_in;
257  if (this->do_debug()) std::cout << "Coefficients: " << mat_to_string(Eigen::MatrixXd(tmpCoefficients)) << std::endl;
258  Eigen::PolynomialSolver<double, Eigen::Dynamic> polySolver(tmpCoefficients);
259  std::vector<double> rootsVec;
260  polySolver.realRoots(rootsVec);
261  if (this->do_debug()) std::cout << "Real roots: " << vec_to_string(rootsVec) << std::endl;
262  return vec_to_eigen(rootsVec);
263  //return rootsVec[0]; // TODO: implement root selection algorithm
264 }
265 
269 double Polynomial2D::solve_limits(const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const double& min, const double& max,
270  const int& axis) {
271  Poly2DResidual res = Poly2DResidual(*this, coefficients, in, z_in, axis);
272  return solve_limits(&res, min, max);
273 }
274 
279 double Polynomial2D::solve_guess(const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const double& guess, const int& axis) {
280  Poly2DResidual res = Poly2DResidual(*this, coefficients, in, z_in, axis);
281  return solve_guess(&res, guess);
282 }
283 
285 
288 double Polynomial2D::simplePolynomial(std::vector<double> const& coefficients, double x) {
289  double result = 0.;
290  for (unsigned int i = 0; i < coefficients.size(); i++) {
291  result += coefficients[i] * pow(x, (int)i);
292  }
293  if (this->do_debug())
294  std::cout << "Running simplePolynomial(" << vec_to_string(coefficients) << ", " << vec_to_string(x) << "): " << result << std::endl;
295  return result;
296 }
297 double Polynomial2D::simplePolynomial(std::vector<std::vector<double>> const& coefficients, double x, double y) {
298  double result = 0;
299  for (unsigned int i = 0; i < coefficients.size(); i++) {
300  result += pow(x, (int)i) * simplePolynomial(coefficients[i], y);
301  }
302  if (this->do_debug())
303  std::cout << "Running simplePolynomial(" << vec_to_string(coefficients) << ", " << vec_to_string(x) << ", " << vec_to_string(y)
304  << "): " << result << std::endl;
305  return result;
306 }
307 
309 
313 double Polynomial2D::baseHorner(std::vector<double> const& coefficients, double x) {
314  double result = 0;
315  for (int i = static_cast<int>(coefficients.size()) - 1; i >= 0; i--) {
316  result *= x;
317  result += coefficients[i];
318  }
319  if (this->do_debug())
320  std::cout << "Running baseHorner(" << vec_to_string(coefficients) << ", " << vec_to_string(x) << "): " << result << std::endl;
321  return result;
322 }
323 
324 double Polynomial2D::baseHorner(std::vector<std::vector<double>> const& coefficients, double x, double y) {
325  double result = 0;
326  for (int i = static_cast<int>(coefficients.size() - 1); i >= 0; i--) {
327  result *= x;
328  result += baseHorner(coefficients[i], y);
329  }
330  if (this->do_debug())
331  std::cout << "Running baseHorner(" << vec_to_string(coefficients) << ", " << vec_to_string(x) << ", " << vec_to_string(y)
332  << "): " << result << std::endl;
333  return result;
334 }
335 
336 Poly2DResidual::Poly2DResidual(Polynomial2D& poly, const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const int& axis) {
337  switch (axis) {
338  case iX:
339  case iY:
340  this->axis = axis;
341  break;
342  default:
343  throw ValueError(format("%s (%d): You have to provide a dimension to the solver, %d is not valid. ", __FILE__, __LINE__, axis));
344  break;
345  }
346 
347  this->poly = poly;
348  this->coefficients = coefficients;
349  this->derIsSet = false;
350  this->in = in;
351  this->z_in = z_in;
352 }
353 
354 double Poly2DResidual::call(double target) {
355  if (axis == iX) return poly.evaluate(coefficients, target, in) - z_in;
356  if (axis == iY) return poly.evaluate(coefficients, in, target) - z_in;
357  return -_HUGE;
358 }
359 
360 double Poly2DResidual::deriv(double target) {
361  if (!this->derIsSet) {
362  this->coefficientsDer = poly.deriveCoeffs(coefficients, axis);
363  this->derIsSet = true;
364  }
365  return poly.evaluate(coefficientsDer, target, in);
366 }
367 
368 // /// Integration functions
369 // /** Integrating coefficients for polynomials is done by dividing the
370 // * original coefficients by (i+1) and elevating the order by 1
371 // * through adding a zero as first coefficient.
372 // * Some reslicing needs to be applied to integrate along the x-axis.
373 // * In the brine/solution equations, reordering of the parameters
374 // * avoids this expensive operation. However, it is included for the
375 // * sake of completeness.
376 // */
377 // /// @param coefficients matrix containing the ordered coefficients
378 // /// @param axis unsigned integer value that represents the desired direction of integration
379 // /// @param times integer value that represents the desired order of integration
380 // /// @param firstExponent integer value that represents the first exponent of the polynomial in axis direction
381 // Eigen::MatrixXd integrateCoeffs(const Eigen::MatrixXd &coefficients, const int &axis, const int &times, const int &firstExponent);
382 //
384 
397 Eigen::MatrixXd Polynomial2DFrac::deriveCoeffs(const Eigen::MatrixXd& coefficients, const int& axis, const int& times, const int& firstExponent) {
398  if (times < 0)
399  throw ValueError(format("%s (%d): You have to provide a positive order for derivation, %d is not valid. ", __FILE__, __LINE__, times));
400  if (times == 0) return Eigen::MatrixXd(coefficients);
401  // Recursion is also possible, but not recommended
402  //Eigen::MatrixXd newCoefficients;
403  //if (times > 1) newCoefficients = deriveCoeffs(coefficients, axis, times-1);
404  //else newCoefficients = Eigen::MatrixXd(coefficients);
405  Eigen::MatrixXd newCoefficients;
406 
407  switch (axis) {
408  case 0:
409  newCoefficients = Eigen::MatrixXd(coefficients);
410  break;
411  case 1:
412  newCoefficients = Eigen::MatrixXd(coefficients.transpose());
413  break;
414  default:
415  throw ValueError(
416  format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ", __FILE__, __LINE__, axis));
417  break;
418  }
419 
420  std::size_t r = newCoefficients.rows(), c = newCoefficients.cols();
421  std::size_t i, j;
422  for (int k = 0; k < times; k++) {
423  for (i = 0; i < r; ++i) {
424  for (j = 0; j < c; ++j) {
425  newCoefficients(i, j) *= double(i) + double(firstExponent);
426  }
427  }
428  }
429 
430  switch (axis) {
431  case 0:
432  break;
433  case 1:
434  newCoefficients.transposeInPlace();
435  break;
436  default:
437  throw ValueError(
438  format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ", __FILE__, __LINE__, axis));
439  break;
440  }
441 
442  return newCoefficients;
443 }
444 
446 
460 double Polynomial2DFrac::evaluate(const Eigen::MatrixXd& coefficients, const double& x_in, const int& firstExponent, const double& x_base) {
461  size_t r = coefficients.rows();
462  size_t c = coefficients.cols();
463 
464  if ((r != 1) && (c != 1)) {
465  throw ValueError(format("%s (%d): You have a 2D coefficient matrix (%d,%d), please use the 2D functions. ", __FILE__, __LINE__,
466  coefficients.rows(), coefficients.cols()));
467  }
468  if ((firstExponent < 0) && (std::abs(x_in - x_base) < CPPOLY_EPSILON)) {
469  //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));
470  const double x_lo = x_base - CPPOLY_DELTA;
471  const double x_hi = x_base + CPPOLY_DELTA;
472  const double y_lo = evaluate(coefficients, x_lo, firstExponent, x_base);
473  const double y_hi = evaluate(coefficients, x_hi, firstExponent, x_base);
474  return (y_hi - y_lo)/(x_hi - x_lo) * (x_in - x_lo) + y_lo;
475  }
476 
477  Eigen::MatrixXd tmpCoeffs(coefficients);
478  if (c == 1) {
479  tmpCoeffs.transposeInPlace();
480  c = r;
481  r = 1;
482  }
483  Eigen::MatrixXd newCoeffs;
484  double negExp = 0; // First we treat the negative exponents
485  double posExp = 0; // then the positive exponents
486 
487  for (int i = 0; i > firstExponent; i--) { // only for firstExponent<0
488  if (c > 0) {
489  negExp += tmpCoeffs(0, 0);
490  removeColumn(tmpCoeffs, 0);
491  }
492  negExp /= x_in - x_base;
493  c = tmpCoeffs.cols();
494  }
495 
496  for (int i = 0; i < firstExponent; i++) { // only for firstExponent>0
497  newCoeffs = Eigen::MatrixXd::Zero(r, c + 1);
498  newCoeffs.block(0, 1, r, c) = tmpCoeffs.block(0, 0, r, c);
499  tmpCoeffs = Eigen::MatrixXd(newCoeffs);
500  c = tmpCoeffs.cols();
501  }
502 
503  if (c > 0) posExp += Eigen::poly_eval(Eigen::RowVectorXd(tmpCoeffs), x_in - x_base);
504  return negExp + posExp;
505 }
506 
514 double Polynomial2DFrac::evaluate(const Eigen::MatrixXd& coefficients, const double& x_in, const double& y_in, const int& x_exp, const int& y_exp,
515  const double& x_base, const double& y_base) {
516  if ((x_exp < 0) && (std::abs(x_in - x_base) < CPPOLY_EPSILON)) {
517  // 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));
518  if (this->do_debug()) std::cout << "Interpolating in x-direction for base " << x_base << " and input " << x_in << std::endl;
519  const double x_lo = x_base - CPPOLY_DELTA;
520  const double x_hi = x_base + CPPOLY_DELTA;
521  const double z_lo = evaluate(coefficients, x_lo, y_in, x_exp, y_exp, x_base, y_base);
522  const double z_hi = evaluate(coefficients, x_hi, y_in, x_exp, y_exp, x_base, y_base);
523  return (z_hi - z_lo)/(x_hi - x_lo) * (x_in - x_lo) + z_lo;
524  }
525  if ((y_exp < 0) && (std::abs(y_in - y_base) < CPPOLY_EPSILON)) {
526  // throw ValueError(format("%s (%d): A fraction cannot be evaluated with zero as denominator, y_in-y_base=%f ", __FILE__, __LINE__, y_in - y_base));
527  if (this->do_debug()) std::cout << "Interpolating in y-direction for base " << y_base << " and input " << y_in << std::endl;
528  const double y_lo = y_base - CPPOLY_DELTA;
529  const double y_hi = y_base + CPPOLY_DELTA;
530  const double z_lo = evaluate(coefficients, x_in, y_lo, x_exp, y_exp, x_base, y_base);
531  const double z_hi = evaluate(coefficients, x_in, y_hi, x_exp, y_exp, x_base, y_base);
532  return (z_hi - z_lo)/(y_hi - y_lo) * (y_in - y_lo) + z_lo;
533  }
534 
535  Eigen::MatrixXd tmpCoeffs(coefficients);
536  Eigen::MatrixXd newCoeffs;
537  size_t r = tmpCoeffs.rows();
538  size_t c = tmpCoeffs.cols();
539  double negExp = 0; // First we treat the negative exponents
540  double posExp = 0; // then the positive exponents
541 
542  for (int i = 0; i > x_exp; i--) { // only for x_exp<0
543  r = tmpCoeffs.rows();
544  if (r > 0) {
545  negExp += evaluate(tmpCoeffs.row(0), y_in, y_exp, y_base);
546  removeRow(tmpCoeffs, 0);
547  }
548  negExp /= x_in - x_base;
549  }
550 
551  r = tmpCoeffs.rows();
552  for (int i = 0; i < x_exp; i++) { // only for x_exp>0
553  newCoeffs = Eigen::MatrixXd::Zero(r + 1, c);
554  newCoeffs.block(1, 0, r, c) = tmpCoeffs.block(0, 0, r, c);
555  tmpCoeffs = Eigen::MatrixXd(newCoeffs);
556  r += 1; // r = tmpCoeffs.rows();
557  }
558 
559  //r = tmpCoeffs.rows();
560  if (r > 0) posExp += evaluate(tmpCoeffs.row(r - 1), y_in, y_exp, y_base);
561  for (int i = static_cast<int>(r) - 2; i >= 0; i--) {
562  posExp *= x_in - x_base;
563  posExp += evaluate(tmpCoeffs.row(i), y_in, y_exp, y_base);
564  }
565  if (this->do_debug()) std::cout << "Running 2D evaluate(" << mat_to_string(coefficients) << ", " << std::endl;
566  if (this->do_debug())
567  std::cout << "x_in:" << vec_to_string(x_in) << ", y_in:" << vec_to_string(y_in) << ", x_exp:" << vec_to_string(x_exp)
568  << ", y_exp:" << vec_to_string(y_exp) << ", x_base:" << vec_to_string(x_base) << ", y_base:" << vec_to_string(y_base)
569  << "): " << negExp + posExp << std::endl;
570  return negExp + posExp;
571 }
572 
581 double Polynomial2DFrac::derivative(const Eigen::MatrixXd& coefficients, const double& x_in, const double& y_in, const int& axis, const int& x_exp,
582  const int& y_exp, const double& x_base, const double& y_base) {
583  Eigen::MatrixXd newCoefficients;
584  int der_exp, other_exp;
585  double der_val, other_val;
586  double int_base, other_base;
587 
588  switch (axis) {
589  case 0:
590  newCoefficients = Eigen::MatrixXd(coefficients);
591  der_exp = x_exp;
592  other_exp = y_exp;
593  der_val = x_in;
594  other_val = y_in;
595  int_base = x_base;
596  other_base = y_base;
597  break;
598  case 1:
599  newCoefficients = Eigen::MatrixXd(coefficients.transpose());
600  der_exp = y_exp;
601  other_exp = x_exp;
602  der_val = y_in;
603  other_val = x_in;
604  int_base = y_base;
605  other_base = x_base;
606  break;
607  default:
608  throw ValueError(
609  format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ", __FILE__, __LINE__, axis));
610  break;
611  }
612 
613  const int times = 1;
614  newCoefficients = deriveCoeffs(newCoefficients, 0, times, der_exp);
615  der_exp -= times;
616 
617  return evaluate(newCoefficients, der_val, other_val, der_exp, other_exp, int_base, other_base);
618 }
619 
629 double Polynomial2DFrac::integral(const Eigen::MatrixXd& coefficients, const double& x_in, const double& y_in, const int& axis, const int& x_exp,
630  const int& y_exp, const double& x_base, const double& y_base, const double& ax_val) {
631 
632  Eigen::MatrixXd newCoefficients;
633  int int_exp, other_exp;
634  double int_val, other_val;
635  double int_base, other_base;
636 
637  switch (axis) {
638  case 0:
639  newCoefficients = Eigen::MatrixXd(coefficients);
640  int_exp = x_exp;
641  other_exp = y_exp;
642  int_val = x_in;
643  other_val = y_in;
644  int_base = x_base;
645  other_base = y_base;
646  break;
647  case 1:
648  newCoefficients = Eigen::MatrixXd(coefficients.transpose());
649  int_exp = y_exp;
650  other_exp = x_exp;
651  int_val = y_in;
652  other_val = x_in;
653  int_base = y_base;
654  other_base = x_base;
655  break;
656  default:
657  throw ValueError(
658  format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ", __FILE__, __LINE__, axis));
659  break;
660  }
661 
662  if (int_exp < -1)
663  throw NotImplementedError(
664  format("%s (%d): This function is only implemented for lowest exponents >= -1, int_exp=%d ", __FILE__, __LINE__, int_exp));
665  // TODO: Fix this and allow the direct calculation of integrals
666  if (std::abs(ax_val) > DBL_EPSILON)
667  throw NotImplementedError(
668  format("%s (%d): This function is only implemented for indefinite integrals, ax_val=%d ", __FILE__, __LINE__, ax_val));
669 
670  size_t r = newCoefficients.rows();
671  size_t c = newCoefficients.cols();
672 
673  if (int_exp == -1) {
674  if (std::abs(int_base) < DBL_EPSILON) {
675  Eigen::MatrixXd tmpCoefficients = newCoefficients.row(0) * log(int_val - int_base);
676  newCoefficients = integrateCoeffs(newCoefficients.block(1, 0, r - 1, c), 0, 1);
677  newCoefficients.row(0) = tmpCoefficients;
678  return evaluate(newCoefficients, int_val, other_val, 0, other_exp, int_base, other_base);
679  } else {
680  // Reduce the coefficients to the integration dimension:
681  newCoefficients = Eigen::MatrixXd(r, 1);
682  for (std::size_t i = 0; i < r; i++) {
683  newCoefficients(i, 0) = evaluate(coefficients.row(i), other_val, other_exp, other_base);
684  }
685  return fracIntCentral(newCoefficients.transpose(), int_val, int_base);
686  }
687  }
688 
689  Eigen::MatrixXd tmpCoeffs;
690  r = newCoefficients.rows();
691  for (int i = 0; i < int_exp; i++) { // only for x_exp>0
692  tmpCoeffs = Eigen::MatrixXd::Zero(r + 1, c);
693  tmpCoeffs.block(1, 0, r, c) = newCoefficients.block(0, 0, r, c);
694  newCoefficients = Eigen::MatrixXd(tmpCoeffs);
695  r += 1; // r = newCoefficients.rows();
696  }
697 
698  return evaluate(integrateCoeffs(newCoefficients, 0, 1), int_val, other_val, 0, other_exp, int_base, other_base);
699 }
700 
710 Eigen::VectorXd Polynomial2DFrac::solve(const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const int& axis, const int& x_exp,
711  const int& y_exp, const double& x_base, const double& y_base) {
712 
713  Eigen::MatrixXd newCoefficients;
714  Eigen::VectorXd tmpCoefficients;
715  int solve_exp, other_exp;
716  double input;
717 
718  switch (axis) {
719  case 0:
720  newCoefficients = Eigen::MatrixXd(coefficients);
721  solve_exp = x_exp;
722  other_exp = y_exp;
723  input = in - y_base;
724  break;
725  case 1:
726  newCoefficients = Eigen::MatrixXd(coefficients.transpose());
727  solve_exp = y_exp;
728  other_exp = x_exp;
729  input = in - x_base;
730  break;
731  default:
732  throw ValueError(format("%s (%d): You have to provide a dimension, 0 or 1, for the solver, %d is not valid. ", __FILE__, __LINE__, axis));
733  break;
734  }
735 
736  if (this->do_debug()) std::cout << "Coefficients: " << mat_to_string(Eigen::MatrixXd(newCoefficients)) << std::endl;
737 
738  const size_t r = newCoefficients.rows();
739  for (size_t i = 0; i < r; i++) {
740  newCoefficients(i, 0) = evaluate(newCoefficients.row(i), input, other_exp);
741  }
742 
743  //Eigen::VectorXd tmpCoefficients;
744  if (solve_exp >= 0) { // extend vector to zero exponent
745  tmpCoefficients = Eigen::VectorXd::Zero(r + solve_exp);
746  tmpCoefficients.block(solve_exp, 0, r, 1) = newCoefficients.block(0, 0, r, 1);
747  tmpCoefficients(0, 0) -= z_in;
748  } else { // check if vector reaches to zero exponent
749  int diff = 1 - static_cast<int>(r) - solve_exp; // How many entries have to be added
750  tmpCoefficients = Eigen::VectorXd::Zero(r + std::max(diff, 0));
751  tmpCoefficients.block(0, 0, r, 1) = newCoefficients.block(0, 0, r, 1);
752  tmpCoefficients(r + diff - 1, 0) -= z_in;
753  throw NotImplementedError(format("%s (%d): Currently, there is no solver for the generalised polynomial, an exponent of %d is not valid. ",
754  __FILE__, __LINE__, solve_exp));
755  }
756 
757  if (this->do_debug()) std::cout << "Coefficients: " << mat_to_string(Eigen::MatrixXd(tmpCoefficients)) << std::endl;
758  Eigen::PolynomialSolver<double, -1> polySolver(tmpCoefficients);
759  std::vector<double> rootsVec;
760  polySolver.realRoots(rootsVec);
761  if (this->do_debug()) std::cout << "Real roots: " << vec_to_string(rootsVec) << std::endl;
762  return vec_to_eigen(rootsVec);
763  //return rootsVec[0]; // TODO: implement root selection algorithm
764 }
765 
777 double Polynomial2DFrac::solve_limits(const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const double& min, const double& max,
778  const int& axis, const int& x_exp, const int& y_exp, const double& x_base, const double& y_base) {
779  if (do_debug())
780  std::cout << format("Called solve_limits with: %f, %f, %f, %f, %d, %d, %d, %f, %f", in, z_in, min, max, axis, x_exp, y_exp, x_base, y_base)
781  << std::endl;
782  Poly2DFracResidual res = Poly2DFracResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base);
783  return Polynomial2D::solve_limits(&res, min, max);
784 } //TODO: Implement tests for this solver
785 
796 double Polynomial2DFrac::solve_guess(const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const double& guess, const int& axis,
797  const int& x_exp, const int& y_exp, const double& x_base, const double& y_base) {
798  if (do_debug())
799  std::cout << format("Called solve_guess with: %f, %f, %f, %d, %d, %d, %f, %f", in, z_in, guess, axis, x_exp, y_exp, x_base, y_base)
800  << std::endl;
801  Poly2DFracResidual res = Poly2DFracResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base);
802  return Polynomial2D::solve_guess(&res, guess);
803 } //TODO: Implement tests for this solver
804 
817 double Polynomial2DFrac::solve_limitsInt(const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const double& min,
818  const double& max, const int& axis, const int& x_exp, const int& y_exp, const double& x_base,
819  const double& y_base, const int& int_axis) {
820  Poly2DFracIntResidual res = Poly2DFracIntResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base, int_axis);
821  return Polynomial2D::solve_limits(&res, min, max);
822 } //TODO: Implement tests for this solver
823 
835 double Polynomial2DFrac::solve_guessInt(const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const double& guess,
836  const int& axis, const int& x_exp, const int& y_exp, const double& x_base, const double& y_base,
837  const int& int_axis) {
838  Poly2DFracIntResidual res = Poly2DFracIntResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base, int_axis);
839  return Polynomial2D::solve_guess(&res, guess);
840 } //TODO: Implement tests for this solver
841 
850 //Helper functions to calculate binomial coefficients:
851 //http://rosettacode.org/wiki/Evaluate_binomial_coefficients#C.2B.2B
853 double Polynomial2DFrac::factorial(const int& nValue) {
854  double value = 1;
855  for (int i = 2; i <= nValue; i++)
856  value = value * i;
857  return value;
858 }
861 double Polynomial2DFrac::binom(const int& nValue, const int& nValue2) {
862  if (nValue2 == 1) return nValue * 1.0;
863  return (factorial(nValue)) / (factorial(nValue2) * factorial((nValue - nValue2)));
864 }
869 Eigen::MatrixXd Polynomial2DFrac::fracIntCentralDvector(const int& m, const double& x_in, const double& x_base) {
870  if (m < 1) throw ValueError(format("%s (%d): You have to provide coefficients, a vector length of %d is not a valid. ", __FILE__, __LINE__, m));
871 
872  Eigen::MatrixXd D = Eigen::MatrixXd::Zero(1, m);
873  double tmp;
874  // TODO: This can be optimized using the Horner scheme!
875  for (int j = 0; j < m; j++) { // loop through row
876  tmp = pow(-1.0, j) * log(x_in) * pow(x_base, j);
877  for (int k = 0; k < j; k++) { // internal loop for every entry
878  tmp += binom(j, k) * pow(-1.0, k) / (j - k) * pow(x_in, j - k) * pow(x_base, k);
879  }
880  D(0, j) = tmp;
881  }
882  return D;
883 }
888 double Polynomial2DFrac::fracIntCentral(const Eigen::MatrixXd& coefficients, const double& x_in, const double& x_base) {
889  if (coefficients.rows() != 1) {
890  throw ValueError(format("%s (%d): You have a 2D coefficient matrix (%d,%d), please use the 2D functions. ", __FILE__, __LINE__,
891  coefficients.rows(), coefficients.cols()));
892  }
893  int m = static_cast<int>(coefficients.cols());
894  Eigen::MatrixXd D = fracIntCentralDvector(m, x_in, x_base);
895  double result = 0;
896  for (int j = 0; j < m; j++) {
897  result += coefficients(0, j) * D(0, j);
898  }
899  if (this->do_debug())
900  std::cout << "Running fracIntCentral(" << mat_to_string(coefficients) << ", " << vec_to_string(x_in) << ", " << vec_to_string(x_base)
901  << "): " << result << std::endl;
902  return result;
903 }
904 
905 Poly2DFracResidual::Poly2DFracResidual(Polynomial2DFrac& poly, const Eigen::MatrixXd& coefficients, const double& in, const double& z_in,
906  const int& axis, const int& x_exp, const int& y_exp, const double& x_base, const double& y_base)
907  : Poly2DResidual(poly, coefficients, in, z_in, axis) {
908  this->x_exp = x_exp;
909  this->y_exp = y_exp;
910  this->x_base = x_base;
911  this->y_base = y_base;
912 }
913 
914 double Poly2DFracResidual::call(double target) {
915  if (axis == iX) return poly.evaluate(coefficients, target, in, x_exp, y_exp, x_base, y_base) - z_in;
916  if (axis == iY) return poly.evaluate(coefficients, in, target, x_exp, y_exp, x_base, y_base) - z_in;
917  return _HUGE;
918 }
919 
920 double Poly2DFracResidual::deriv(double target) {
921  if (axis == iX) return poly.derivative(coefficients, target, in, axis, x_exp, y_exp, x_base, y_base);
922  if (axis == iY) return poly.derivative(coefficients, in, target, axis, x_exp, y_exp, x_base, y_base);
923  return _HUGE;
924 }
925 
926 Poly2DFracIntResidual::Poly2DFracIntResidual(Polynomial2DFrac& poly, const Eigen::MatrixXd& coefficients, const double& in, const double& z_in,
927  const int& axis, const int& x_exp, const int& y_exp, const double& x_base, const double& y_base,
928  const int& int_axis)
929  : Poly2DFracResidual(poly, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base) {
930  this->int_axis = int_axis;
931 }
932 
933 double Poly2DFracIntResidual::call(double target) {
934  if (axis == iX) return poly.integral(coefficients, target, in, int_axis, x_exp, y_exp, x_base, y_base) - z_in;
935  if (axis == iY) return poly.integral(coefficients, in, target, int_axis, x_exp, y_exp, x_base, y_base) - z_in;
936  return _HUGE;
937 }
938 
939 double Poly2DFracIntResidual::deriv(double target) {
940  if (axis == iX) return poly.evaluate(coefficients, target, in, x_exp, y_exp, x_base, y_base);
941  if (axis == iY) return poly.evaluate(coefficients, in, target, x_exp, y_exp, x_base, y_base);
942  return _HUGE;
943 }
944 
945 } // namespace CoolProp
946 
947 #ifdef ENABLE_CATCH
948 # include <math.h>
949 # include <iostream>
950 # include <catch2/catch_all.hpp>
951 
952 TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp", "[PolyMath]") {
953  bool PRINT = false;
954  std::string tmpStr;
955 
957  std::vector<double> cHeat;
958  cHeat.clear();
959  cHeat.push_back(+1.1562261074E+03);
960  cHeat.push_back(+2.0994549103E+00);
961  cHeat.push_back(+7.7175381057E-07);
962  cHeat.push_back(-3.7008444051E-20);
963 
964  double deltaT = 0.1;
965  double Tmin = 273.15 - 50;
966  double Tmax = 273.15 + 250;
967  double Tinc = 200;
968 
969  std::vector<std::vector<double>> cHeat2D;
970  cHeat2D.push_back(cHeat);
971  cHeat2D.push_back(cHeat);
972  cHeat2D.push_back(cHeat);
973 
974  Eigen::MatrixXd matrix2D = CoolProp::vec_to_eigen(cHeat2D);
975 
976  Eigen::MatrixXd matrix2Dtmp;
977  std::vector<std::vector<double>> vec2Dtmp;
978 
979  SECTION("Coefficient parsing") {
980  CoolProp::Polynomial2D poly;
981  CHECK_THROWS(poly.checkCoefficients(matrix2D, 4, 5));
982  CHECK(poly.checkCoefficients(matrix2D, 3, 4));
983  }
984 
985  SECTION("Coefficient operations") {
986  Eigen::MatrixXd matrix;
987  CoolProp::Polynomial2D poly;
988 
989  CHECK_THROWS(poly.integrateCoeffs(matrix2D));
990 
991  CHECK_NOTHROW(matrix = poly.integrateCoeffs(matrix2D, 0));
992  tmpStr = CoolProp::mat_to_string(matrix2D);
993  if (PRINT) std::cout << tmpStr << std::endl;
994  tmpStr = CoolProp::mat_to_string(matrix);
995  if (PRINT) std::cout << tmpStr << std::endl << std::endl;
996 
997  CHECK_NOTHROW(matrix = poly.integrateCoeffs(matrix2D, 1));
998  tmpStr = CoolProp::mat_to_string(matrix2D);
999  if (PRINT) std::cout << tmpStr << std::endl;
1000  tmpStr = CoolProp::mat_to_string(matrix);
1001  if (PRINT) std::cout << tmpStr << std::endl << std::endl;
1002 
1003  CHECK_THROWS(poly.deriveCoeffs(matrix2D));
1004 
1005  CHECK_NOTHROW(matrix = poly.deriveCoeffs(matrix2D, 0));
1006  tmpStr = CoolProp::mat_to_string(matrix2D);
1007  if (PRINT) std::cout << tmpStr << std::endl;
1008  tmpStr = CoolProp::mat_to_string(matrix);
1009  if (PRINT) std::cout << tmpStr << std::endl << std::endl;
1010 
1011  CHECK_NOTHROW(matrix = poly.deriveCoeffs(matrix2D, 1));
1012  tmpStr = CoolProp::mat_to_string(matrix2D);
1013  if (PRINT) std::cout << tmpStr << std::endl;
1014  tmpStr = CoolProp::mat_to_string(matrix);
1015  if (PRINT) std::cout << tmpStr << std::endl << std::endl;
1016  }
1017 
1018  SECTION("Evaluation and test values") {
1019 
1020  Eigen::MatrixXd matrix = CoolProp::vec_to_eigen(cHeat);
1021  CoolProp::Polynomial2D poly;
1022 
1023  double acc = 0.0001;
1024 
1025  double T = 273.15 + 50;
1026  double c = poly.evaluate(matrix, T, 0.0);
1027  double d = 1834.746;
1028 
1029  {
1030  CAPTURE(T);
1031  CAPTURE(c);
1032  CAPTURE(d);
1033  tmpStr = CoolProp::mat_to_string(matrix);
1034  CAPTURE(tmpStr);
1035  CHECK(check_abs(c, d, acc));
1036  }
1037 
1038  c = 2.0;
1039  c = poly.solve(matrix, 0.0, d, 0)[0];
1040  {
1041  CAPTURE(T);
1042  CAPTURE(c);
1043  CAPTURE(d);
1044  CHECK(check_abs(c, T, acc));
1045  }
1046 
1047  c = 2.0;
1048  c = poly.solve_limits(matrix, 0.0, d, -50, 750, 0);
1049  {
1050  CAPTURE(T);
1051  CAPTURE(c);
1052  CAPTURE(d);
1053  CHECK(check_abs(c, T, acc));
1054  }
1055 
1056  c = 2.0;
1057  c = poly.solve_guess(matrix, 0.0, d, 350, 0);
1058  {
1059  CAPTURE(T);
1060  CAPTURE(c);
1061  CAPTURE(d);
1062  CHECK(check_abs(c, T, acc));
1063  }
1064 
1065  // T = 0.0;
1066  // solve.setGuess(75+273.15);
1067  // T = solve.polyval(cHeat,c);
1068  // printf("Should be : T = %3.3f \t K \n",273.15+50.0);
1069  // printf("From object: T = %3.3f \t K \n",T);
1070  //
1071  // T = 0.0;
1072  // solve.setLimits(273.15+10,273.15+100);
1073  // T = solve.polyval(cHeat,c);
1074  // printf("Should be : T = %3.3f \t K \n",273.15+50.0);
1075  // printf("From object: T = %3.3f \t K \n",T);
1076  }
1077 
1078  SECTION("Integration and derivation tests") {
1079 
1080  CoolProp::Polynomial2D poly;
1081 
1082  Eigen::MatrixXd matrix(matrix2D);
1083  Eigen::MatrixXd matrixInt = poly.integrateCoeffs(matrix, 1);
1084  Eigen::MatrixXd matrixDer = poly.deriveCoeffs(matrix, 1);
1085  Eigen::MatrixXd matrixInt2 = poly.integrateCoeffs(matrix, 1, 2);
1086  Eigen::MatrixXd matrixDer2 = poly.deriveCoeffs(matrix, 1, 2);
1087 
1088  CHECK_THROWS(poly.evaluate(matrix, 0.0));
1089 
1090  double x = 0.3, y = 255.3, val1, val2, val3, val4;
1091 
1092  //CHECK( std::abs( polyInt.derivative(x,y,0)-poly2D.evaluate(x,y) ) <= 1e-10 );
1093 
1094  std::string tmpStr;
1095 
1096  double acc = 0.001;
1097 
1098  for (double T = Tmin; T < Tmax; T += Tinc) {
1099  val1 = poly.evaluate(matrix, x, T - deltaT);
1100  val2 = poly.evaluate(matrix, x, T + deltaT);
1101  val3 = (val2 - val1) / 2 / deltaT;
1102  val4 = poly.evaluate(matrixDer, x, T);
1103  CAPTURE(T);
1104  CAPTURE(val3);
1105  CAPTURE(val4);
1106  tmpStr = CoolProp::mat_to_string(matrix);
1107  CAPTURE(tmpStr);
1108  tmpStr = CoolProp::mat_to_string(matrixDer);
1109  CAPTURE(tmpStr);
1110  CHECK(check_abs(val3, val4, acc));
1111  }
1112 
1113  for (double T = Tmin; T < Tmax; T += Tinc) {
1114  val1 = poly.evaluate(matrixDer, x, T - deltaT);
1115  val2 = poly.evaluate(matrixDer, x, T + deltaT);
1116  val3 = (val2 - val1) / 2 / deltaT;
1117  val4 = poly.evaluate(matrixDer2, x, T);
1118  CAPTURE(T);
1119  CAPTURE(val3);
1120  CAPTURE(val4);
1121  tmpStr = CoolProp::mat_to_string(matrixDer);
1122  CAPTURE(tmpStr);
1123  tmpStr = CoolProp::mat_to_string(matrixDer2);
1124  CAPTURE(tmpStr);
1125  CHECK(check_abs(val3, val4, acc));
1126  }
1127 
1128  for (double T = Tmin; T < Tmax; T += Tinc) {
1129  val1 = poly.evaluate(matrixInt, x, T - deltaT);
1130  val2 = poly.evaluate(matrixInt, x, T + deltaT);
1131  val3 = (val2 - val1) / 2 / deltaT;
1132  val4 = poly.evaluate(matrix, x, T);
1133  CAPTURE(T);
1134  CAPTURE(val3);
1135  CAPTURE(val4);
1136  tmpStr = CoolProp::mat_to_string(matrixInt);
1137  CAPTURE(tmpStr);
1138  tmpStr = CoolProp::mat_to_string(matrix);
1139  CAPTURE(tmpStr);
1140  CHECK(check_abs(val3, val4, acc));
1141  }
1142 
1143  for (double T = Tmin; T < Tmax; T += Tinc) {
1144  val1 = poly.evaluate(matrixInt2, x, T - deltaT);
1145  val2 = poly.evaluate(matrixInt2, x, T + deltaT);
1146  val3 = (val2 - val1) / 2 / deltaT;
1147  val4 = poly.evaluate(matrixInt, x, T);
1148  CAPTURE(T);
1149  CAPTURE(val3);
1150  CAPTURE(val4);
1151  tmpStr = CoolProp::mat_to_string(matrixInt2);
1152  CAPTURE(tmpStr);
1153  tmpStr = CoolProp::mat_to_string(matrixInt);
1154  CAPTURE(tmpStr);
1155  CHECK(check_abs(val3, val4, acc));
1156  }
1157 
1158  for (double T = Tmin; T < Tmax; T += Tinc) {
1159  val1 = poly.evaluate(matrix, x, T);
1160  val2 = poly.derivative(matrixInt, x, T, 1);
1161  CAPTURE(T);
1162  CAPTURE(val1);
1163  CAPTURE(val2);
1164  CHECK(check_abs(val1, val2, acc));
1165  }
1166 
1167  for (double T = Tmin; T < Tmax; T += Tinc) {
1168  val1 = poly.derivative(matrix, x, T, 1);
1169  val2 = poly.evaluate(matrixDer, x, T);
1170  CAPTURE(T);
1171  CAPTURE(val1);
1172  CAPTURE(val2);
1173  CHECK(check_abs(val1, val2, acc));
1174  }
1175  }
1176 
1177  SECTION("Testing Polynomial2DFrac") {
1178 
1179  Eigen::MatrixXd matrix = CoolProp::vec_to_eigen(cHeat);
1180  CoolProp::Polynomial2D poly;
1181  CoolProp::Polynomial2DFrac frac;
1182 
1183  double acc = 0.0001;
1184 
1185  double T = 273.15 + 50;
1186  double a, b;
1187  double c = poly.evaluate(matrix, T, 0.0);
1188  double d = frac.evaluate(matrix, T, 0.0, 0, 0);
1189 
1190  {
1191  CAPTURE(T);
1192  CAPTURE(c);
1193  CAPTURE(d);
1194  tmpStr = CoolProp::mat_to_string(matrix);
1195  CAPTURE(tmpStr);
1196  CHECK(check_abs(c, d, acc));
1197  }
1198 
1199  c = poly.evaluate(matrix, T, 0.0) / T / T;
1200  d = frac.evaluate(matrix, T, 0.0, -2, 0);
1201 
1202  {
1203  CAPTURE(T);
1204  CAPTURE(c);
1205  CAPTURE(d);
1206  tmpStr = CoolProp::mat_to_string(matrix);
1207  CAPTURE(tmpStr);
1208  CHECK(check_abs(c, d, acc));
1209  }
1210 
1211  matrix = CoolProp::vec_to_eigen(cHeat2D);
1212  double y = 0.1;
1213  c = poly.evaluate(matrix, T, y) / T / T;
1214  d = frac.evaluate(matrix, T, y, -2, 0);
1215 
1216  {
1217  CAPTURE(T);
1218  CAPTURE(c);
1219  CAPTURE(d);
1220  tmpStr = CoolProp::mat_to_string(matrix);
1221  CAPTURE(tmpStr);
1222  CHECK(check_abs(c, d, acc));
1223  }
1224 
1225  c = poly.evaluate(matrix, T, y) / y / y;
1226  d = frac.evaluate(matrix, T, y, 0, -2);
1227 
1228  {
1229  CAPTURE(T);
1230  CAPTURE(c);
1231  CAPTURE(d);
1232  tmpStr = CoolProp::mat_to_string(matrix);
1233  CAPTURE(tmpStr);
1234  CHECK(check_abs(c, d, acc));
1235  }
1236 
1237  c = poly.evaluate(matrix, T, y) / T / T / y / y;
1238  d = frac.evaluate(matrix, T, y, -2, -2);
1239 
1240  {
1241  CAPTURE(T);
1242  CAPTURE(c);
1243  CAPTURE(d);
1244  tmpStr = CoolProp::mat_to_string(matrix);
1245  CAPTURE(tmpStr);
1246  CHECK(check_abs(c, d, acc));
1247  }
1248 
1249  c = poly.evaluate(matrix, T, y) / T / T * y * y;
1250  d = frac.evaluate(matrix, T, y, -2, 2);
1251 
1252  {
1253  CAPTURE(T);
1254  CAPTURE(c);
1255  CAPTURE(d);
1256  tmpStr = CoolProp::mat_to_string(matrix);
1257  CAPTURE(tmpStr);
1258  CHECK(check_abs(c, d, acc));
1259  }
1260 
1261  matrix = CoolProp::vec_to_eigen(cHeat);
1262  T = 273.15 + 50;
1263  c = 145.59157247249246;
1264  d = frac.integral(matrix, T, 0.0, 0, -1, 0) - frac.integral(matrix, 273.15 + 25, 0.0, 0, -1, 0);
1265 
1266  {
1267  CAPTURE(T);
1268  CAPTURE(c);
1269  CAPTURE(d);
1270  tmpStr = CoolProp::mat_to_string(matrix);
1271  CAPTURE(tmpStr);
1272  CHECK(check_abs(c, d, acc));
1273  }
1274 
1275  T = 423.15;
1276  c = 3460.895272;
1277  d = frac.integral(matrix, T, 0.0, 0, -1, 0, 348.15, 0.0);
1278 
1279  {
1280  CAPTURE(T);
1281  CAPTURE(c);
1282  CAPTURE(d);
1283  tmpStr = CoolProp::mat_to_string(matrix);
1284  CAPTURE(tmpStr);
1285  CHECK(check_abs(c, d, acc));
1286  }
1287 
1288  deltaT = 0.01;
1289  for (T = Tmin; T < Tmax; T += Tinc) {
1290  a = poly.evaluate(matrix, T - deltaT, y);
1291  b = poly.evaluate(matrix, T + deltaT, y);
1292  c = (b - a) / 2.0 / deltaT;
1293  d = frac.derivative(matrix, T, y, 0, 0, 0);
1294  CAPTURE(a);
1295  CAPTURE(b);
1296  CAPTURE(T);
1297  CAPTURE(c);
1298  CAPTURE(d);
1299  tmpStr = CoolProp::mat_to_string(matrix);
1300  CAPTURE(tmpStr);
1301  CHECK(check_abs(c, d, acc));
1302  }
1303 
1304  T = 273.15 + 150;
1305  c = -2.100108045;
1306  d = frac.derivative(matrix, T, 0.0, 0, 0, 0);
1307  {
1308  CAPTURE(T);
1309  CAPTURE(c);
1310  CAPTURE(d);
1311  tmpStr = CoolProp::mat_to_string(matrix);
1312  CAPTURE(tmpStr);
1313  CHECK(check_abs(c, d, acc));
1314  }
1315 
1316  c = -0.006456574589;
1317  d = frac.derivative(matrix, T, 0.0, 0, -1, 0);
1318  {
1319  CAPTURE(T);
1320  CAPTURE(c);
1321  CAPTURE(d);
1322  tmpStr = CoolProp::mat_to_string(matrix);
1323  CAPTURE(tmpStr);
1324  CHECK(check_abs(c, d, acc));
1325  }
1326 
1327  c = frac.evaluate(matrix, T, 0.0, 2, 0);
1328  d = frac.solve(matrix, 0.0, c, 0, 2, 0)[0];
1329  {
1330  CAPTURE(T);
1331  CAPTURE(c);
1332  CAPTURE(d);
1333  tmpStr = CoolProp::mat_to_string(matrix);
1334  CAPTURE(tmpStr);
1335  CHECK(check_abs(T, d, acc));
1336  }
1337 
1338  c = frac.evaluate(matrix, T, 0.0, 0, 0);
1339  d = frac.solve(matrix, 0.0, c, 0, 0, 0)[0];
1340  {
1341  CAPTURE(T);
1342  CAPTURE(c);
1343  CAPTURE(d);
1344  tmpStr = CoolProp::mat_to_string(matrix);
1345  CAPTURE(tmpStr);
1346  CHECK(check_abs(T, d, acc));
1347  }
1348 
1349  c = frac.evaluate(matrix, T, 0.0, -1, 0);
1350  CHECK_THROWS(d = frac.solve(matrix, 0.0, c, 0, -1, 0)[0]);
1351  // {
1352  // CAPTURE(T);
1353  // CAPTURE(c);
1354  // CAPTURE(d);
1355  // tmpStr = CoolProp::mat_to_string(matrix);
1356  // CAPTURE(tmpStr);
1357  // CHECK( check_abs(T,d,acc) );
1358  // }
1359 
1360  d = frac.solve_limits(matrix, 0.0, c, T - 10, T + 10, 0, -1, 0);
1361  {
1362  CAPTURE(T);
1363  CAPTURE(c);
1364  CAPTURE(d);
1365  tmpStr = CoolProp::mat_to_string(matrix);
1366  CAPTURE(tmpStr);
1367  CHECK(check_abs(T, d, acc));
1368  }
1369 
1370  d = frac.solve_guess(matrix, 0.0, c, T - 10, 0, -1, 0);
1371  {
1372  CAPTURE(T);
1373  CAPTURE(c);
1374  CAPTURE(d);
1375  tmpStr = CoolProp::mat_to_string(matrix);
1376  CAPTURE(tmpStr);
1377  CHECK(check_abs(T, d, acc));
1378  }
1379 
1380  c = -0.00004224550082;
1381  d = frac.derivative(matrix, T, 0.0, 0, -2, 0);
1382  {
1383  CAPTURE(T);
1384  CAPTURE(c);
1385  CAPTURE(d);
1386  tmpStr = CoolProp::mat_to_string(matrix);
1387  CAPTURE(tmpStr);
1388  CHECK(check_abs(c, d, acc));
1389  }
1390 
1391  c = frac.evaluate(matrix, T, 0.0, 0, 0, 0.0, 0.0);
1392  d = frac.solve(matrix, 0.0, c, 0, 0, 0, 0.0, 0.0)[0];
1393  {
1394  CAPTURE(T);
1395  CAPTURE(c);
1396  CAPTURE(d);
1397  tmpStr = CoolProp::mat_to_string(matrix);
1398  CAPTURE(tmpStr);
1399  tmpStr = CoolProp::mat_to_string(Eigen::MatrixXd(frac.solve(matrix, 0.0, c, 0, 0, 0, 250, 0.0)));
1400  CAPTURE(tmpStr);
1401  CHECK(check_abs(T, d, acc));
1402  }
1403  }
1404 }
1405 
1406 #endif /* ENABLE_CATCH */
Generated on Wed Nov 29 2023 14:26:24 for CoolProp by doxygen 1.9.1