1 #ifndef COOLPROP_NUMERICS_H
2 #define COOLPROP_NUMERICS_H
15 #if defined(HUGE_VAL) && !defined(_HUGE)
16 # define _HUGE HUGE_VAL
19 # if defined(HUGE) && !defined(_HUGE)
26 return (x <= DBL_MAX && x >= -DBL_MAX);
30 # define M_PI 3.14159265358979323846
34 # define COOLPROP_OK 1
51 template <
typename X,
typename Y>
59 Spline(
const std::vector<X>& x,
const std::vector<Y>& y) {
60 if (x.size() != y.size()) {
61 std::cerr <<
"X and Y must be the same size " << std::endl;
66 std::cerr <<
"Must have at least three points for interpolation" << std::endl;
70 typedef typename std::vector<X>::difference_type size_type;
72 size_type n = y.size() - 1;
74 std::vector<Y> b(n), d(n), a(n), c(n + 1), l(n + 1), u(n + 1), z(n + 1);
75 std::vector<X> h(n + 1);
82 for (size_type i = 1; i < n; i++) {
83 h[i] = x[i + 1] - x[i];
84 l[i] = Y(2 * (x[i + 1] - x[i - 1])) - Y(h[i - 1]) * u[i - 1];
85 u[i] = Y(h[i]) / l[i];
86 a[i] = (Y(3) / Y(h[i])) * (y[i + 1] - y[i]) - (Y(3) / Y(h[i - 1])) * (y[i] - y[i - 1]);
87 z[i] = (a[i] - Y(h[i - 1]) * z[i - 1]) / l[i];
93 for (size_type j = n - 1; j >= 0; j--) {
94 c[j] = z[j] - u[j] * c[j + 1];
95 b[j] = (y[j + 1] - y[j]) / Y(h[j]) - (Y(h[j]) * (c[j + 1] + Y(2) * c[j])) / Y(3);
96 d[j] = (c[j + 1] - c[j]) / Y(3 * h[j]);
99 for (size_type i = 0; i < n; i++) {
112 typename std::vector<element_type>::const_iterator it;
123 if (
mElements.size() == 0)
return std::vector<Y>(xx.size());
125 typename std::vector<X>::const_iterator it;
126 typename std::vector<element_type>::const_iterator it2;
129 for (it = xx.begin(); it != xx.end(); it++) {
135 ys.push_back(it2->eval(*it));
146 Element(
X _x, Y _a, Y _b, Y _c, Y _d) :
x(_x),
a(_a),
b(_b),
c(_c),
d(_d) {}
150 return a +
b * xix +
c * (xix * xix) +
d * (xix * xix * xix);
169 template <
typename T>
172 for (std::size_t i = 0; i < z1.size(); ++i) {
173 T diff = std::abs(z1[i] - z2[i]);
174 if (std::abs(diff) > maxvecdiff) {
182 template <
typename T>
183 std::vector<T>
linspace(T xmin, T xmax, std::size_t n) {
184 std::vector<T> x(n, 0.0);
186 for (std::size_t i = 0; i < n; ++i) {
187 x[i] = (xmax - xmin) / (n - 1) * i + xmin;
192 template <
typename T>
194 std::vector<T> x(n, 0.0);
195 T logxmin = log10(xmin), logxmax = log10(xmax);
197 for (std::size_t i = 0; i < n; ++i) {
198 x[i] = exp((logxmax - logxmin) / (n - 1) * i + logxmin);
203 template <
typename T>
204 std::vector<T>
logspace(T xmin, T xmax, std::size_t n) {
205 std::vector<T> x(n, 0.0);
206 T logxmin = log(xmin), logxmax = log(xmax);
208 for (std::size_t i = 0; i < n; ++i) {
209 x[i] = exp((logxmax - logxmin) / (n - 1) * i + logxmin);
222 template <
typename T>
225 std::size_t N = vec.size(), L = 0, R = N - 1, M = (L + R) / 2;
235 if (L == vec.size() - 1) {
244 std::size_t MR = M, ML = M;
247 if (MR == vec.size() - 1) {
259 T rML = vec[ML] - val;
260 T rMR = vec[MR] - val;
262 if (rR * rML > 0 && rL * rML < 0) {
266 }
else if (rR * rMR < 0 && rL * rMR > 0) {
272 format(
"Unable to bisect segmented vector; neither chunk contains the solution val:%g left:(%g,%g) right:(%g,%g)", val, vec[L],
273 vec[ML], vec[MR], vec[R]));
278 if (rR * rM > 0 && rL * rM < 0) {
302 template <
typename T>
305 std::size_t N =
mat[j].size(), L = 0, R = N - 1, M = (L + R) / 2;
313 rR =
mat[R][j] - val;
316 if (L ==
mat.size() - 1) {
321 rL =
mat[L][j] - val;
324 std::size_t MR = M, ML = M;
327 if (MR ==
mat.size() - 1) {
339 T rML =
mat[ML][j] - val;
340 T rMR =
mat[MR][j] - val;
342 if (rR * rMR > 0 && rL * rML < 0) {
345 rR =
mat[ML][j] - val;
346 }
else if (rR * rMR < 0 && rL * rML > 0) {
349 rL =
mat[MR][j] - val;
352 format(
"Unable to bisect segmented vector slice; neither chunk contains the solution %g lef:(%g,%g) right:(%g,%g)", val,
mat[L][j],
357 rM =
mat[M][j] - val;
358 if (rR * rM > 0 && rL * rM < 0) {
361 rR =
mat[R][j] - val;
365 rL =
mat[L][j] - val;
379 inline std::set<std::set<std::size_t>>
powerset(std::set<std::size_t>
const& set) {
380 std::set<std::set<std::size_t>> result;
381 std::vector<std::set<std::size_t>::const_iterator> elements;
383 std::set<std::size_t> tmp;
384 std::transform(elements.begin(), elements.end(), std::inserter(tmp, tmp.end()),
powerset_dereference);
386 if (!elements.empty() && ++elements.back() == set.end()) {
389 std::set<std::size_t>::const_iterator iter;
390 if (elements.empty()) {
393 iter = elements.back();
396 for (; iter != set.end(); ++iter) {
397 elements.push_back(iter);
400 }
while (!elements.empty());
407 double max = std::abs(A);
408 double min = std::abs(B);
414 return ((1.0 - min / max * 1e0) < D);
417 format(
"Too small numbers: %f cannot be tested with an accepted error of %f for a machine precision of %f. ", max, D,
DBL_EPSILON));
426 T sumx = std::accumulate(x.begin(), x.end(),
static_cast<T
>(0));
428 for (std::size_t i = 0; i < x.size(); ++i) {
440 std::vector<std::vector<double>>
A;
441 std::vector<double>
B;
448 void add_4value_constraints(
double x1,
double x2,
double x3,
double x4,
double y1,
double y2,
double y3,
double y4);
456 if (n == 0)
return 1;
461 template <
class T1,
class T2>
470 void MatInv_2(
double A[2][2],
double B[2][2]);
473 double interp1d(
const std::vector<double>* x,
const std::vector<double>* y,
double x0);
474 double powInt(
double x,
int y);
488 #define POW5(x) ((x) * (x) * (x) * (x) * (x))
489 #define POW6(x) ((x) * (x) * (x) * (x) * (x) * (x))
490 #define POW7(x) ((x) * (x) * (x) * (x) * (x) * (x) * (x))
494 return (y1 - y0) / (x1 - x0) * (x - x0) + y0;
496 template <
class T1,
class T2>
497 T2
LinearInterp(
const std::vector<T1>& x,
const std::vector<T1>& y, std::size_t i0, std::size_t
i1, T2 val) {
508 L0 = ((x - x1) * (x - x2)) / ((x0 - x1) * (x0 - x2));
509 L1 = ((x - x0) * (x - x2)) / ((x1 - x0) * (x1 - x2));
510 L2 = ((x - x0) * (x - x1)) / ((x2 - x0) * (x2 - x1));
511 return L0 * f0 + L1 * f1 + L2 * f2;
513 template <
class T1,
class T2>
514 T2
QuadInterp(
const std::vector<T1>& x,
const std::vector<T1>& y, std::size_t i0, std::size_t
i1, std::size_t i2, T2 val) {
515 return QuadInterp(x[i0], x[
i1], x[i2], y[i0], y[
i1], y[i2],
static_cast<T1
>(val));
519 T
CubicInterp(T x0, T x1, T x2, T x3, T f0, T f1, T f2, T f3, T x) {
525 L0 = ((x - x1) * (x - x2) * (x - x3)) / ((x0 - x1) * (x0 - x2) * (x0 - x3));
526 L1 = ((x - x0) * (x - x2) * (x - x3)) / ((x1 - x0) * (x1 - x2) * (x1 - x3));
527 L2 = ((x - x0) * (x - x1) * (x - x3)) / ((x2 - x0) * (x2 - x1) * (x2 - x3));
528 L3 = ((x - x0) * (x - x1) * (x - x2)) / ((x3 - x0) * (x3 - x1) * (x3 - x2));
529 return L0 * f0 + L1 * f1 + L2 * f2 + L3 * f3;
536 T L0 = ((x - x1) * (x - x2) * (x - x3)) / ((x0 - x1) * (x0 - x2) * (x0 - x3));
537 T dL0_dx = (1 / (x - x1) + 1 / (x - x2) + 1 / (x - x3)) * L0;
538 T L1 = ((x - x0) * (x - x2) * (x - x3)) / ((x1 - x0) * (x1 - x2) * (x1 - x3));
539 T dL1_dx = (1 / (x - x0) + 1 / (x - x2) + 1 / (x - x3)) * L1;
540 T L2 = ((x - x0) * (x - x1) * (x - x3)) / ((x2 - x0) * (x2 - x1) * (x2 - x3));
541 T dL2_dx = (1 / (x - x0) + 1 / (x - x1) + 1 / (x - x3)) * L2;
542 T L3 = ((x - x0) * (x - x1) * (x - x2)) / ((x3 - x0) * (x3 - x1) * (x3 - x2));
543 T dL3_dx = (1 / (x - x0) + 1 / (x - x1) + 1 / (x - x2)) * L3;
544 return dL0_dx * f0 + dL1_dx * f1 + dL2_dx * f2 + dL3_dx * f3;
546 template <
class T1,
class T2>
547 T2
CubicInterp(
const std::vector<T1>& x,
const std::vector<T1>& y, std::size_t i0, std::size_t
i1, std::size_t i2, std::size_t i3, T2 val) {
548 return CubicInterp(x[i0], x[
i1], x[i2], x[i3], y[i0], y[
i1], y[i2], y[i3],
static_cast<T1
>(val));
553 return (x >= std::min(x1, x2) && x <= std::max(x1, x2));
569 void solve_cubic(
double a,
double b,
double c,
double d,
int& N,
double& x0,
double& x1,
double& x2);
571 void solve_quartic(
double a,
double b,
double c,
double d,
double e,
int& N,
double& x0,
double& x1,
double& x2,
double& x3);
574 inline T
min3(T x1, T x2, T x3) {
575 return std::min(std::min(x1, x2), x3);
578 inline T
max3(T x1, T x2, T x3) {
579 return std::max(std::max(x1, x2), x3);
582 inline T
min4(T x1, T x2, T x3, T x4) {
583 return std::min(std::min(std::min(x1, x2), x3), x4);
586 inline T
max4(T x1, T x2, T x3, T x4) {
587 return std::max(std::max(std::max(x1, x2), x3), x4);
591 return std::abs(a - b) <= 1 *
DBL_EPSILON * std::max(std::abs(a), std::abs(b));
597 std::size_t N = x.size();
598 for (std::size_t i = 0; i < N; ++i) {
599 T axi = std::abs(x[i]);
610 std::size_t N = x.size();
611 for (std::size_t i = 0; i < N; ++i) {
612 T axi = std::abs(x[i]);
622 return static_cast<int>(1);
624 return static_cast<int>(0);
636 template <
typename T>
669 #if defined(_MSC_VER)
673 double acosh(
double x) {
674 return log(x + sqrt(x * x - 1.0));
676 double asinh(
double value) {
678 return log(value + sqrt(value * value + 1));
680 return -log(-value + sqrt(value * value + 1));
686 #if defined(__ISPOWERPC__)
688 double acosh(
double x) {
689 return log(x + sqrt(x * x - 1.0));
691 double asinh(
double value) {
693 return log(value + sqrt(value * value + 1));
695 return -log(-value + sqrt(value * value + 1));
700 #if defined(__ISPOWERPC__)