CPnumerics.cpp
#include "CPnumerics.h"
#include "MatrixMath.h"
#include <unsupported/Eigen/Polynomials>
double root_sum_square(const std::vector<double>& x) {
double sum = 0;
for (unsigned int i = 0; i < x.size(); i++) {
sum += pow(x[i], 2);
}
return sqrt(sum);
}
double interp1d(const std::vector<double>* x, const std::vector<double>* y, double x0) {
std::size_t i, L, R, M;
L = 0;
R = (*x).size() - 1;
M = (L + R) / 2;
// Use interval halving to find the indices which bracket the density of interest
while (R - L > 1) {
if (x0 >= (*x)[M]) {
L = M;
M = (L + R) / 2;
continue;
}
if (x0 < (*x)[M]) {
R = M;
M = (L + R) / 2;
continue;
}
}
i = L;
if (i < (*x).size() - 2) {
// Go "forwards" with the interpolation range
return QuadInterp((*x)[i], (*x)[i + 1], (*x)[i + 2], (*y)[i], (*y)[i + 1], (*y)[i + 2], x0);
} else {
// Go "backwards" with the interpolation range
return QuadInterp((*x)[i], (*x)[i - 1], (*x)[i - 2], (*y)[i], (*y)[i - 1], (*y)[i - 2], x0);
}
}
double powInt(double x, int y) {
// Raise a double to an integer power
// Overload not provided in math.h
int i;
double product = 1.0;
double x_in;
int y_in;
if (y == 0) {
return 1.0;
}
if (y < 0) {
x_in = 1 / x;
y_in = -y;
} else {
x_in = x;
y_in = y;
}
if (y_in == 1) {
return x_in;
}
product = x_in;
for (i = 1; i < y_in; i++) {
product = product * x_in;
}
return product;
}
void MatInv_2(double A[2][2], double B[2][2]) {
double Det;
//Using Cramer's Rule to solve
Det = A[0][0] * A[1][1] - A[1][0] * A[0][1];
B[0][0] = 1.0 / Det * A[1][1];
B[1][1] = 1.0 / Det * A[0][0];
B[1][0] = -1.0 / Det * A[1][0];
B[0][1] = -1.0 / Det * A[0][1];
}
void solve_cubic(double a, double b, double c, double d, int& N, double& x0, double& x1, double& x2) {
// 0 = ax^3 + b*x^2 + c*x + d
// First check if the "cubic" is actually a second order or first order curve
if (std::abs(a) < 10 * DBL_EPSILON) {
if (std::abs(b) < 10 * DBL_EPSILON) {
// Linear solution if a = 0 and b = 0
x0 = -d / c;
N = 1;
return;
} else {
// Quadratic solution(s) if a = 0 and b != 0
x0 = (-c + sqrt(c * c - 4 * b * d)) / (2 * b);
x1 = (-c - sqrt(c * c - 4 * b * d)) / (2 * b);
N = 2;
return;
}
}
// Ok, it is really a cubic
// Discriminant
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;
// Coefficients for the depressed cubic t^3+p*t+q = 0
double p = (3 * a * c - b * b) / (3 * a * a);
double q = (2 * b * b * b - 9 * a * b * c + 27 * a * a * d) / (27 * a * a * a);
if (DELTA < 0) {
// One real root
double t0;
if (4 * p * p * p + 27 * q * q > 0 && p < 0) {
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)));
} else {
t0 = -2.0 * sqrt(p / 3.0) * sinh(1.0 / 3.0 * asinh(3.0 * q / (2.0 * p) * sqrt(3.0 / p)));
}
N = 1;
x0 = t0 - b / (3 * a);
x1 = t0 - b / (3 * a);
x2 = t0 - b / (3 * a);
} else //(DELTA>0)
{
// Three real roots
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);
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);
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);
N = 3;
x0 = t0 - b / (3 * a);
x1 = t1 - b / (3 * a);
x2 = t2 - b / (3 * a);
}
}
void solve_quartic(double a, double b, double c, double d, double e, int& N, double& x0, double& x1, double& x2, double& x3) {
// 0 = ax^4 + b*x^3 + c*x^2 + d*x + e
Eigen::PolynomialSolver<double, Eigen::Dynamic> solver;
Eigen::VectorXd coeff(5);
coeff << e, d, c, b, a;
solver.compute(coeff);
std::vector<double> realRoots;
solver.realRoots(realRoots);
N = static_cast<int>(realRoots.size());
if (N > 0) {
x0 = realRoots[0];
}
if (N > 1) {
x1 = realRoots[1];
}
if (N > 2) {
x2 = realRoots[2];
}
if (N > 3) {
x3 = realRoots[3];
}
}
bool SplineClass::build() {
if (Nconstraints == 4) {
std::vector<double> abcd = CoolProp::linsolve(A, B);
a = abcd[0];
b = abcd[1];
c = abcd[2];
d = abcd[3];
return true;
} else {
throw CoolProp::ValueError(format("Number of constraints[%d] is not equal to 4", Nconstraints));
}
}
bool SplineClass::add_value_constraint(double x, double y) {
int i = Nconstraints;
if (i == 4) return false;
A[i][0] = x * x * x;
A[i][1] = x * x;
A[i][2] = x;
A[i][3] = 1;
B[i] = y;
Nconstraints++;
return true;
}
void SplineClass::add_4value_constraints(double x1, double x2, double x3, double x4, double y1, double y2, double y3, double y4) {
add_value_constraint(x1, y1);
add_value_constraint(x2, y2);
add_value_constraint(x3, y3);
add_value_constraint(x4, y4);
}
bool SplineClass::add_derivative_constraint(double x, double dydx) {
int i = Nconstraints;
if (i == 4) return false;
A[i][0] = 3 * x * x;
A[i][1] = 2 * x;
A[i][2] = 1;
A[i][3] = 0;
B[i] = dydx;
Nconstraints++;
return true;
}
double SplineClass::evaluate(double x) {
return a * x * x * x + b * x * x + c * x + d;
}