CoolProp  6.6.0
An open-source fluid property and humid air property database
CPnumerics.cpp
Go to the documentation of this file.
1 #include "CPnumerics.h"
2 #include "MatrixMath.h"
3 #include <unsupported/Eigen/Polynomials>
4 
5 double root_sum_square(const std::vector<double>& x) {
6  double sum = 0;
7  for (unsigned int i = 0; i < x.size(); i++) {
8  sum += pow(x[i], 2);
9  }
10  return sqrt(sum);
11 }
12 double interp1d(const std::vector<double>* x, const std::vector<double>* y, double x0) {
13  std::size_t i, L, R, M;
14  L = 0;
15  R = (*x).size() - 1;
16  M = (L + R) / 2;
17  // Use interval halving to find the indices which bracket the density of interest
18  while (R - L > 1) {
19  if (x0 >= (*x)[M]) {
20  L = M;
21  M = (L + R) / 2;
22  continue;
23  }
24  if (x0 < (*x)[M]) {
25  R = M;
26  M = (L + R) / 2;
27  continue;
28  }
29  }
30  i = L;
31  if (i < (*x).size() - 2) {
32  // Go "forwards" with the interpolation range
33  return QuadInterp((*x)[i], (*x)[i + 1], (*x)[i + 2], (*y)[i], (*y)[i + 1], (*y)[i + 2], x0);
34  } else {
35  // Go "backwards" with the interpolation range
36  return QuadInterp((*x)[i], (*x)[i - 1], (*x)[i - 2], (*y)[i], (*y)[i - 1], (*y)[i - 2], x0);
37  }
38 }
39 double powInt(double x, int y) {
40  // Raise a double to an integer power
41  // Overload not provided in math.h
42  int i;
43  double product = 1.0;
44  double x_in;
45  int y_in;
46 
47  if (y == 0) {
48  return 1.0;
49  }
50 
51  if (y < 0) {
52  x_in = 1 / x;
53  y_in = -y;
54  } else {
55  x_in = x;
56  y_in = y;
57  }
58 
59  if (y_in == 1) {
60  return x_in;
61  }
62 
63  product = x_in;
64  for (i = 1; i < y_in; i++) {
65  product = product * x_in;
66  }
67  return product;
68 }
69 
70 void MatInv_2(double A[2][2], double B[2][2]) {
71  double Det;
72  //Using Cramer's Rule to solve
73 
74  Det = A[0][0] * A[1][1] - A[1][0] * A[0][1];
75  B[0][0] = 1.0 / Det * A[1][1];
76  B[1][1] = 1.0 / Det * A[0][0];
77  B[1][0] = -1.0 / Det * A[1][0];
78  B[0][1] = -1.0 / Det * A[0][1];
79 }
80 
81 void solve_cubic(double a, double b, double c, double d, int& N, double& x0, double& x1, double& x2) {
82  // 0 = ax^3 + b*x^2 + c*x + d
83 
84  // First check if the "cubic" is actually a second order or first order curve
85  if (std::abs(a) < 10 * DBL_EPSILON) {
86  if (std::abs(b) < 10 * DBL_EPSILON) {
87  // Linear solution if a = 0 and b = 0
88  x0 = -d / c;
89  N = 1;
90  return;
91  } else {
92  // Quadratic solution(s) if a = 0 and b != 0
93  x0 = (-c + sqrt(c * c - 4 * b * d)) / (2 * b);
94  x1 = (-c - sqrt(c * c - 4 * b * d)) / (2 * b);
95  N = 2;
96  return;
97  }
98  }
99 
100  // Ok, it is really a cubic
101 
102  // Discriminant
103  double DELTA = 18 * a * b * c * d - 4 * b * b * b * d + b * b * c * c - 4 * a * c * c * c - 27 * a * a * d * d;
104  // Coefficients for the depressed cubic t^3+p*t+q = 0
105  double p = (3 * a * c - b * b) / (3 * a * a);
106  double q = (2 * b * b * b - 9 * a * b * c + 27 * a * a * d) / (27 * a * a * a);
107 
108  if (DELTA < 0) {
109  // One real root
110  double t0;
111  if (4 * p * p * p + 27 * q * q > 0 && p < 0) {
112  t0 = -2.0 * std::abs(q) / q * sqrt(-p / 3.0) * cosh(1.0 / 3.0 * acosh(-3.0 * std::abs(q) / (2.0 * p) * sqrt(-3.0 / p)));
113  } else {
114  t0 = -2.0 * sqrt(p / 3.0) * sinh(1.0 / 3.0 * asinh(3.0 * q / (2.0 * p) * sqrt(3.0 / p)));
115  }
116  N = 1;
117  x0 = t0 - b / (3 * a);
118  x1 = t0 - b / (3 * a);
119  x2 = t0 - b / (3 * a);
120  } else //(DELTA>0)
121  {
122  // Three real roots
123  double t0 = 2.0 * sqrt(-p / 3.0) * cos(1.0 / 3.0 * acos(3.0 * q / (2.0 * p) * sqrt(-3.0 / p)) - 0 * 2.0 * M_PI / 3.0);
124  double t1 = 2.0 * sqrt(-p / 3.0) * cos(1.0 / 3.0 * acos(3.0 * q / (2.0 * p) * sqrt(-3.0 / p)) - 1 * 2.0 * M_PI / 3.0);
125  double t2 = 2.0 * sqrt(-p / 3.0) * cos(1.0 / 3.0 * acos(3.0 * q / (2.0 * p) * sqrt(-3.0 / p)) - 2 * 2.0 * M_PI / 3.0);
126 
127  N = 3;
128  x0 = t0 - b / (3 * a);
129  x1 = t1 - b / (3 * a);
130  x2 = t2 - b / (3 * a);
131  }
132 }
133 void solve_quartic(double a, double b, double c, double d, double e, int& N, double& x0, double& x1, double& x2, double& x3) {
134 
135  // 0 = ax^4 + b*x^3 + c*x^2 + d*x + e
136 
137  Eigen::PolynomialSolver<double, Eigen::Dynamic> solver;
138  Eigen::VectorXd coeff(5);
139  coeff << e, d, c, b, a;
140  solver.compute(coeff);
141 
142  std::vector<double> realRoots;
143  solver.realRoots(realRoots);
144  N = static_cast<int>(realRoots.size());
145 
146  if (N > 0) {
147  x0 = realRoots[0];
148  }
149  if (N > 1) {
150  x1 = realRoots[1];
151  }
152  if (N > 2) {
153  x2 = realRoots[2];
154  }
155  if (N > 3) {
156  x3 = realRoots[3];
157  }
158 }
159 
161  if (Nconstraints == 4) {
162  std::vector<double> abcd = CoolProp::linsolve(A, B);
163  a = abcd[0];
164  b = abcd[1];
165  c = abcd[2];
166  d = abcd[3];
167  return true;
168  } else {
169  throw CoolProp::ValueError(format("Number of constraints[%d] is not equal to 4", Nconstraints));
170  }
171 }
172 bool SplineClass::add_value_constraint(double x, double y) {
173  int i = Nconstraints;
174  if (i == 4) return false;
175  A[i][0] = x * x * x;
176  A[i][1] = x * x;
177  A[i][2] = x;
178  A[i][3] = 1;
179  B[i] = y;
180  Nconstraints++;
181  return true;
182 }
183 void SplineClass::add_4value_constraints(double x1, double x2, double x3, double x4, double y1, double y2, double y3, double y4) {
184  add_value_constraint(x1, y1);
185  add_value_constraint(x2, y2);
186  add_value_constraint(x3, y3);
187  add_value_constraint(x4, y4);
188 }
189 bool SplineClass::add_derivative_constraint(double x, double dydx) {
190  int i = Nconstraints;
191  if (i == 4) return false;
192  A[i][0] = 3 * x * x;
193  A[i][1] = 2 * x;
194  A[i][2] = 1;
195  A[i][3] = 0;
196  B[i] = dydx;
197  Nconstraints++;
198  return true;
199 }
200 double SplineClass::evaluate(double x) {
201  return a * x * x * x + b * x * x + c * x + d;
202 }