CoolProp 6.8.1dev
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
18namespace CoolProp {
19
20constexpr double CPPOLY_EPSILON = DBL_EPSILON * 100.0;
21constexpr double CPPOLY_DELTA = CPPOLY_EPSILON * 10.0;
22
24
30bool 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
55Eigen::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
109Eigen::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
167double 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}
176double 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
193double 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
201double 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
209double 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
222double 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
236Eigen::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
269double 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
279double 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
288double 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}
297double 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
313double 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
324double 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
336Poly2DResidual::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
354double 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
360double 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
397Eigen::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
460double 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
514double 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;
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
581double 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
629double 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
710Eigen::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
777double 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
796double 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
817double 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
835double 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
853double 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}
861double 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}
869Eigen::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}
888double 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
905Poly2DFracResidual::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
914double 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
920double 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
926Poly2DFracIntResidual::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
933double 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
939double 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
952TEST_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 = poly.solve(matrix, 0.0, d, 0)[0];
1039 {
1040 CAPTURE(T);
1041 CAPTURE(c);
1042 CAPTURE(d);
1043 CHECK(check_abs(c, T, acc));
1044 }
1045
1046 c = poly.solve_limits(matrix, 0.0, d, -50, 750, 0);
1047 {
1048 CAPTURE(T);
1049 CAPTURE(c);
1050 CAPTURE(d);
1051 CHECK(check_abs(c, T, acc));
1052 }
1053
1054 c = poly.solve_guess(matrix, 0.0, d, 350, 0);
1055 {
1056 CAPTURE(T);
1057 CAPTURE(c);
1058 CAPTURE(d);
1059 CHECK(check_abs(c, T, acc));
1060 }
1061
1062 // T = 0.0;
1063 // solve.setGuess(75+273.15);
1064 // T = solve.polyval(cHeat,c);
1065 // printf("Should be : T = %3.3f \t K \n",273.15+50.0);
1066 // printf("From object: T = %3.3f \t K \n",T);
1067 //
1068 // T = 0.0;
1069 // solve.setLimits(273.15+10,273.15+100);
1070 // T = solve.polyval(cHeat,c);
1071 // printf("Should be : T = %3.3f \t K \n",273.15+50.0);
1072 // printf("From object: T = %3.3f \t K \n",T);
1073 }
1074
1075 SECTION("Integration and derivation tests") {
1076
1077 CoolProp::Polynomial2D poly;
1078
1079 Eigen::MatrixXd matrix(matrix2D);
1080 Eigen::MatrixXd matrixInt = poly.integrateCoeffs(matrix, 1);
1081 Eigen::MatrixXd matrixDer = poly.deriveCoeffs(matrix, 1);
1082 Eigen::MatrixXd matrixInt2 = poly.integrateCoeffs(matrix, 1, 2);
1083 Eigen::MatrixXd matrixDer2 = poly.deriveCoeffs(matrix, 1, 2);
1084
1085 CHECK_THROWS(poly.evaluate(matrix, 0.0));
1086
1087 double x = 0.3, y = 255.3, val1, val2, val3, val4;
1088
1089 //CHECK( std::abs( polyInt.derivative(x,y,0)-poly2D.evaluate(x,y) ) <= 1e-10 );
1090
1091 std::string tmpStr;
1092
1093 double acc = 0.001;
1094
1095 for (double T = Tmin; T < Tmax; T += Tinc) {
1096 val1 = poly.evaluate(matrix, x, T - deltaT);
1097 val2 = poly.evaluate(matrix, x, T + deltaT);
1098 val3 = (val2 - val1) / 2 / deltaT;
1099 val4 = poly.evaluate(matrixDer, x, T);
1100 CAPTURE(T);
1101 CAPTURE(val3);
1102 CAPTURE(val4);
1103 tmpStr = CoolProp::mat_to_string(matrix);
1104 CAPTURE(tmpStr);
1105 tmpStr = CoolProp::mat_to_string(matrixDer);
1106 CAPTURE(tmpStr);
1107 CHECK(check_abs(val3, val4, acc));
1108 }
1109
1110 for (double T = Tmin; T < Tmax; T += Tinc) {
1111 val1 = poly.evaluate(matrixDer, x, T - deltaT);
1112 val2 = poly.evaluate(matrixDer, x, T + deltaT);
1113 val3 = (val2 - val1) / 2 / deltaT;
1114 val4 = poly.evaluate(matrixDer2, x, T);
1115 CAPTURE(T);
1116 CAPTURE(val3);
1117 CAPTURE(val4);
1118 tmpStr = CoolProp::mat_to_string(matrixDer);
1119 CAPTURE(tmpStr);
1120 tmpStr = CoolProp::mat_to_string(matrixDer2);
1121 CAPTURE(tmpStr);
1122 CHECK(check_abs(val3, val4, acc));
1123 }
1124
1125 for (double T = Tmin; T < Tmax; T += Tinc) {
1126 val1 = poly.evaluate(matrixInt, x, T - deltaT);
1127 val2 = poly.evaluate(matrixInt, x, T + deltaT);
1128 val3 = (val2 - val1) / 2 / deltaT;
1129 val4 = poly.evaluate(matrix, x, T);
1130 CAPTURE(T);
1131 CAPTURE(val3);
1132 CAPTURE(val4);
1133 tmpStr = CoolProp::mat_to_string(matrixInt);
1134 CAPTURE(tmpStr);
1135 tmpStr = CoolProp::mat_to_string(matrix);
1136 CAPTURE(tmpStr);
1137 CHECK(check_abs(val3, val4, acc));
1138 }
1139
1140 for (double T = Tmin; T < Tmax; T += Tinc) {
1141 val1 = poly.evaluate(matrixInt2, x, T - deltaT);
1142 val2 = poly.evaluate(matrixInt2, x, T + deltaT);
1143 val3 = (val2 - val1) / 2 / deltaT;
1144 val4 = poly.evaluate(matrixInt, x, T);
1145 CAPTURE(T);
1146 CAPTURE(val3);
1147 CAPTURE(val4);
1148 tmpStr = CoolProp::mat_to_string(matrixInt2);
1149 CAPTURE(tmpStr);
1150 tmpStr = CoolProp::mat_to_string(matrixInt);
1151 CAPTURE(tmpStr);
1152 CHECK(check_abs(val3, val4, acc));
1153 }
1154
1155 for (double T = Tmin; T < Tmax; T += Tinc) {
1156 val1 = poly.evaluate(matrix, x, T);
1157 val2 = poly.derivative(matrixInt, x, T, 1);
1158 CAPTURE(T);
1159 CAPTURE(val1);
1160 CAPTURE(val2);
1161 CHECK(check_abs(val1, val2, acc));
1162 }
1163
1164 for (double T = Tmin; T < Tmax; T += Tinc) {
1165 val1 = poly.derivative(matrix, x, T, 1);
1166 val2 = poly.evaluate(matrixDer, x, T);
1167 CAPTURE(T);
1168 CAPTURE(val1);
1169 CAPTURE(val2);
1170 CHECK(check_abs(val1, val2, acc));
1171 }
1172 }
1173
1174 SECTION("Testing Polynomial2DFrac") {
1175
1176 Eigen::MatrixXd matrix = CoolProp::vec_to_eigen(cHeat);
1177 CoolProp::Polynomial2D poly;
1178 CoolProp::Polynomial2DFrac frac;
1179
1180 double acc = 0.0001;
1181
1182 double T = 273.15 + 50;
1183 double a, b;
1184 double c = poly.evaluate(matrix, T, 0.0);
1185 double d = frac.evaluate(matrix, T, 0.0, 0, 0);
1186
1187 {
1188 CAPTURE(T);
1189 CAPTURE(c);
1190 CAPTURE(d);
1191 tmpStr = CoolProp::mat_to_string(matrix);
1192 CAPTURE(tmpStr);
1193 CHECK(check_abs(c, d, acc));
1194 }
1195
1196 c = poly.evaluate(matrix, T, 0.0) / T / T;
1197 d = frac.evaluate(matrix, T, 0.0, -2, 0);
1198
1199 {
1200 CAPTURE(T);
1201 CAPTURE(c);
1202 CAPTURE(d);
1203 tmpStr = CoolProp::mat_to_string(matrix);
1204 CAPTURE(tmpStr);
1205 CHECK(check_abs(c, d, acc));
1206 }
1207
1208 matrix = CoolProp::vec_to_eigen(cHeat2D);
1209 double y = 0.1;
1210 c = poly.evaluate(matrix, T, y) / T / T;
1211 d = frac.evaluate(matrix, T, y, -2, 0);
1212
1213 {
1214 CAPTURE(T);
1215 CAPTURE(c);
1216 CAPTURE(d);
1217 tmpStr = CoolProp::mat_to_string(matrix);
1218 CAPTURE(tmpStr);
1219 CHECK(check_abs(c, d, acc));
1220 }
1221
1222 c = poly.evaluate(matrix, T, y) / y / y;
1223 d = frac.evaluate(matrix, T, y, 0, -2);
1224
1225 {
1226 CAPTURE(T);
1227 CAPTURE(c);
1228 CAPTURE(d);
1229 tmpStr = CoolProp::mat_to_string(matrix);
1230 CAPTURE(tmpStr);
1231 CHECK(check_abs(c, d, acc));
1232 }
1233
1234 c = poly.evaluate(matrix, T, y) / T / T / y / y;
1235 d = frac.evaluate(matrix, T, y, -2, -2);
1236
1237 {
1238 CAPTURE(T);
1239 CAPTURE(c);
1240 CAPTURE(d);
1241 tmpStr = CoolProp::mat_to_string(matrix);
1242 CAPTURE(tmpStr);
1243 CHECK(check_abs(c, d, acc));
1244 }
1245
1246 c = poly.evaluate(matrix, T, y) / T / T * y * y;
1247 d = frac.evaluate(matrix, T, y, -2, 2);
1248
1249 {
1250 CAPTURE(T);
1251 CAPTURE(c);
1252 CAPTURE(d);
1253 tmpStr = CoolProp::mat_to_string(matrix);
1254 CAPTURE(tmpStr);
1255 CHECK(check_abs(c, d, acc));
1256 }
1257
1258 matrix = CoolProp::vec_to_eigen(cHeat);
1259 T = 273.15 + 50;
1260 c = 145.59157247249246;
1261 d = frac.integral(matrix, T, 0.0, 0, -1, 0) - frac.integral(matrix, 273.15 + 25, 0.0, 0, -1, 0);
1262
1263 {
1264 CAPTURE(T);
1265 CAPTURE(c);
1266 CAPTURE(d);
1267 tmpStr = CoolProp::mat_to_string(matrix);
1268 CAPTURE(tmpStr);
1269 CHECK(check_abs(c, d, acc));
1270 }
1271
1272 T = 423.15;
1273 c = 3460.895272;
1274 d = frac.integral(matrix, T, 0.0, 0, -1, 0, 348.15, 0.0);
1275
1276 {
1277 CAPTURE(T);
1278 CAPTURE(c);
1279 CAPTURE(d);
1280 tmpStr = CoolProp::mat_to_string(matrix);
1281 CAPTURE(tmpStr);
1282 CHECK(check_abs(c, d, acc));
1283 }
1284
1285 deltaT = 0.01;
1286 for (T = Tmin; T < Tmax; T += Tinc) {
1287 a = poly.evaluate(matrix, T - deltaT, y);
1288 b = poly.evaluate(matrix, T + deltaT, y);
1289 c = (b - a) / 2.0 / deltaT;
1290 d = frac.derivative(matrix, T, y, 0, 0, 0);
1291 CAPTURE(a);
1292 CAPTURE(b);
1293 CAPTURE(T);
1294 CAPTURE(c);
1295 CAPTURE(d);
1296 tmpStr = CoolProp::mat_to_string(matrix);
1297 CAPTURE(tmpStr);
1298 CHECK(check_abs(c, d, acc));
1299 }
1300
1301 T = 273.15 + 150;
1302 c = -2.100108045;
1303 d = frac.derivative(matrix, T, 0.0, 0, 0, 0);
1304 {
1305 CAPTURE(T);
1306 CAPTURE(c);
1307 CAPTURE(d);
1308 tmpStr = CoolProp::mat_to_string(matrix);
1309 CAPTURE(tmpStr);
1310 CHECK(check_abs(c, d, acc));
1311 }
1312
1313 c = -0.006456574589;
1314 d = frac.derivative(matrix, T, 0.0, 0, -1, 0);
1315 {
1316 CAPTURE(T);
1317 CAPTURE(c);
1318 CAPTURE(d);
1319 tmpStr = CoolProp::mat_to_string(matrix);
1320 CAPTURE(tmpStr);
1321 CHECK(check_abs(c, d, acc));
1322 }
1323
1324 c = frac.evaluate(matrix, T, 0.0, 2, 0);
1325 d = frac.solve(matrix, 0.0, c, 0, 2, 0)[0];
1326 {
1327 CAPTURE(T);
1328 CAPTURE(c);
1329 CAPTURE(d);
1330 tmpStr = CoolProp::mat_to_string(matrix);
1331 CAPTURE(tmpStr);
1332 CHECK(check_abs(T, d, acc));
1333 }
1334
1335 c = frac.evaluate(matrix, T, 0.0, 0, 0);
1336 d = frac.solve(matrix, 0.0, c, 0, 0, 0)[0];
1337 {
1338 CAPTURE(T);
1339 CAPTURE(c);
1340 CAPTURE(d);
1341 tmpStr = CoolProp::mat_to_string(matrix);
1342 CAPTURE(tmpStr);
1343 CHECK(check_abs(T, d, acc));
1344 }
1345
1346 c = frac.evaluate(matrix, T, 0.0, -1, 0);
1347 CHECK_THROWS(d = frac.solve(matrix, 0.0, c, 0, -1, 0)[0]);
1348 // {
1349 // CAPTURE(T);
1350 // CAPTURE(c);
1351 // CAPTURE(d);
1352 // tmpStr = CoolProp::mat_to_string(matrix);
1353 // CAPTURE(tmpStr);
1354 // CHECK( check_abs(T,d,acc) );
1355 // }
1356
1357 d = frac.solve_limits(matrix, 0.0, c, T - 10, T + 10, 0, -1, 0);
1358 {
1359 CAPTURE(T);
1360 CAPTURE(c);
1361 CAPTURE(d);
1362 tmpStr = CoolProp::mat_to_string(matrix);
1363 CAPTURE(tmpStr);
1364 CHECK(check_abs(T, d, acc));
1365 }
1366
1367 d = frac.solve_guess(matrix, 0.0, c, T - 10, 0, -1, 0);
1368 {
1369 CAPTURE(T);
1370 CAPTURE(c);
1371 CAPTURE(d);
1372 tmpStr = CoolProp::mat_to_string(matrix);
1373 CAPTURE(tmpStr);
1374 CHECK(check_abs(T, d, acc));
1375 }
1376
1377 c = -0.00004224550082;
1378 d = frac.derivative(matrix, T, 0.0, 0, -2, 0);
1379 {
1380 CAPTURE(T);
1381 CAPTURE(c);
1382 CAPTURE(d);
1383 tmpStr = CoolProp::mat_to_string(matrix);
1384 CAPTURE(tmpStr);
1385 CHECK(check_abs(c, d, acc));
1386 }
1387
1388 c = frac.evaluate(matrix, T, 0.0, 0, 0, 0.0, 0.0);
1389 d = frac.solve(matrix, 0.0, c, 0, 0, 0, 0.0, 0.0)[0];
1390 {
1391 CAPTURE(T);
1392 CAPTURE(c);
1393 CAPTURE(d);
1394 tmpStr = CoolProp::mat_to_string(matrix);
1395 CAPTURE(tmpStr);
1396 tmpStr = CoolProp::mat_to_string(Eigen::MatrixXd(frac.solve(matrix, 0.0, c, 0, 0, 0, 250, 0.0)));
1397 CAPTURE(tmpStr);
1398 CHECK(check_abs(T, d, acc));
1399 }
1400 }
1401}
1402
1403#endif /* ENABLE_CATCH */
Generated on Mon Jun 2 2025 05:21:52 for CoolProp by doxygen 1.9.4