CoolProp 8.0.0
An open-source fluid property and humid air property database
PolyMath.cpp
Go to the documentation of this file.
1
2#include "CoolProp/numerics/PolyMath.h"
3
4#include "CoolProp/Exceptions.h"
5#include "CoolProp/numerics/MatrixMath.h"
6
7#include <cmath>
8#include <vector>
9#include <string>
10//#include <sstream>
11//#include <numeric>
12#include <iostream>
13
14#include "CoolProp/numerics/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
52Eigen::MatrixXd Polynomial2D::integrateCoeffs(const Eigen::MatrixXd& coefficients, const int& axis = -1, const int& times = 1) {
53 if (times < 0)
54 throw ValueError(format("%s (%d): You have to provide a positive order for integration, %d is not valid. ", __FILE__, __LINE__, times));
55 if (times == 0) return Eigen::MatrixXd(coefficients);
56 Eigen::MatrixXd oldCoefficients;
57 Eigen::MatrixXd newCoefficients(coefficients);
58
59 switch (axis) {
60 case 0:
61 newCoefficients = Eigen::MatrixXd(coefficients);
62 break;
63 case 1:
64 newCoefficients = Eigen::MatrixXd(coefficients.transpose());
65 break;
66 default:
67 throw ValueError(
68 format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ", __FILE__, __LINE__, axis));
69 break;
70 }
71
72 std::size_t r = 0, c = 0;
73 for (int k = 0; k < times; k++) {
74 oldCoefficients = Eigen::MatrixXd(newCoefficients);
75 r = oldCoefficients.rows(), c = oldCoefficients.cols();
76 newCoefficients = Eigen::MatrixXd::Zero(r + 1, c);
77 newCoefficients.block(1, 0, r, c) = oldCoefficients.block(0, 0, r, c);
78 for (size_t i = 0; i < r; ++i) {
79 for (size_t j = 0; j < c; ++j)
80 newCoefficients(i + 1, j) /= (i + 1.);
81 }
82 }
83
84 switch (axis) {
85 case 0:
86 break;
87 case 1:
88 newCoefficients.transposeInPlace();
89 break;
90 default:
91 throw ValueError(
92 format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ", __FILE__, __LINE__, axis));
93 break;
94 }
95
96 return newCoefficients;
97}
98
100
103Eigen::MatrixXd Polynomial2D::deriveCoeffs(const Eigen::MatrixXd& coefficients, const int& axis, const int& times) {
104 if (times < 0)
105 throw ValueError(format("%s (%d): You have to provide a positive order for derivation, %d is not valid. ", __FILE__, __LINE__, times));
106 if (times == 0) return Eigen::MatrixXd(coefficients);
107 // Recursion is also possible, but not recommended
108 //Eigen::MatrixXd newCoefficients;
109 //if (times > 1) newCoefficients = deriveCoeffs(coefficients, axis, times-1);
110 //else newCoefficients = Eigen::MatrixXd(coefficients);
111 Eigen::MatrixXd newCoefficients;
112
113 switch (axis) {
114 case 0:
115 newCoefficients = Eigen::MatrixXd(coefficients);
116 break;
117 case 1:
118 newCoefficients = Eigen::MatrixXd(coefficients.transpose());
119 break;
120 default:
121 throw ValueError(
122 format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ", __FILE__, __LINE__, axis));
123 break;
124 }
125
126 std::size_t r = 0, c = 0, i = 0, j = 0;
127 for (int k = 0; k < times; k++) {
128 r = newCoefficients.rows(), c = newCoefficients.cols();
129 for (i = 1; i < r; ++i) {
130 for (j = 0; j < c; ++j)
131 newCoefficients(i, j) *= i;
132 }
133 removeRow(newCoefficients, 0);
134 }
135
136 switch (axis) {
137 case 0:
138 break;
139 case 1:
140 newCoefficients.transposeInPlace();
141 break;
142 default:
143 throw ValueError(
144 format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ", __FILE__, __LINE__, axis));
145 break;
146 }
147
148 return newCoefficients;
149}
150
152
159double Polynomial2D::evaluate(const Eigen::MatrixXd& coefficients, const double& x_in) {
160 double result = Eigen::poly_eval(makeVector(coefficients), x_in);
161 if (this->do_debug())
162 std::cout << "Running 1D evaluate(" << mat_to_string(coefficients) << ", x_in:" << vec_to_string(x_in) << "): " << result << '\n';
163 return result;
164}
165double Polynomial2D::evaluate(const Eigen::MatrixXd& coefficients, const double& x_in, const double& y_in) {
166 size_t r = coefficients.rows();
167 double result = evaluate(coefficients.row(r - 1), y_in);
168 for (int i = static_cast<int>(r) - 2; i >= 0; i--) {
169 result *= x_in;
170 result += evaluate(coefficients.row(i), y_in);
171 }
172 if (this->do_debug())
173 std::cout << "Running 2D evaluate(" << mat_to_string(coefficients) << ", x_in:" << vec_to_string(x_in)
174 << ", y_in:" << vec_to_string(y_in) << "): " << result << '\n';
175 return result;
176}
177
178double Polynomial2D::derivative(const Eigen::MatrixXd& coefficients, const double& x_in, const double& y_in, const int& axis = -1) {
179 return this->evaluate(this->deriveCoeffs(coefficients, axis, 1), x_in, y_in);
180}
181
182double Polynomial2D::integral(const Eigen::MatrixXd& coefficients, const double& x_in, const double& y_in, const int& axis = -1) {
183 return this->evaluate(this->integrateCoeffs(coefficients, axis, 1), x_in, y_in);
184}
185
190double Polynomial2D::solve_limits(Poly2DResidual* res, const double& min, const double& max) {
191 if (do_debug()) std::cout << format("Called solve_limits with: min=%f and max=%f", min, max) << '\n';
192 double macheps = DBL_EPSILON;
193 double tol = DBL_EPSILON * 1e3;
194 int maxiter = 10;
195 double result = Brent(res, min, max, macheps, tol, maxiter);
196 if (this->do_debug()) std::cout << "Brent solver message: " << res->errstring << '\n';
197 return result;
198}
199
203double Polynomial2D::solve_guess(Poly2DResidual* res, const double& guess) {
204 if (do_debug()) std::cout << format("Called solve_guess with: guess=%f ", guess) << '\n';
205 //set_debug_level(1000);
206 double tol = DBL_EPSILON * 1e3;
207 int maxiter = 10;
208 double result = Newton(res, guess, tol, maxiter);
209 if (this->do_debug()) std::cout << "Newton solver message: " << res->errstring << '\n';
210 return result;
211}
212
213Eigen::VectorXd Polynomial2D::solve(const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const int& axis = -1) {
214 std::size_t r = coefficients.rows(), c = coefficients.cols();
215 Eigen::VectorXd tmpCoefficients;
216 switch (axis) {
217 case 0:
218 tmpCoefficients = Eigen::VectorXd::Zero(r);
219 for (size_t i = 0; i < r; i++) {
220 tmpCoefficients(i, 0) = evaluate(coefficients.row(i), in);
221 }
222 break;
223 case 1:
224 tmpCoefficients = Eigen::VectorXd::Zero(c);
225 for (size_t i = 0; i < c; i++) {
226 tmpCoefficients(i, 0) = evaluate(coefficients.col(i), in);
227 }
228 break;
229 default:
230 throw ValueError(format("%s (%d): You have to provide a dimension, 0 or 1, for the solver, %d is not valid. ", __FILE__, __LINE__, axis));
231 break;
232 }
233 tmpCoefficients(0, 0) -= z_in;
234 if (this->do_debug()) std::cout << "Coefficients: " << mat_to_string(Eigen::MatrixXd(tmpCoefficients)) << '\n';
235 Eigen::PolynomialSolver<double, Eigen::Dynamic> polySolver(tmpCoefficients);
236 std::vector<double> rootsVec;
237 polySolver.realRoots(rootsVec);
238 if (this->do_debug()) std::cout << "Real roots: " << vec_to_string(rootsVec) << '\n';
239 return vec_to_eigen(rootsVec);
240 //return rootsVec[0]; // TODO: implement root selection algorithm
241}
242
243double Polynomial2D::solve_limits(const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const double& min, const double& max,
244 const int& axis) {
245 Poly2DResidual res = Poly2DResidual(*this, coefficients, in, z_in, axis);
246 return solve_limits(&res, min, max);
247}
248
249double Polynomial2D::solve_guess(const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const double& guess, const int& axis) {
250 Poly2DResidual res = Poly2DResidual(*this, coefficients, in, z_in, axis);
251 return solve_guess(&res, guess);
252}
253
255
258double Polynomial2D::simplePolynomial(std::vector<double> const& coefficients, double x) {
259 double result = 0.;
260 for (unsigned int i = 0; i < coefficients.size(); i++) {
261 result += coefficients[i] * pow(x, (int)i);
262 }
263 if (this->do_debug())
264 std::cout << "Running simplePolynomial(" << vec_to_string(coefficients) << ", " << vec_to_string(x) << "): " << result << '\n';
265 return result;
266}
267double Polynomial2D::simplePolynomial(std::vector<std::vector<double>> const& coefficients, double x, double y) {
268 double result = 0;
269 for (unsigned int i = 0; i < coefficients.size(); i++) {
270 result += pow(x, (int)i) * simplePolynomial(coefficients[i], y);
271 }
272 if (this->do_debug())
273 std::cout << "Running simplePolynomial(" << vec_to_string(coefficients) << ", " << vec_to_string(x) << ", " << vec_to_string(y)
274 << "): " << result << '\n';
275 return result;
276}
277
279
283double Polynomial2D::baseHorner(std::vector<double> const& coefficients, double x) {
284 double result = 0;
285 for (int i = static_cast<int>(coefficients.size()) - 1; i >= 0; i--) {
286 result *= x;
287 result += coefficients[i];
288 }
289 if (this->do_debug())
290 std::cout << "Running baseHorner(" << vec_to_string(coefficients) << ", " << vec_to_string(x) << "): " << result << '\n';
291 return result;
292}
293
294double Polynomial2D::baseHorner(std::vector<std::vector<double>> const& coefficients, double x, double y) {
295 double result = 0;
296 for (int i = static_cast<int>(coefficients.size() - 1); i >= 0; i--) {
297 result *= x;
298 result += baseHorner(coefficients[i], y);
299 }
300 if (this->do_debug())
301 std::cout << "Running baseHorner(" << vec_to_string(coefficients) << ", " << vec_to_string(x) << ", " << vec_to_string(y)
302 << "): " << result << '\n';
303 return result;
304}
305
306Poly2DResidual::Poly2DResidual(Polynomial2D& poly, const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const int& axis) {
307 switch (axis) {
308 case iX:
309 case iY:
310 this->axis = axis;
311 break;
312 default:
313 throw ValueError(format("%s (%d): You have to provide a dimension to the solver, %d is not valid. ", __FILE__, __LINE__, axis));
314 break;
315 }
316
317 this->poly = poly;
318 this->coefficients = coefficients;
319 this->derIsSet = false;
320 this->in = in;
321 this->z_in = z_in;
322}
323
324double Poly2DResidual::call(double target) {
325 if (axis == iX) return poly.evaluate(coefficients, target, in) - z_in;
326 if (axis == iY) return poly.evaluate(coefficients, in, target) - z_in;
327 return -_HUGE;
328}
329
330double Poly2DResidual::deriv(double target) {
331 if (!this->derIsSet) {
332 this->coefficientsDer = poly.deriveCoeffs(coefficients, axis);
333 this->derIsSet = true;
334 }
335 return poly.evaluate(coefficientsDer, target, in);
336}
337
338// /// Integration functions
339// /** Integrating coefficients for polynomials is done by dividing the
340// * original coefficients by (i+1) and elevating the order by 1
341// * through adding a zero as first coefficient.
342// * Some reslicing needs to be applied to integrate along the x-axis.
343// * In the brine/solution equations, reordering of the parameters
344// * avoids this expensive operation. However, it is included for the
345// * sake of completeness.
346// */
347// /// @param coefficients matrix containing the ordered coefficients
348// /// @param axis unsigned integer value that represents the desired direction of integration
349// /// @param times integer value that represents the desired order of integration
350// /// @param firstExponent integer value that represents the first exponent of the polynomial in axis direction
351// Eigen::MatrixXd integrateCoeffs(const Eigen::MatrixXd &coefficients, const int &axis, const int &times, const int &firstExponent);
352//
354
363Eigen::MatrixXd Polynomial2DFrac::deriveCoeffs(const Eigen::MatrixXd& coefficients, const int& axis, const int& times, const int& firstExponent) {
364 if (times < 0)
365 throw ValueError(format("%s (%d): You have to provide a positive order for derivation, %d is not valid. ", __FILE__, __LINE__, times));
366 if (times == 0) return Eigen::MatrixXd(coefficients);
367 // Recursion is also possible, but not recommended
368 //Eigen::MatrixXd newCoefficients;
369 //if (times > 1) newCoefficients = deriveCoeffs(coefficients, axis, times-1);
370 //else newCoefficients = Eigen::MatrixXd(coefficients);
371 Eigen::MatrixXd newCoefficients;
372
373 switch (axis) {
374 case 0:
375 newCoefficients = Eigen::MatrixXd(coefficients);
376 break;
377 case 1:
378 newCoefficients = Eigen::MatrixXd(coefficients.transpose());
379 break;
380 default:
381 throw ValueError(
382 format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ", __FILE__, __LINE__, axis));
383 break;
384 }
385
386 std::size_t r = newCoefficients.rows(), c = newCoefficients.cols();
387 std::size_t i = 0, j = 0;
388 for (int k = 0; k < times; k++) {
389 for (i = 0; i < r; ++i) {
390 for (j = 0; j < c; ++j) {
391 newCoefficients(i, j) *= double(i) + double(firstExponent);
392 }
393 }
394 }
395
396 switch (axis) {
397 case 0:
398 break;
399 case 1:
400 newCoefficients.transposeInPlace();
401 break;
402 default:
403 throw ValueError(
404 format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ", __FILE__, __LINE__, axis));
405 break;
406 }
407
408 return newCoefficients;
409}
410
412
422double Polynomial2DFrac::evaluate(const Eigen::MatrixXd& coefficients, const double& x_in, const int& firstExponent, const double& x_base) {
423 size_t r = coefficients.rows();
424 size_t c = coefficients.cols();
425
426 if ((r != 1) && (c != 1)) {
427 throw ValueError(format("%s (%d): You have a 2D coefficient matrix (%d,%d), please use the 2D functions. ", __FILE__, __LINE__,
428 coefficients.rows(), coefficients.cols()));
429 }
430 if ((firstExponent < 0) && (std::abs(x_in - x_base) < CPPOLY_EPSILON)) {
431 //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));
432 const double x_lo = x_base - CPPOLY_DELTA;
433 const double x_hi = x_base + CPPOLY_DELTA;
434 const double y_lo = evaluate(coefficients, x_lo, firstExponent, x_base);
435 const double y_hi = evaluate(coefficients, x_hi, firstExponent, x_base);
436 return (y_hi - y_lo) / (x_hi - x_lo) * (x_in - x_lo) + y_lo;
437 }
438
439 Eigen::MatrixXd tmpCoeffs(coefficients);
440 if (c == 1) {
441 tmpCoeffs.transposeInPlace();
442 c = r;
443 r = 1;
444 }
445 Eigen::MatrixXd newCoeffs;
446 double negExp = 0; // First we treat the negative exponents
447 double posExp = 0; // then the positive exponents
448
449 for (int i = 0; i > firstExponent; i--) { // only for firstExponent<0
450 if (c > 0) {
451 negExp += tmpCoeffs(0, 0);
452 removeColumn(tmpCoeffs, 0);
453 }
454 negExp /= x_in - x_base;
455 c = tmpCoeffs.cols();
456 }
457
458 for (int i = 0; i < firstExponent; i++) { // only for firstExponent>0
459 newCoeffs = Eigen::MatrixXd::Zero(r, c + 1);
460 newCoeffs.block(0, 1, r, c) = tmpCoeffs.block(0, 0, r, c);
461 tmpCoeffs = Eigen::MatrixXd(newCoeffs);
462 c = tmpCoeffs.cols();
463 }
464
465 if (c > 0) posExp += Eigen::poly_eval(Eigen::RowVectorXd(tmpCoeffs), x_in - x_base);
466 return negExp + posExp;
467}
468
469double Polynomial2DFrac::evaluate(const Eigen::MatrixXd& coefficients, const double& x_in, const double& y_in, const int& x_exp, const int& y_exp,
470 const double& x_base, const double& y_base) {
471 if ((x_exp < 0) && (std::abs(x_in - x_base) < CPPOLY_EPSILON)) {
472 // 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));
473 if (this->do_debug()) std::cout << "Interpolating in x-direction for base " << x_base << " and input " << x_in << '\n';
474 const double x_lo = x_base - CPPOLY_DELTA;
475 const double x_hi = x_base + CPPOLY_DELTA;
476 const double z_lo = evaluate(coefficients, x_lo, y_in, x_exp, y_exp, x_base, y_base);
477 const double z_hi = evaluate(coefficients, x_hi, y_in, x_exp, y_exp, x_base, y_base);
478 return (z_hi - z_lo) / (x_hi - x_lo) * (x_in - x_lo) + z_lo;
479 }
480 if ((y_exp < 0) && (std::abs(y_in - y_base) < CPPOLY_EPSILON)) {
481 // 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));
482 if (this->do_debug()) std::cout << "Interpolating in y-direction for base " << y_base << " and input " << y_in << '\n';
483 const double y_lo = y_base - CPPOLY_DELTA;
484 const double y_hi = y_base + CPPOLY_DELTA;
485 const double z_lo = evaluate(coefficients, x_in, y_lo, x_exp, y_exp, x_base, y_base);
486 const double z_hi = evaluate(coefficients, x_in, y_hi, x_exp, y_exp, x_base, y_base);
487 return (z_hi - z_lo) / (y_hi - y_lo) * (y_in - y_lo) + z_lo;
488 }
489
490 Eigen::MatrixXd tmpCoeffs(coefficients);
491 Eigen::MatrixXd newCoeffs;
492 size_t r = 0;
493 size_t c = tmpCoeffs.cols();
494 double negExp = 0; // First we treat the negative exponents
495 double posExp = 0; // then the positive exponents
496
497 for (int i = 0; i > x_exp; i--) { // only for x_exp<0
498 r = tmpCoeffs.rows();
499 if (r > 0) {
500 negExp += evaluate(tmpCoeffs.row(0), y_in, y_exp, y_base);
501 removeRow(tmpCoeffs, 0);
502 }
503 negExp /= x_in - x_base;
504 }
505
506 r = tmpCoeffs.rows();
507 for (int i = 0; i < x_exp; i++) { // only for x_exp>0
508 newCoeffs = Eigen::MatrixXd::Zero(r + 1, c);
509 newCoeffs.block(1, 0, r, c) = tmpCoeffs.block(0, 0, r, c);
510 tmpCoeffs = Eigen::MatrixXd(newCoeffs);
511 r += 1; // r = tmpCoeffs.rows();
512 }
513
514 //r = tmpCoeffs.rows();
515 if (r > 0) posExp += evaluate(tmpCoeffs.row(r - 1), y_in, y_exp, y_base);
516 for (int i = static_cast<int>(r) - 2; i >= 0; i--) {
517 posExp *= x_in - x_base;
518 posExp += evaluate(tmpCoeffs.row(i), y_in, y_exp, y_base);
519 }
520 if (this->do_debug()) std::cout << "Running 2D evaluate(" << mat_to_string(coefficients) << ", " << '\n';
521 if (this->do_debug())
522 std::cout << "x_in:" << vec_to_string(x_in) << ", y_in:" << vec_to_string(y_in) << ", x_exp:" << vec_to_string(x_exp)
523 << ", y_exp:" << vec_to_string(y_exp) << ", x_base:" << vec_to_string(x_base) << ", y_base:" << vec_to_string(y_base)
524 << "): " << negExp + posExp << '\n';
525 return negExp + posExp;
526}
527
528double Polynomial2DFrac::derivative(const Eigen::MatrixXd& coefficients, const double& x_in, const double& y_in, const int& axis, const int& x_exp,
529 const int& y_exp, const double& x_base, const double& y_base) {
530 Eigen::MatrixXd newCoefficients;
531 int der_exp = 0, other_exp = 0;
532 double der_val = NAN, other_val = NAN;
533 double int_base = NAN, other_base = NAN;
534
535 switch (axis) {
536 case 0:
537 newCoefficients = Eigen::MatrixXd(coefficients);
538 der_exp = x_exp;
539 other_exp = y_exp;
540 der_val = x_in;
541 other_val = y_in;
542 int_base = x_base;
543 other_base = y_base;
544 break;
545 case 1:
546 newCoefficients = Eigen::MatrixXd(coefficients.transpose());
547 der_exp = y_exp;
548 other_exp = x_exp;
549 der_val = y_in;
550 other_val = x_in;
551 int_base = y_base;
552 other_base = x_base;
553 break;
554 default:
555 throw ValueError(
556 format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ", __FILE__, __LINE__, axis));
557 break;
558 }
559
560 const int times = 1;
561 newCoefficients = deriveCoeffs(newCoefficients, 0, times, der_exp);
562 der_exp -= times;
563
564 return evaluate(newCoefficients, der_val, other_val, der_exp, other_exp, int_base, other_base);
565}
566
567double Polynomial2DFrac::integral(const Eigen::MatrixXd& coefficients, const double& x_in, const double& y_in, const int& axis, const int& x_exp,
568 const int& y_exp, const double& x_base, const double& y_base, const double& ax_val) {
569
570 Eigen::MatrixXd newCoefficients;
571 int int_exp = 0, other_exp = 0;
572 double int_val = NAN, other_val = NAN;
573 double int_base = NAN, other_base = NAN;
574
575 switch (axis) {
576 case 0:
577 newCoefficients = Eigen::MatrixXd(coefficients);
578 int_exp = x_exp;
579 other_exp = y_exp;
580 int_val = x_in;
581 other_val = y_in;
582 int_base = x_base;
583 other_base = y_base;
584 break;
585 case 1:
586 newCoefficients = Eigen::MatrixXd(coefficients.transpose());
587 int_exp = y_exp;
588 other_exp = x_exp;
589 int_val = y_in;
590 other_val = x_in;
591 int_base = y_base;
592 other_base = x_base;
593 break;
594 default:
595 throw ValueError(
596 format("%s (%d): You have to provide a dimension, 0 or 1, for integration, %d is not valid. ", __FILE__, __LINE__, axis));
597 break;
598 }
599
600 if (int_exp < -1)
601 throw NotImplementedError(
602 format("%s (%d): This function is only implemented for lowest exponents >= -1, int_exp=%d ", __FILE__, __LINE__, int_exp));
603 // TODO: Fix this and allow the direct calculation of integrals
604 if (std::abs(ax_val) > DBL_EPSILON)
605 throw NotImplementedError(
606 format("%s (%d): This function is only implemented for indefinite integrals, ax_val=%d ", __FILE__, __LINE__, ax_val));
607
608 size_t r = newCoefficients.rows();
609 size_t c = newCoefficients.cols();
610
611 if (int_exp == -1) {
612 if (std::abs(int_base) < DBL_EPSILON) {
613 Eigen::MatrixXd tmpCoefficients = newCoefficients.row(0) * log(int_val - int_base);
614 newCoefficients = integrateCoeffs(newCoefficients.block(1, 0, r - 1, c), 0, 1);
615 newCoefficients.row(0) = tmpCoefficients;
616 return evaluate(newCoefficients, int_val, other_val, 0, other_exp, int_base, other_base);
617 } else {
618 // Reduce the coefficients to the integration dimension:
619 newCoefficients = Eigen::MatrixXd(r, 1);
620 for (std::size_t i = 0; i < r; i++) {
621 newCoefficients(i, 0) = evaluate(coefficients.row(i), other_val, other_exp, other_base);
622 }
623 return fracIntCentral(newCoefficients.transpose(), int_val, int_base);
624 }
625 }
626
627 Eigen::MatrixXd tmpCoeffs;
628 r = newCoefficients.rows();
629 for (int i = 0; i < int_exp; i++) { // only for x_exp>0
630 tmpCoeffs = Eigen::MatrixXd::Zero(r + 1, c);
631 tmpCoeffs.block(1, 0, r, c) = newCoefficients.block(0, 0, r, c);
632 newCoefficients = Eigen::MatrixXd(tmpCoeffs);
633 r += 1; // r = newCoefficients.rows();
634 }
635
636 return evaluate(integrateCoeffs(newCoefficients, 0, 1), int_val, other_val, 0, other_exp, int_base, other_base);
637}
638
640Eigen::VectorXd Polynomial2DFrac::solve(const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const int& axis, const int& x_exp,
641 const int& y_exp, const double& x_base, const double& y_base) {
642
643 Eigen::MatrixXd newCoefficients;
644 Eigen::VectorXd tmpCoefficients;
645 int solve_exp = 0, other_exp = 0;
646 double input = NAN;
647
648 switch (axis) {
649 case 0:
650 newCoefficients = Eigen::MatrixXd(coefficients);
651 solve_exp = x_exp;
652 other_exp = y_exp;
653 input = in - y_base;
654 break;
655 case 1:
656 newCoefficients = Eigen::MatrixXd(coefficients.transpose());
657 solve_exp = y_exp;
658 other_exp = x_exp;
659 input = in - x_base;
660 break;
661 default:
662 throw ValueError(format("%s (%d): You have to provide a dimension, 0 or 1, for the solver, %d is not valid. ", __FILE__, __LINE__, axis));
663 break;
664 }
665
666 if (this->do_debug()) std::cout << "Coefficients: " << mat_to_string(Eigen::MatrixXd(newCoefficients)) << '\n';
667
668 const size_t r = newCoefficients.rows();
669 for (size_t i = 0; i < r; i++) {
670 newCoefficients(i, 0) = evaluate(newCoefficients.row(i), input, other_exp);
671 }
672
673 //Eigen::VectorXd tmpCoefficients;
674 if (solve_exp >= 0) { // extend vector to zero exponent
675 tmpCoefficients = Eigen::VectorXd::Zero(r + solve_exp);
676 tmpCoefficients.block(solve_exp, 0, r, 1) = newCoefficients.block(0, 0, r, 1);
677 tmpCoefficients(0, 0) -= z_in;
678 } else { // check if vector reaches to zero exponent
679 int diff = 1 - static_cast<int>(r) - solve_exp; // How many entries have to be added
680 tmpCoefficients = Eigen::VectorXd::Zero(r + std::max(diff, 0));
681 tmpCoefficients.block(0, 0, r, 1) = newCoefficients.block(0, 0, r, 1);
682 tmpCoefficients(r + diff - 1, 0) -= z_in;
683 throw NotImplementedError(format("%s (%d): Currently, there is no solver for the generalised polynomial, an exponent of %d is not valid. ",
684 __FILE__, __LINE__, solve_exp));
685 }
686
687 if (this->do_debug()) std::cout << "Coefficients: " << mat_to_string(Eigen::MatrixXd(tmpCoefficients)) << '\n';
688 Eigen::PolynomialSolver<double, -1> polySolver(tmpCoefficients);
689 std::vector<double> rootsVec;
690 polySolver.realRoots(rootsVec);
691 if (this->do_debug()) std::cout << "Real roots: " << vec_to_string(rootsVec) << '\n';
692 return vec_to_eigen(rootsVec);
693 //return rootsVec[0]; // TODO: implement root selection algorithm
694}
695
697double Polynomial2DFrac::solve_limits(const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const double& min, const double& max,
698 const int& axis, const int& x_exp, const int& y_exp, const double& x_base, const double& y_base) {
699 if (do_debug())
700 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)
701 << '\n';
702 Poly2DFracResidual res = Poly2DFracResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base);
703 return Polynomial2D::solve_limits(&res, min, max);
704} //TODO: Implement tests for this solver
705
707double Polynomial2DFrac::solve_guess(const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const double& guess, const int& axis,
708 const int& x_exp, const int& y_exp, const double& x_base, const double& y_base) {
709 if (do_debug())
710 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) << '\n';
711 Poly2DFracResidual res = Poly2DFracResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base);
712 return Polynomial2D::solve_guess(&res, guess);
713} //TODO: Implement tests for this solver
714
716double Polynomial2DFrac::solve_limitsInt(const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const double& min,
717 const double& max, const int& axis, const int& x_exp, const int& y_exp, const double& x_base,
718 const double& y_base, const int& int_axis) {
719 Poly2DFracIntResidual res = Poly2DFracIntResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base, int_axis);
720 return Polynomial2D::solve_limits(&res, min, max);
721} //TODO: Implement tests for this solver
722
724double Polynomial2DFrac::solve_guessInt(const Eigen::MatrixXd& coefficients, const double& in, const double& z_in, const double& guess,
725 const int& axis, const int& x_exp, const int& y_exp, const double& x_base, const double& y_base,
726 const int& int_axis) {
727 Poly2DFracIntResidual res = Poly2DFracIntResidual(*this, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base, int_axis);
728 return Polynomial2D::solve_guess(&res, guess);
729} //TODO: Implement tests for this solver
730
739//Helper functions to calculate binomial coefficients:
740//http://rosettacode.org/wiki/Evaluate_binomial_coefficients#C.2B.2B
741double Polynomial2DFrac::factorial(const int& nValue) {
742 double value = 1;
743 for (int i = 2; i <= nValue; i++)
744 value = value * i;
745 return value;
746}
749double Polynomial2DFrac::binom(const int& nValue, const int& nValue2) {
750 if (nValue2 == 1) return nValue * 1.0;
751 return (factorial(nValue)) / (factorial(nValue2) * factorial((nValue - nValue2)));
752}
757Eigen::MatrixXd Polynomial2DFrac::fracIntCentralDvector(const int& m, const double& x_in, const double& x_base) {
758 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));
759
760 Eigen::MatrixXd D = Eigen::MatrixXd::Zero(1, m);
761 double tmp = NAN;
762 // TODO: This can be optimized using the Horner scheme!
763 for (int j = 0; j < m; j++) { // loop through row
764 tmp = pow(-1.0, j) * log(x_in) * pow(x_base, j);
765 for (int k = 0; k < j; k++) { // internal loop for every entry
766 tmp += binom(j, k) * pow(-1.0, k) / (j - k) * pow(x_in, j - k) * pow(x_base, k);
767 }
768 D(0, j) = tmp;
769 }
770 return D;
771}
776double Polynomial2DFrac::fracIntCentral(const Eigen::MatrixXd& coefficients, const double& x_in, const double& x_base) {
777 if (coefficients.rows() != 1) {
778 throw ValueError(format("%s (%d): You have a 2D coefficient matrix (%d,%d), please use the 2D functions. ", __FILE__, __LINE__,
779 coefficients.rows(), coefficients.cols()));
780 }
781 int m = static_cast<int>(coefficients.cols());
782 Eigen::MatrixXd D = fracIntCentralDvector(m, x_in, x_base);
783 double result = 0;
784 for (int j = 0; j < m; j++) {
785 result += coefficients(0, j) * D(0, j);
786 }
787 if (this->do_debug())
788 std::cout << "Running fracIntCentral(" << mat_to_string(coefficients) << ", " << vec_to_string(x_in) << ", " << vec_to_string(x_base)
789 << "): " << result << '\n';
790 return result;
791}
792
793Poly2DFracResidual::Poly2DFracResidual(Polynomial2DFrac& poly, const Eigen::MatrixXd& coefficients, const double& in, const double& z_in,
794 const int& axis, const int& x_exp, const int& y_exp, const double& x_base, const double& y_base)
795 : Poly2DResidual(poly, coefficients, in, z_in, axis), x_exp(x_exp), y_exp(y_exp), x_base(x_base), y_base(y_base) {}
796
797double Poly2DFracResidual::call(double target) {
798 if (axis == iX) return poly.evaluate(coefficients, target, in, x_exp, y_exp, x_base, y_base) - z_in;
799 if (axis == iY) return poly.evaluate(coefficients, in, target, x_exp, y_exp, x_base, y_base) - z_in;
800 return _HUGE;
801}
802
803double Poly2DFracResidual::deriv(double target) {
804 if (axis == iX) return poly.derivative(coefficients, target, in, axis, x_exp, y_exp, x_base, y_base);
805 if (axis == iY) return poly.derivative(coefficients, in, target, axis, x_exp, y_exp, x_base, y_base);
806 return _HUGE;
807}
808
809Poly2DFracIntResidual::Poly2DFracIntResidual(Polynomial2DFrac& poly, const Eigen::MatrixXd& coefficients, const double& in, const double& z_in,
810 const int& axis, const int& x_exp, const int& y_exp, const double& x_base, const double& y_base,
811 const int& int_axis)
812 : Poly2DFracResidual(poly, coefficients, in, z_in, axis, x_exp, y_exp, x_base, y_base), int_axis(int_axis) {}
813
814double Poly2DFracIntResidual::call(double target) {
815 if (axis == iX) return poly.integral(coefficients, target, in, int_axis, x_exp, y_exp, x_base, y_base) - z_in;
816 if (axis == iY) return poly.integral(coefficients, in, target, int_axis, x_exp, y_exp, x_base, y_base) - z_in;
817 return _HUGE;
818}
819
820double Poly2DFracIntResidual::deriv(double target) {
821 if (axis == iX) return poly.evaluate(coefficients, target, in, x_exp, y_exp, x_base, y_base);
822 if (axis == iY) return poly.evaluate(coefficients, in, target, x_exp, y_exp, x_base, y_base);
823 return _HUGE;
824}
825
826} // namespace CoolProp
827
828#ifdef ENABLE_CATCH
829# include <math.h>
830# include <iostream>
831# include <catch2/catch_all.hpp>
832
833TEST_CASE("Internal consistency checks and example use cases for PolyMath.cpp", "[PolyMath]") {
834 bool PRINT = false;
835 std::string tmpStr;
836
838 std::vector<double> cHeat;
839 cHeat.clear();
840 cHeat.push_back(+1.1562261074E+03);
841 cHeat.push_back(+2.0994549103E+00);
842 cHeat.push_back(+7.7175381057E-07);
843 cHeat.push_back(-3.7008444051E-20);
844
845 double deltaT = 0.1;
846 double Tmin = 273.15 - 50;
847 double Tmax = 273.15 + 250;
848 double Tinc = 200;
849
850 std::vector<std::vector<double>> cHeat2D;
851 cHeat2D.push_back(cHeat);
852 cHeat2D.push_back(cHeat);
853 cHeat2D.push_back(cHeat);
854
855 Eigen::MatrixXd matrix2D = CoolProp::vec_to_eigen(cHeat2D);
856
857 Eigen::MatrixXd matrix2Dtmp;
858 std::vector<std::vector<double>> vec2Dtmp;
859
860 SECTION("Coefficient parsing") {
861 CoolProp::Polynomial2D poly;
862 CHECK_THROWS(poly.checkCoefficients(matrix2D, 4, 5));
863 CHECK(poly.checkCoefficients(matrix2D, 3, 4));
864 }
865
866 SECTION("Coefficient operations") {
867 Eigen::MatrixXd matrix;
868 CoolProp::Polynomial2D poly;
869
870 CHECK_THROWS(poly.integrateCoeffs(matrix2D));
871
872 CHECK_NOTHROW(matrix = poly.integrateCoeffs(matrix2D, 0));
873 tmpStr = CoolProp::mat_to_string(matrix2D);
874 if (PRINT) std::cout << tmpStr << '\n';
875 tmpStr = CoolProp::mat_to_string(matrix);
876 if (PRINT) std::cout << tmpStr << '\n' << '\n';
877
878 CHECK_NOTHROW(matrix = poly.integrateCoeffs(matrix2D, 1));
879 tmpStr = CoolProp::mat_to_string(matrix2D);
880 if (PRINT) std::cout << tmpStr << '\n';
881 tmpStr = CoolProp::mat_to_string(matrix);
882 if (PRINT) std::cout << tmpStr << '\n' << '\n';
883
884 CHECK_THROWS(poly.deriveCoeffs(matrix2D));
885
886 CHECK_NOTHROW(matrix = poly.deriveCoeffs(matrix2D, 0));
887 tmpStr = CoolProp::mat_to_string(matrix2D);
888 if (PRINT) std::cout << tmpStr << '\n';
889 tmpStr = CoolProp::mat_to_string(matrix);
890 if (PRINT) std::cout << tmpStr << '\n' << '\n';
891
892 CHECK_NOTHROW(matrix = poly.deriveCoeffs(matrix2D, 1));
893 tmpStr = CoolProp::mat_to_string(matrix2D);
894 if (PRINT) std::cout << tmpStr << '\n';
895 tmpStr = CoolProp::mat_to_string(matrix);
896 if (PRINT) std::cout << tmpStr << '\n' << '\n';
897 }
898
899 SECTION("Evaluation and test values") {
900
901 Eigen::MatrixXd matrix = CoolProp::vec_to_eigen(cHeat);
902 CoolProp::Polynomial2D poly;
903
904 double acc = 0.0001;
905
906 double T = 273.15 + 50;
907 double c = poly.evaluate(matrix, T, 0.0);
908 double d = 1834.746;
909
910 {
911 CAPTURE(T);
912 CAPTURE(c);
913 CAPTURE(d);
914 tmpStr = CoolProp::mat_to_string(matrix);
915 CAPTURE(tmpStr);
916 CHECK(check_abs(c, d, acc));
917 }
918
919 c = poly.solve(matrix, 0.0, d, 0)[0];
920 {
921 CAPTURE(T);
922 CAPTURE(c);
923 CAPTURE(d);
924 CHECK(check_abs(c, T, acc));
925 }
926
927 c = poly.solve_limits(matrix, 0.0, d, -50, 750, 0);
928 {
929 CAPTURE(T);
930 CAPTURE(c);
931 CAPTURE(d);
932 CHECK(check_abs(c, T, acc));
933 }
934
935 c = poly.solve_guess(matrix, 0.0, d, 350, 0);
936 {
937 CAPTURE(T);
938 CAPTURE(c);
939 CAPTURE(d);
940 CHECK(check_abs(c, T, acc));
941 }
942
943 // T = 0.0;
944 // solve.setGuess(75+273.15);
945 // T = solve.polyval(cHeat,c);
946 // printf("Should be : T = %3.3f \t K \n",273.15+50.0);
947 // printf("From object: T = %3.3f \t K \n",T);
948 //
949 // T = 0.0;
950 // solve.setLimits(273.15+10,273.15+100);
951 // T = solve.polyval(cHeat,c);
952 // printf("Should be : T = %3.3f \t K \n",273.15+50.0);
953 // printf("From object: T = %3.3f \t K \n",T);
954 }
955
956 SECTION("Integration and derivation tests") {
957
958 CoolProp::Polynomial2D poly;
959
960 const Eigen::MatrixXd& matrix(matrix2D);
961 Eigen::MatrixXd matrixInt = poly.integrateCoeffs(matrix, 1);
962 Eigen::MatrixXd matrixDer = poly.deriveCoeffs(matrix, 1);
963 Eigen::MatrixXd matrixInt2 = poly.integrateCoeffs(matrix, 1, 2);
964 Eigen::MatrixXd matrixDer2 = poly.deriveCoeffs(matrix, 1, 2);
965
966 CHECK_THROWS(poly.evaluate(matrix, 0.0));
967
968 double x = 0.3, y = 255.3, val1, val2, val3, val4;
969
970 //CHECK( std::abs( polyInt.derivative(x,y,0)-poly2D.evaluate(x,y) ) <= 1e-10 );
971
972 std::string tmpStr;
973
974 double acc = 0.001;
975
976 // Integer-indexed grid: preserves the original `T < Tmax` exit
977 // semantics without floating-point loop-counter accumulation
978 // (cert-flp30-c). N_T = ceil((Tmax-Tmin)/Tinc) matches the
979 // count the original `for (T = Tmin; T < Tmax; T += Tinc)` loop
980 // produces for the current Tmin/Tmax/Tinc (223.15/523.15/200 -> 2).
981 const auto N_T = static_cast<std::size_t>(std::ceil((Tmax - Tmin) / Tinc));
982
983 for (std::size_t i = 0; i < N_T; ++i) {
984 const double T = Tmin + i * Tinc;
985 val1 = poly.evaluate(matrix, x, T - deltaT);
986 val2 = poly.evaluate(matrix, x, T + deltaT);
987 val3 = (val2 - val1) / 2 / deltaT;
988 val4 = poly.evaluate(matrixDer, x, T);
989 CAPTURE(T);
990 CAPTURE(val3);
991 CAPTURE(val4);
992 tmpStr = CoolProp::mat_to_string(matrix);
993 CAPTURE(tmpStr);
994 tmpStr = CoolProp::mat_to_string(matrixDer);
995 CAPTURE(tmpStr);
996 CHECK(check_abs(val3, val4, acc));
997 }
998
999 for (std::size_t i = 0; i < N_T; ++i) {
1000 const double T = Tmin + i * Tinc;
1001 val1 = poly.evaluate(matrixDer, x, T - deltaT);
1002 val2 = poly.evaluate(matrixDer, x, T + deltaT);
1003 val3 = (val2 - val1) / 2 / deltaT;
1004 val4 = poly.evaluate(matrixDer2, x, T);
1005 CAPTURE(T);
1006 CAPTURE(val3);
1007 CAPTURE(val4);
1008 tmpStr = CoolProp::mat_to_string(matrixDer);
1009 CAPTURE(tmpStr);
1010 tmpStr = CoolProp::mat_to_string(matrixDer2);
1011 CAPTURE(tmpStr);
1012 CHECK(check_abs(val3, val4, acc));
1013 }
1014
1015 for (std::size_t i = 0; i < N_T; ++i) {
1016 const double T = Tmin + i * Tinc;
1017 val1 = poly.evaluate(matrixInt, x, T - deltaT);
1018 val2 = poly.evaluate(matrixInt, x, T + deltaT);
1019 val3 = (val2 - val1) / 2 / deltaT;
1020 val4 = poly.evaluate(matrix, x, T);
1021 CAPTURE(T);
1022 CAPTURE(val3);
1023 CAPTURE(val4);
1024 tmpStr = CoolProp::mat_to_string(matrixInt);
1025 CAPTURE(tmpStr);
1026 tmpStr = CoolProp::mat_to_string(matrix);
1027 CAPTURE(tmpStr);
1028 CHECK(check_abs(val3, val4, acc));
1029 }
1030
1031 for (std::size_t i = 0; i < N_T; ++i) {
1032 const double T = Tmin + i * Tinc;
1033 val1 = poly.evaluate(matrixInt2, x, T - deltaT);
1034 val2 = poly.evaluate(matrixInt2, x, T + deltaT);
1035 val3 = (val2 - val1) / 2 / deltaT;
1036 val4 = poly.evaluate(matrixInt, x, T);
1037 CAPTURE(T);
1038 CAPTURE(val3);
1039 CAPTURE(val4);
1040 tmpStr = CoolProp::mat_to_string(matrixInt2);
1041 CAPTURE(tmpStr);
1042 tmpStr = CoolProp::mat_to_string(matrixInt);
1043 CAPTURE(tmpStr);
1044 CHECK(check_abs(val3, val4, acc));
1045 }
1046
1047 for (std::size_t i = 0; i < N_T; ++i) {
1048 const double T = Tmin + i * Tinc;
1049 val1 = poly.evaluate(matrix, x, T);
1050 val2 = poly.derivative(matrixInt, x, T, 1);
1051 CAPTURE(T);
1052 CAPTURE(val1);
1053 CAPTURE(val2);
1054 CHECK(check_abs(val1, val2, acc));
1055 }
1056
1057 for (std::size_t i = 0; i < N_T; ++i) {
1058 const double T = Tmin + i * Tinc;
1059 val1 = poly.derivative(matrix, x, T, 1);
1060 val2 = poly.evaluate(matrixDer, x, T);
1061 CAPTURE(T);
1062 CAPTURE(val1);
1063 CAPTURE(val2);
1064 CHECK(check_abs(val1, val2, acc));
1065 }
1066 }
1067
1068 SECTION("Testing Polynomial2DFrac") {
1069
1070 Eigen::MatrixXd matrix = CoolProp::vec_to_eigen(cHeat);
1071 CoolProp::Polynomial2D poly;
1072 CoolProp::Polynomial2DFrac frac;
1073
1074 double acc = 0.0001;
1075
1076 double T = 273.15 + 50;
1077 double a, b;
1078 double c = poly.evaluate(matrix, T, 0.0);
1079 double d = frac.evaluate(matrix, T, 0.0, 0, 0);
1080
1081 {
1082 CAPTURE(T);
1083 CAPTURE(c);
1084 CAPTURE(d);
1085 tmpStr = CoolProp::mat_to_string(matrix);
1086 CAPTURE(tmpStr);
1087 CHECK(check_abs(c, d, acc));
1088 }
1089
1090 c = poly.evaluate(matrix, T, 0.0) / T / T;
1091 d = frac.evaluate(matrix, T, 0.0, -2, 0);
1092
1093 {
1094 CAPTURE(T);
1095 CAPTURE(c);
1096 CAPTURE(d);
1097 tmpStr = CoolProp::mat_to_string(matrix);
1098 CAPTURE(tmpStr);
1099 CHECK(check_abs(c, d, acc));
1100 }
1101
1102 matrix = CoolProp::vec_to_eigen(cHeat2D);
1103 double y = 0.1;
1104 c = poly.evaluate(matrix, T, y) / T / T;
1105 d = frac.evaluate(matrix, T, y, -2, 0);
1106
1107 {
1108 CAPTURE(T);
1109 CAPTURE(c);
1110 CAPTURE(d);
1111 tmpStr = CoolProp::mat_to_string(matrix);
1112 CAPTURE(tmpStr);
1113 CHECK(check_abs(c, d, acc));
1114 }
1115
1116 c = poly.evaluate(matrix, T, y) / y / y;
1117 d = frac.evaluate(matrix, T, y, 0, -2);
1118
1119 {
1120 CAPTURE(T);
1121 CAPTURE(c);
1122 CAPTURE(d);
1123 tmpStr = CoolProp::mat_to_string(matrix);
1124 CAPTURE(tmpStr);
1125 CHECK(check_abs(c, d, acc));
1126 }
1127
1128 c = poly.evaluate(matrix, T, y) / T / T / y / y;
1129 d = frac.evaluate(matrix, T, y, -2, -2);
1130
1131 {
1132 CAPTURE(T);
1133 CAPTURE(c);
1134 CAPTURE(d);
1135 tmpStr = CoolProp::mat_to_string(matrix);
1136 CAPTURE(tmpStr);
1137 CHECK(check_abs(c, d, acc));
1138 }
1139
1140 c = poly.evaluate(matrix, T, y) / T / T * y * y;
1141 d = frac.evaluate(matrix, T, y, -2, 2);
1142
1143 {
1144 CAPTURE(T);
1145 CAPTURE(c);
1146 CAPTURE(d);
1147 tmpStr = CoolProp::mat_to_string(matrix);
1148 CAPTURE(tmpStr);
1149 CHECK(check_abs(c, d, acc));
1150 }
1151
1152 matrix = CoolProp::vec_to_eigen(cHeat);
1153 T = 273.15 + 50;
1154 c = 145.59157247249246;
1155 d = frac.integral(matrix, T, 0.0, 0, -1, 0) - frac.integral(matrix, 273.15 + 25, 0.0, 0, -1, 0);
1156
1157 {
1158 CAPTURE(T);
1159 CAPTURE(c);
1160 CAPTURE(d);
1161 tmpStr = CoolProp::mat_to_string(matrix);
1162 CAPTURE(tmpStr);
1163 CHECK(check_abs(c, d, acc));
1164 }
1165
1166 T = 423.15;
1167 c = 3460.895272;
1168 d = frac.integral(matrix, T, 0.0, 0, -1, 0, 348.15, 0.0);
1169
1170 {
1171 CAPTURE(T);
1172 CAPTURE(c);
1173 CAPTURE(d);
1174 tmpStr = CoolProp::mat_to_string(matrix);
1175 CAPTURE(tmpStr);
1176 CHECK(check_abs(c, d, acc));
1177 }
1178
1179 deltaT = 0.01;
1180 // Integer-indexed grid (cert-flp30-c) — same N_T derivation as
1181 // the loops in the previous SECTION; T is the outer-scope test
1182 // variable, reuse it.
1183 const auto N_T_frac = static_cast<std::size_t>(std::ceil((Tmax - Tmin) / Tinc));
1184 for (std::size_t i = 0; i < N_T_frac; ++i) {
1185 T = Tmin + i * Tinc;
1186 a = poly.evaluate(matrix, T - deltaT, y);
1187 b = poly.evaluate(matrix, T + deltaT, y);
1188 c = (b - a) / 2.0 / deltaT;
1189 d = frac.derivative(matrix, T, y, 0, 0, 0);
1190 CAPTURE(a);
1191 CAPTURE(b);
1192 CAPTURE(T);
1193 CAPTURE(c);
1194 CAPTURE(d);
1195 tmpStr = CoolProp::mat_to_string(matrix);
1196 CAPTURE(tmpStr);
1197 CHECK(check_abs(c, d, acc));
1198 }
1199
1200 T = 273.15 + 150;
1201 c = -2.100108045;
1202 d = frac.derivative(matrix, T, 0.0, 0, 0, 0);
1203 {
1204 CAPTURE(T);
1205 CAPTURE(c);
1206 CAPTURE(d);
1207 tmpStr = CoolProp::mat_to_string(matrix);
1208 CAPTURE(tmpStr);
1209 CHECK(check_abs(c, d, acc));
1210 }
1211
1212 c = -0.006456574589;
1213 d = frac.derivative(matrix, T, 0.0, 0, -1, 0);
1214 {
1215 CAPTURE(T);
1216 CAPTURE(c);
1217 CAPTURE(d);
1218 tmpStr = CoolProp::mat_to_string(matrix);
1219 CAPTURE(tmpStr);
1220 CHECK(check_abs(c, d, acc));
1221 }
1222
1223 c = frac.evaluate(matrix, T, 0.0, 2, 0);
1224 d = frac.solve(matrix, 0.0, c, 0, 2, 0)[0];
1225 {
1226 CAPTURE(T);
1227 CAPTURE(c);
1228 CAPTURE(d);
1229 tmpStr = CoolProp::mat_to_string(matrix);
1230 CAPTURE(tmpStr);
1231 CHECK(check_abs(T, d, acc));
1232 }
1233
1234 c = frac.evaluate(matrix, T, 0.0, 0, 0);
1235 d = frac.solve(matrix, 0.0, c, 0, 0, 0)[0];
1236 {
1237 CAPTURE(T);
1238 CAPTURE(c);
1239 CAPTURE(d);
1240 tmpStr = CoolProp::mat_to_string(matrix);
1241 CAPTURE(tmpStr);
1242 CHECK(check_abs(T, d, acc));
1243 }
1244
1245 c = frac.evaluate(matrix, T, 0.0, -1, 0);
1246 CHECK_THROWS(d = frac.solve(matrix, 0.0, c, 0, -1, 0)[0]);
1247 // {
1248 // CAPTURE(T);
1249 // CAPTURE(c);
1250 // CAPTURE(d);
1251 // tmpStr = CoolProp::mat_to_string(matrix);
1252 // CAPTURE(tmpStr);
1253 // CHECK( check_abs(T,d,acc) );
1254 // }
1255
1256 d = frac.solve_limits(matrix, 0.0, c, T - 10, T + 10, 0, -1, 0);
1257 {
1258 CAPTURE(T);
1259 CAPTURE(c);
1260 CAPTURE(d);
1261 tmpStr = CoolProp::mat_to_string(matrix);
1262 CAPTURE(tmpStr);
1263 CHECK(check_abs(T, d, acc));
1264 }
1265
1266 d = frac.solve_guess(matrix, 0.0, c, T - 10, 0, -1, 0);
1267 {
1268 CAPTURE(T);
1269 CAPTURE(c);
1270 CAPTURE(d);
1271 tmpStr = CoolProp::mat_to_string(matrix);
1272 CAPTURE(tmpStr);
1273 CHECK(check_abs(T, d, acc));
1274 }
1275
1276 c = -0.00004224550082;
1277 d = frac.derivative(matrix, T, 0.0, 0, -2, 0);
1278 {
1279 CAPTURE(T);
1280 CAPTURE(c);
1281 CAPTURE(d);
1282 tmpStr = CoolProp::mat_to_string(matrix);
1283 CAPTURE(tmpStr);
1284 CHECK(check_abs(c, d, acc));
1285 }
1286
1287 c = frac.evaluate(matrix, T, 0.0, 0, 0, 0.0, 0.0);
1288 d = frac.solve(matrix, 0.0, c, 0, 0, 0, 0.0, 0.0)[0];
1289 {
1290 CAPTURE(T);
1291 CAPTURE(c);
1292 CAPTURE(d);
1293 tmpStr = CoolProp::mat_to_string(matrix);
1294 CAPTURE(tmpStr);
1295 tmpStr = CoolProp::mat_to_string(Eigen::MatrixXd(frac.solve(matrix, 0.0, c, 0, 0, 0, 250, 0.0)));
1296 CAPTURE(tmpStr);
1297 CHECK(check_abs(T, d, acc));
1298 }
1299 }
1300}
1301
1302#endif /* ENABLE_CATCH */
Generated by doxygen 1.9.4