CoolProp  6.6.1dev
An open-source fluid property and humid air property database
CPnumerics.h
Go to the documentation of this file.
1 #ifndef COOLPROP_NUMERICS_H
2 #define COOLPROP_NUMERICS_H
3 
4 #include <vector>
5 #include <set>
6 #include <cfloat>
7 #include <stdlib.h> // For abs
8 #include <algorithm> // For max
9 #include <numeric>
10 #include <cmath>
11 #include "PlatformDetermination.h"
12 #include "CPstrings.h"
13 #include "Exceptions.h"
14 
15 #if defined(HUGE_VAL) && !defined(_HUGE)
16 # define _HUGE HUGE_VAL
17 #else
18 // GCC Version of huge value macro
19 # if defined(HUGE) && !defined(_HUGE)
20 # define _HUGE HUGE
21 # endif
22 #endif
23 
24 inline bool ValidNumber(double x) {
25  // Idea from http://www.johndcook.com/IEEE_exceptions_in_cpp.html
26  return (x <= DBL_MAX && x >= -DBL_MAX);
27 };
28 
29 #ifndef M_PI
30 # define M_PI 3.14159265358979323846
31 #endif
32 
33 #ifndef COOLPROP_OK
34 # define COOLPROP_OK 1
35 #endif
36 
37 // Undefine these terrible macros defined in windows header
38 #undef min
39 #undef max
40 
41 /* "THE BEER-WARE LICENSE" (Revision 42): Devin Lane wrote this file. As long as you retain
42  * this notice you can do whatever you want with this stuff. If we meet some day, and you
43  * think this stuff is worth it, you can buy me a beer in return.
44  *
45  * From http://shiftedbits.org/2011/01/30/cubic-spline-interpolation/
46  *
47  * IHB(05/01/2016): Removed overload and renamed the interpolate function (cython cannot disambiguate the functions)
48  *
49  * Templated on type of X, Y. X and Y must have operator +, -, *, /. Y must have defined
50  * a constructor that takes a scalar. */
51 template <typename X, typename Y>
52 class Spline
53 {
54  public:
56  Spline() {}
57 
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;
62  return;
63  }
64 
65  if (x.size() < 3) {
66  std::cerr << "Must have at least three points for interpolation" << std::endl;
67  return;
68  }
69 
70  typedef typename std::vector<X>::difference_type size_type;
71 
72  size_type n = y.size() - 1;
73 
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);
76 
77  l[0] = Y(1);
78  u[0] = Y(0);
79  z[0] = Y(0);
80  h[0] = x[1] - x[0];
81 
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];
88  }
89 
90  l[n] = Y(1);
91  z[n] = c[n] = Y(0);
92 
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]);
97  }
98 
99  for (size_type i = 0; i < n; i++) {
100  mElements.push_back(Element(x[i], y[i], b[i], c[i], d[i]));
101  }
102  }
103  virtual ~Spline() {}
104 
105  Y operator[](const X& x) const {
106  return interpolate(x);
107  }
108 
109  Y interpolate(const X& x) const {
110  if (mElements.size() == 0) return Y();
111 
112  typename std::vector<element_type>::const_iterator it;
113  it = std::lower_bound(mElements.begin(), mElements.end(), element_type(x));
114  if (it != mElements.begin()) {
115  it--;
116  }
117 
118  return it->eval(x);
119  }
120 
121  /* Evaluate at multiple locations, assuming xx is sorted ascending */
122  std::vector<Y> interpolate_vec(const std::vector<X>& xx) const {
123  if (mElements.size() == 0) return std::vector<Y>(xx.size());
124 
125  typename std::vector<X>::const_iterator it;
126  typename std::vector<element_type>::const_iterator it2;
127  it2 = mElements.begin();
128  std::vector<Y> ys;
129  for (it = xx.begin(); it != xx.end(); it++) {
130  it2 = std::lower_bound(it2, mElements.end(), element_type(*it));
131  if (it2 != mElements.begin()) {
132  it2--;
133  }
134 
135  ys.push_back(it2->eval(*it));
136  }
137 
138  return ys;
139  }
140 
141  protected:
142  class Element
143  {
144  public:
145  Element(X _x) : x(_x) {}
146  Element(X _x, Y _a, Y _b, Y _c, Y _d) : x(_x), a(_a), b(_b), c(_c), d(_d) {}
147 
148  Y eval(const X& xx) const {
149  X xix(xx - x);
150  return a + b * xix + c * (xix * xix) + d * (xix * xix * xix);
151  }
152 
153  bool operator<(const Element& e) const {
154  return x < e.x;
155  }
156  bool operator<(const X& xx) const {
157  return x < xx;
158  }
159 
160  X x;
161  Y a, b, c, d;
162  };
163 
165  std::vector<element_type> mElements;
166 };
167 
169 template <typename T>
170 T maxvectordiff(const std::vector<T>& z1, const std::vector<T>& z2) {
171  T maxvecdiff = 0;
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) {
175  maxvecdiff = diff;
176  }
177  }
178  return maxvecdiff;
179 }
180 
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);
185 
186  for (std::size_t i = 0; i < n; ++i) {
187  x[i] = (xmax - xmin) / (n - 1) * i + xmin;
188  }
189  return x;
190 }
192 template <typename T>
193 std::vector<T> log10space(T xmin, T xmax, std::size_t n) {
194  std::vector<T> x(n, 0.0);
195  T logxmin = log10(xmin), logxmax = log10(xmax);
196 
197  for (std::size_t i = 0; i < n; ++i) {
198  x[i] = exp((logxmax - logxmin) / (n - 1) * i + logxmin);
199  }
200  return x;
201 }
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);
207 
208  for (std::size_t i = 0; i < n; ++i) {
209  x[i] = exp((logxmax - logxmin) / (n - 1) * i + logxmin);
210  }
211  return x;
212 }
213 
222 template <typename T>
223 void bisect_vector(const std::vector<T>& vec, T val, std::size_t& i) {
224  T rL, rM, rR;
225  std::size_t N = vec.size(), L = 0, R = N - 1, M = (L + R) / 2;
226  // Move the right limits in until they are good
227  while (!ValidNumber(vec[R])) {
228  if (R == 1) {
229  throw CoolProp::ValueError("All the values in bisection vector are invalid");
230  }
231  R--;
232  }
233  // Move the left limits in until they are good
234  while (!ValidNumber(vec[L])) {
235  if (L == vec.size() - 1) {
236  throw CoolProp::ValueError("All the values in bisection vector are invalid");
237  }
238  L++;
239  }
240  rL = vec[L] - val;
241  rR = vec[R] - val;
242  while (R - L > 1) {
243  if (!ValidNumber(vec[M])) {
244  std::size_t MR = M, ML = M;
245  // Move middle-right to the right until it is ok
246  while (!ValidNumber(vec[MR])) {
247  if (MR == vec.size() - 1) {
248  throw CoolProp::ValueError("All the values in bisection vector are invalid");
249  }
250  MR++;
251  }
252  // Move middle-left to the left until it is ok
253  while (!ValidNumber(vec[ML])) {
254  if (ML == 1) {
255  throw CoolProp::ValueError("All the values in bisection vector are invalid");
256  }
257  ML--;
258  }
259  T rML = vec[ML] - val;
260  T rMR = vec[MR] - val;
261  // Figure out which chunk is the good part
262  if (rR * rML > 0 && rL * rML < 0) {
263  // solution is between L and ML
264  R = ML;
265  rR = vec[ML] - val;
266  } else if (rR * rMR < 0 && rL * rMR > 0) {
267  // solution is between R and MR
268  L = MR;
269  rL = vec[MR] - val;
270  } else {
271  throw CoolProp::ValueError(
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]));
274  }
275  M = (L + R) / 2;
276  } else {
277  rM = vec[M] - val;
278  if (rR * rM > 0 && rL * rM < 0) {
279  // solution is between L and M
280  R = M;
281  rR = vec[R] - val;
282  } else {
283  // solution is between R and M
284  L = M;
285  rL = vec[L] - val;
286  }
287  M = (L + R) / 2;
288  }
289  }
290  i = L;
291 }
292 
302 template <typename T>
303 void bisect_segmented_vector_slice(const std::vector<std::vector<T>>& mat, std::size_t j, T val, std::size_t& i) {
304  T rL, rM, rR;
305  std::size_t N = mat[j].size(), L = 0, R = N - 1, M = (L + R) / 2;
306  // Move the right limits in until they are good
307  while (!ValidNumber(mat[R][j])) {
308  if (R == 1) {
309  throw CoolProp::ValueError("All the values in bisection vector are invalid");
310  }
311  R--;
312  }
313  rR = mat[R][j] - val;
314  // Move the left limits in until they are good
315  while (!ValidNumber(mat[L][j])) {
316  if (L == mat.size() - 1) {
317  throw CoolProp::ValueError("All the values in bisection vector are invalid");
318  }
319  L++;
320  }
321  rL = mat[L][j] - val;
322  while (R - L > 1) {
323  if (!ValidNumber(mat[M][j])) {
324  std::size_t MR = M, ML = M;
325  // Move middle-right to the right until it is ok
326  while (!ValidNumber(mat[MR][j])) {
327  if (MR == mat.size() - 1) {
328  throw CoolProp::ValueError("All the values in bisection vector are invalid");
329  }
330  MR++;
331  }
332  // Move middle-left to the left until it is ok
333  while (!ValidNumber(mat[ML][j])) {
334  if (ML == 1) {
335  throw CoolProp::ValueError("All the values in bisection vector are invalid");
336  }
337  ML--;
338  }
339  T rML = mat[ML][j] - val;
340  T rMR = mat[MR][j] - val;
341  // Figure out which chunk is the good part
342  if (rR * rMR > 0 && rL * rML < 0) {
343  // solution is between L and ML
344  R = ML;
345  rR = mat[ML][j] - val;
346  } else if (rR * rMR < 0 && rL * rML > 0) {
347  // solution is between R and MR
348  L = MR;
349  rL = mat[MR][j] - val;
350  } else {
351  throw CoolProp::ValueError(
352  format("Unable to bisect segmented vector slice; neither chunk contains the solution %g lef:(%g,%g) right:(%g,%g)", val, mat[L][j],
353  mat[ML][j], mat[MR][j], mat[R][j]));
354  }
355  M = (L + R) / 2;
356  } else {
357  rM = mat[M][j] - val;
358  if (rR * rM > 0 && rL * rM < 0) {
359  // solution is between L and M
360  R = M;
361  rR = mat[R][j] - val;
362  } else {
363  // solution is between R and M
364  L = M;
365  rL = mat[L][j] - val;
366  }
367  M = (L + R) / 2;
368  }
369  }
370  i = L;
371 }
372 
373 // From http://rosettacode.org/wiki/Power_set#C.2B.2B
374 inline std::size_t powerset_dereference(std::set<std::size_t>::const_iterator v) {
375  return *v;
376 };
377 
378 // From http://rosettacode.org/wiki/Power_set#C.2B.2B
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;
382  do {
383  std::set<std::size_t> tmp;
384  std::transform(elements.begin(), elements.end(), std::inserter(tmp, tmp.end()), powerset_dereference);
385  result.insert(tmp);
386  if (!elements.empty() && ++elements.back() == set.end()) {
387  elements.pop_back();
388  } else {
389  std::set<std::size_t>::const_iterator iter;
390  if (elements.empty()) {
391  iter = set.begin();
392  } else {
393  iter = elements.back();
394  ++iter;
395  }
396  for (; iter != set.end(); ++iter) {
397  elements.push_back(iter);
398  }
399  }
400  } while (!elements.empty());
401 
402  return result;
403 }
404 
406 bool inline check_abs(double A, double B, double D) {
407  double max = std::abs(A);
408  double min = std::abs(B);
409  if (min > max) {
410  max = min;
411  min = std::abs(A);
412  }
413  if (max > DBL_EPSILON * 1e3)
414  return ((1.0 - min / max * 1e0) < D);
415  else
416  throw CoolProp::ValueError(
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));
418 };
419 bool inline check_abs(double A, double B) {
420  return check_abs(A, B, 1e5 * DBL_EPSILON);
421 };
422 
423 template <class T>
424 void normalize_vector(std::vector<T>& x) {
425  // Sum up all the elements in the vector
426  T sumx = std::accumulate(x.begin(), x.end(), static_cast<T>(0));
427  // Normalize the components by dividing each by the sum
428  for (std::size_t i = 0; i < x.size(); ++i) {
429  x[i] /= sumx;
430  }
431 };
432 
437 {
438  protected:
440  std::vector<std::vector<double>> A;
441  std::vector<double> B;
442 
443  public:
444  double a, b, c, d;
445  SplineClass() : Nconstraints(0), A(4, std::vector<double>(4, 0)), B(4, 0), a(_HUGE), b(_HUGE), c(_HUGE), d(_HUGE) {}
446  bool build(void);
447  bool add_value_constraint(double x, double y);
448  void add_4value_constraints(double x1, double x2, double x3, double x4, double y1, double y2, double y3, double y4);
449  bool add_derivative_constraint(double x, double dydx);
450  double evaluate(double x);
451 };
452 
454 template <class T>
455 T factorial(T n) {
456  if (n == 0) return 1;
457  return n * factorial(n - 1);
458 }
461 template <class T1, class T2>
462 T1 nth_derivative_of_x_to_m(T1 x, T2 n, T2 m) {
463  if (n > m) {
464  return 0;
465  } else {
466  return factorial(m) / factorial(m - n) * pow(x, m - n);
467  }
468 }
469 
470 void MatInv_2(double A[2][2], double B[2][2]);
471 
472 double root_sum_square(const std::vector<double>& x);
473 double interp1d(const std::vector<double>* x, const std::vector<double>* y, double x0);
474 double powInt(double x, int y);
475 
476 template <class T>
477 T POW2(T x) {
478  return x * x;
479 }
480 template <class T>
481 T POW3(T x) {
482  return POW2(x) * x;
483 }
484 template <class T>
485 T POW4(T x) {
486  return POW2(x) * POW2(x);
487 }
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))
491 
492 template <class T>
493 T LinearInterp(T x0, T x1, T y0, T y1, T x) {
494  return (y1 - y0) / (x1 - x0) * (x - x0) + y0;
495 };
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) {
498  return LinearInterp(x[i0], x[i1], y[i0], y[i1], static_cast<T1>(val));
499 };
500 
501 template <class T>
502 T QuadInterp(T x0, T x1, T x2, T f0, T f1, T f2, T x) {
503  /* Quadratic interpolation.
504  Based on method from Kreyszig,
505  Advanced Engineering Mathematics, 9th Edition
506  */
507  T L0, L1, L2;
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;
512 };
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));
516 };
517 
518 template <class T>
519 T CubicInterp(T x0, T x1, T x2, T x3, T f0, T f1, T f2, T f3, T x) {
520  /*
521  Lagrange cubic interpolation as from
522  http://nd.edu/~jjwteach/441/PdfNotes/lecture6.pdf
523  */
524  T L0, L1, L2, L3;
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;
530 };
533 template <class T>
534 T CubicInterpFirstDeriv(T x0, T x1, T x2, T x3, T f0, T f1, T f2, T f3, T x) {
535  // Based on http://math.stackexchange.com/a/809946/66405
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;
545 };
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));
549 };
550 
551 template <class T>
552 T is_in_closed_range(T x1, T x2, T x) {
553  return (x >= std::min(x1, x2) && x <= std::max(x1, x2));
554 };
555 
569 void solve_cubic(double a, double b, double c, double d, int& N, double& x0, double& x1, double& x2);
570 
571 void solve_quartic(double a, double b, double c, double d, double e, int& N, double& x0, double& x1, double& x2, double& x3);
572 
573 template <class T>
574 inline T min3(T x1, T x2, T x3) {
575  return std::min(std::min(x1, x2), x3);
576 };
577 template <class T>
578 inline T max3(T x1, T x2, T x3) {
579  return std::max(std::max(x1, x2), x3);
580 };
581 template <class T>
582 inline T min4(T x1, T x2, T x3, T x4) {
583  return std::min(std::min(std::min(x1, x2), x3), x4);
584 };
585 template <class T>
586 inline T max4(T x1, T x2, T x3, T x4) {
587  return std::max(std::max(std::max(x1, x2), x3), x4);
588 };
589 
590 inline bool double_equal(double a, double b) {
591  return std::abs(a - b) <= 1 * DBL_EPSILON * std::max(std::abs(a), std::abs(b));
592 };
593 
594 template <class T>
595 T max_abs_value(const std::vector<T>& x) {
596  T max = 0;
597  std::size_t N = x.size();
598  for (std::size_t i = 0; i < N; ++i) {
599  T axi = std::abs(x[i]);
600  if (axi > max) {
601  max = axi;
602  }
603  }
604  return max;
605 }
606 
607 template <class T>
608 T min_abs_value(const std::vector<T>& x) {
609  T min = 1e40;
610  std::size_t N = x.size();
611  for (std::size_t i = 0; i < N; ++i) {
612  T axi = std::abs(x[i]);
613  if (axi < min) {
614  min = axi;
615  }
616  }
617  return min;
618 }
619 
620 inline int Kronecker_delta(std::size_t i, std::size_t j) {
621  if (i == j) {
622  return static_cast<int>(1);
623  } else {
624  return static_cast<int>(0);
625  }
626 };
627 inline int Kronecker_delta(int i, int j) {
628  if (i == j) {
629  return 1;
630  } else {
631  return 0;
632  }
633 };
634 
636 template <typename T>
637 void sort3(T& a, T& b, T& c) {
638  if (a > b) {
639  std::swap(a, b);
640  }
641  if (a > c) {
642  std::swap(a, c);
643  }
644  if (b > c) {
645  std::swap(b, c);
646  }
647 }
648 
659 template <class T>
660 T angle_difference(T angle1, T angle2) {
661  return fmod(angle1 - angle2 + M_PI, 2 * M_PI) - M_PI;
662 }
663 
665 inline double get_HUGE() {
666  return _HUGE;
667 }
668 
669 #if defined(_MSC_VER)
670 // Microsoft version of math.h doesn't include acosh or asinh, so we just define them here.
671 // It was included from Visual Studio 2013
672 # if _MSC_VER < 1800
673 double acosh(double x) {
674  return log(x + sqrt(x * x - 1.0));
675 }
676 double asinh(double value) {
677  if (value > 0) {
678  return log(value + sqrt(value * value + 1));
679  } else {
680  return -log(-value + sqrt(value * value + 1));
681  }
682 }
683 # endif
684 #endif
685 
686 #if defined(__ISPOWERPC__)
687 // PPC version of math.h doesn't include acosh or asinh, so we just define them here
688 double acosh(double x) {
689  return log(x + sqrt(x * x - 1.0));
690 }
691 double asinh(double value) {
692  if (value > 0) {
693  return log(value + sqrt(value * value + 1));
694  } else {
695  return -log(-value + sqrt(value * value + 1));
696  }
697 }
698 #endif
699 
700 #if defined(__ISPOWERPC__)
701 # undef EOS
702 #endif
703 
704 #endif