CoolProp 8.0.0
An open-source fluid property and humid air property database
GeneralizedCubic.cpp
Go to the documentation of this file.
1#include "GeneralizedCubic.h"
2#include "CoolProp/numerics/numerics.h"
3#include <array>
4#include <cassert>
5#include <cmath>
6#include <cstring>
7#include <limits>
8#include <utility>
9
10double BasicMathiasCopemanAlphaFunction::term(double tau, std::size_t itau) {
11
12 // If we are not using the full Mathias-Copeman formulation for a_ii,
13 // we just use the simple results from the supplemental information because
14 // they are much more computationally efficient
15
16 // All derivatives have a common bracketed term, so we factor it out
17 // and calculate it here
18 double B = 1 + m * (1 - sqrt_Tr_Tci * sqrt(1 / tau));
19
20 switch (itau) {
21 case 0:
22 return a0 * B * B;
23 case 1:
24 return a0 * m * B / pow(tau, 3.0 / 2.0) * sqrt_Tr_Tci;
25 case 2:
26 return a0 * m / 2.0 * (m / pow(tau, 3) * Tr_over_Tci - 3 * B / pow(tau, 5.0 / 2.0) * sqrt_Tr_Tci);
27 case 3:
28 return (3.0 / 4.0) * a0 * m * (-3.0 * m / pow(tau, 4) * Tr_over_Tci + 5 * B / pow(tau, 7.0 / 2.0) * sqrt_Tr_Tci);
29 case 4:
30 return (3.0 / 8.0) * a0 * m * (29.0 * m / pow(tau, 5) * Tr_over_Tci - 35 * B / pow(tau, 9.0 / 2.0) * sqrt_Tr_Tci);
31 default:
32 throw -1;
33 }
34}
35
36double MathiasCopemanAlphaFunction::term(double tau, std::size_t itau) {
37 // Here we are using the full Mathias-Copeman formulation, introducing
38 // some additional computational effort, so we only evaluate the parameters that
39 // we actually need to evaluate, otherwise we just set their value to zero
40 // See info on the conditional (ternary) operator : http://www.cplusplus.com/articles/1AUq5Di1/
41 // Furthermore, this should help with branch prediction
42 double Di = 1 - sqrt_Tr_Tci / sqrt(tau);
43 double dDi_dtau = (itau >= 1) ? (1.0 / 2.0) * sqrt_Tr_Tci / (pow(tau, 1.5)) : 0;
44 double d2Di_dtau2 = (itau >= 2) ? -(3.0 / 4.0) * sqrt_Tr_Tci / (pow(tau, 2.5)) : 0;
45 double d3Di_dtau3 = (itau >= 3) ? (15.0 / 8.0) * sqrt_Tr_Tci / (pow(tau, 3.5)) : 0;
46 double d4Di_dtau4 = (itau >= 4) ? -(105.0 / 16.0) * sqrt_Tr_Tci / (pow(tau, 4.5)) : 0;
47
48 double Bi = 1, dBi_dtau = 0, d2Bi_dtau2 = 0, d3Bi_dtau3 = 0, d4Bi_dtau4 = 0;
49 for (int n = 1; n <= 3; ++n) {
50 Bi += c[n - 1] * pow(Di, n);
51 dBi_dtau += (itau < 1) ? 0 : (n * c[n - 1] * pow(Di, n - 1) * dDi_dtau);
52 d2Bi_dtau2 += (itau < 2) ? 0 : n * c[n - 1] * ((n - 1) * pow(dDi_dtau, 2) + Di * d2Di_dtau2) * pow(Di, n - 2);
53 d3Bi_dtau3 += (itau < 3)
54 ? 0
55 : n * c[n - 1] * (3 * (n - 1) * Di * dDi_dtau * d2Di_dtau2 + (n * n - 3 * n + 2) * pow(dDi_dtau, 3) + pow(Di, 2) * d3Di_dtau3)
56 * pow(Di, n - 3);
57 d4Bi_dtau4 +=
58 (itau < 4)
59 ? 0
60 : n * c[n - 1]
61 * (6 * (n * n - 3 * n + 2) * Di * pow(dDi_dtau, 2) * d2Di_dtau2 + (n * n * n - 6 * n * n + 11 * n - 6) * pow(dDi_dtau, 4)
62 + (4 * n * dDi_dtau * d3Di_dtau3 + 3 * n * pow(d2Di_dtau2, 2) - 4 * dDi_dtau * d3Di_dtau3 - 3 * pow(d2Di_dtau2, 2)) * pow(Di, 2)
63 + pow(Di, 3) * d4Di_dtau4)
64 * pow(Di, n - 4);
65 }
66 switch (itau) {
67 case 0:
68 return a0 * Bi * Bi;
69 case 1:
70 return 2 * a0 * Bi * dBi_dtau;
71 case 2:
72 return 2 * a0 * (Bi * d2Bi_dtau2 + dBi_dtau * dBi_dtau);
73 case 3:
74 return 2 * a0 * (Bi * d3Bi_dtau3 + 3 * dBi_dtau * d2Bi_dtau2);
75 case 4:
76 return 2 * a0 * (Bi * d4Bi_dtau4 + 4 * dBi_dtau * d3Bi_dtau3 + 3 * pow(d2Bi_dtau2, 2));
77 default:
78 throw -1;
79 }
80}
81
82double TwuAlphaFunction::term(double tau, std::size_t itau) {
83 // Here we are using the Twu formulation, introducing
84 // some additional computational effort, so we only evaluate the parameters that
85 // we actually need to evaluate, otherwise we just set their value to zero
86 // See info on the conditional (ternary) operator : http://www.cplusplus.com/articles/1AUq5Di1/
87 // Furthermore, this should help with branch prediction
88
89 const double L = c[0], M = c[1], N = c[2];
90 double A = pow(Tr_over_Tci / tau, M * N);
91 double B1 = (itau < 1) ? 0 : N / tau * (L * M * A - M + 1);
92 double dB1_dtau = (itau < 2) ? 0 : N / powInt(tau, 2) * (-L * M * M * N * A - L * M * A + M - 1);
93 double d2B1_dtau2 = (itau < 3) ? 0 : N / powInt(tau, 3) * (L * M * M * M * N * N * A + 3 * L * M * M * N * A + 2 * L * M * A - 2 * M + 2);
94 double d3B1_dtau3 =
95 (itau < 4) ? 0
96 : -N / powInt(tau, 4)
97 * (L * powInt(M, 4) * powInt(N, 3) * A + 6 * L * M * M * M * N * N * A + 11 * L * M * M * N * A + 6 * L * M * A - 6 * M + 6);
98
99 double dam_dtau = NAN, d2am_dtau2 = NAN, d3am_dtau3 = NAN, d4am_dtau4 = NAN;
100 double am = a0 * pow(Tr_over_Tci / tau, N * (M - 1)) * exp(L * (1 - A));
101
102 if (itau == 0) {
103 return am;
104 } else {
105 // Calculate terms as needed
106 dam_dtau = am * B1;
107 d2am_dtau2 = (itau < 2) ? 0 : B1 * dam_dtau + am * dB1_dtau;
108 d3am_dtau3 = (itau < 3) ? 0 : B1 * d2am_dtau2 + am * d2B1_dtau2 + 2 * dB1_dtau * dam_dtau;
109 d4am_dtau4 = (itau < 4) ? 0 : B1 * d3am_dtau3 + am * d3B1_dtau3 + 3 * dB1_dtau * d2am_dtau2 + 3 * d2B1_dtau2 * dam_dtau;
110 }
111 switch (itau) {
112 case 1:
113 return dam_dtau;
114 case 2:
115 return d2am_dtau2;
116 case 3:
117 return d3am_dtau3;
118 case 4:
119 return d4am_dtau4;
120 default:
121 throw -1;
122 }
123}
124
125void BasicMathiasCopemanAlphaFunction::calc_all_terms(double tau, std::array<double, 5>& terms) {
126 // Compute B and shared powers of tau once, then fill all 5 derivatives
127 const double sq = sqrt_Tr_Tci / sqrt(tau); // sqrt(Tr/Tci / tau)
128 const double B = 1.0 + m * (1.0 - sq);
129 const double t15 = sqrt_Tr_Tci / (tau * sqrt(tau)); // sqrt_Tr_Tci / tau^1.5
130 const double t25 = t15 / tau; // sqrt_Tr_Tci / tau^2.5
131 const double t35 = t25 / tau; // sqrt_Tr_Tci / tau^3.5
132 const double t45 = t35 / tau; // sqrt_Tr_Tci / tau^4.5
133 const double t3 = Tr_over_Tci / (tau * tau * tau);
134 const double t4 = t3 / tau;
135 const double t5 = t4 / tau;
136 terms[0] = a0 * B * B;
137 terms[1] = a0 * m * B * t15;
138 terms[2] = a0 * m * 0.5 * (m * t3 - 3.0 * B * t25);
139 terms[3] = (3.0 / 4.0) * a0 * m * (-3.0 * m * t4 + 5.0 * B * t35);
140 terms[4] = (3.0 / 8.0) * a0 * m * (29.0 * m * t5 - 35.0 * B * t45);
141}
142
143void TwuAlphaFunction::calc_all_terms(double tau, std::array<double, 5>& terms) {
144 const double L = c[0], M = c[1], N = c[2];
145 // Compute the two pow() and one exp() calls exactly once
146 const double A = pow(Tr_over_Tci / tau, M * N);
147 const double am = a0 * pow(Tr_over_Tci / tau, N * (M - 1)) * exp(L * (1.0 - A));
148 terms[0] = am;
149 const double B1 = N / tau * (L * M * A - M + 1.0);
150 const double dam = am * B1;
151 terms[1] = dam;
152 const double dB1 = N / powInt(tau, 2) * (-L * M * M * N * A - L * M * A + M - 1.0);
153 const double d2am = B1 * dam + am * dB1;
154 terms[2] = d2am;
155 const double d2B1 = N / powInt(tau, 3) * (L * M * M * M * N * N * A + 3.0 * L * M * M * N * A + 2.0 * L * M * A - 2.0 * M + 2.0);
156 const double d3am = B1 * d2am + am * d2B1 + 2.0 * dB1 * dam;
157 terms[3] = d3am;
158 const double d3B1 =
159 -N / powInt(tau, 4)
160 * (L * powInt(M, 4) * powInt(N, 3) * A + 6.0 * L * M * M * M * N * N * A + 11.0 * L * M * M * N * A + 6.0 * L * M * A - 6.0 * M + 6.0);
161 terms[4] = B1 * d3am + am * d3B1 + 3.0 * dB1 * d2am + 3.0 * d2B1 * dam;
162}
163
164void MathiasCopemanAlphaFunction::calc_all_terms(double tau, std::array<double, 5>& terms) {
165 // Compute Di = 1 - sqrt(Tr/Tci / tau) and its tau-derivatives once, reusing
166 // a common power chain. All five derivatives share these values.
167 const double sq = sqrt_Tr_Tci / sqrt(tau); // sqrt(Tr_over_Tci / tau)
168 const double Di = 1.0 - sq;
169 const double p = sq / tau; // sqrt_Tr_Tci / tau^1.5
170 const double d1Di = 0.5 * p;
171 const double d2Di = -0.75 * p / tau;
172 const double d3Di = 1.875 * p / (tau * tau);
173 const double d4Di = -6.5625 * p / (tau * tau * tau);
174
175 const double Di2 = Di * Di;
176 const double Di3 = Di2 * Di;
177 const double c0 = c[0], c1 = c[1], c2 = c[2];
178
179 // Expand the loop over n=1,2,3 analytically.
180 // Many polynomial factors in (n^2-3n+2) and (n^3-6n^2+11n-6) vanish for small n.
181 const double Bi = 1.0 + c0 * Di + c1 * Di2 + c2 * Di3;
182 const double dBi = d1Di * (c0 + 2.0 * c1 * Di + 3.0 * c2 * Di2);
183 const double d2Bi = c0 * d2Di + 2.0 * c1 * (d1Di * d1Di + Di * d2Di) + 3.0 * c2 * Di * (2.0 * d1Di * d1Di + Di * d2Di);
184 const double d3Bi =
185 c0 * d3Di + 2.0 * c1 * (3.0 * d1Di * d2Di + Di * d3Di) + 3.0 * c2 * (6.0 * Di * d1Di * d2Di + 2.0 * d1Di * d1Di * d1Di + Di2 * d3Di);
186 const double d4Bi = c0 * d4Di + 2.0 * c1 * (4.0 * d1Di * d3Di + 3.0 * d2Di * d2Di + Di * d4Di)
187 + 3.0 * c2 * (12.0 * d1Di * d1Di * d2Di + Di * (8.0 * d1Di * d3Di + 6.0 * d2Di * d2Di) + Di2 * d4Di);
188
189 terms[0] = a0 * Bi * Bi;
190 terms[1] = 2.0 * a0 * Bi * dBi;
191 terms[2] = 2.0 * a0 * (Bi * d2Bi + dBi * dBi);
192 terms[3] = 2.0 * a0 * (Bi * d3Bi + 3.0 * dBi * d2Bi);
193 terms[4] = 2.0 * a0 * (Bi * d4Bi + 4.0 * dBi * d3Bi + 3.0 * d2Bi * d2Bi);
194}
195
196AbstractCubic::AbstractCubic(const std::vector<double>& Tc, std::vector<double> pc, std::vector<double> acentric, double R_u, double Delta_1,
197 double Delta_2, const std::vector<double>& C1, const std::vector<double>& C2, const std::vector<double>& C3)
198 : T_r(1.0),
199 rho_r(1.0),
200 Tc(Tc),
201 pc(std::move(pc)),
202 acentric(std::move(acentric)),
203 R_u(R_u),
204 Delta_1(Delta_1),
205 Delta_2(Delta_2),
206 N(static_cast<int>(Tc.size())),
207 cm(0.) {
208
209 k.resize(N, std::vector<double>(N, 0));
210
211 alpha.resize(N);
212 m_tau_cache = std::numeric_limits<double>::quiet_NaN();
213};
214
215void AbstractCubic::set_alpha(const std::vector<double>& C1, const std::vector<double>& C2, const std::vector<double>& C3) {
217 alpha.resize(Tc.size());
218 m_tau_cache = std::numeric_limits<double>::quiet_NaN(); // invalidate cache
220 if (C1.empty() && C2.empty() && C3.empty()) {
221 for (std::size_t i = 0; i < Tc.size(); ++i) {
222 alpha[i] = std::make_shared<BasicMathiasCopemanAlphaFunction>(a0_ii(i), m_ii(i), T_r / Tc[i]);
223 }
224 } else {
226 for (std::size_t i = 0; i < Tc.size(); ++i) {
227 alpha[i] = std::make_shared<MathiasCopemanAlphaFunction>(a0_ii(i), C1[i], C2[i], C3[i], T_r / Tc[i]);
228 }
229 }
230}
231
232double AbstractCubic::am_term(double tau, const std::vector<double>& x, std::size_t itau) {
233 // am_term and its composition-derivatives are the only callers of aij_term/u_term, so we
234 // validate/populate the aii cache for this tau once here. The inner u_term reads then go
235 // straight to m_aii_cache instead of re-validating the cache on every aii lookup, which was
236 // the dominant cost (~38% of self-time) in the pure-fluid hot path.
237 _ensure_aii_cache(tau);
238 double summer = 0;
239 for (int i = N - 1; i >= 0; --i) {
240 for (int j = N - 1; j >= 0; --j) {
241 summer += x[i] * x[j] * aij_term(tau, i, j, itau);
242 }
243 }
244 return summer;
245}
246double AbstractCubic::d_am_term_dxi(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, bool xN_independent) {
247 _ensure_aii_cache(tau); // see am_term: guarantees u_term can read m_aii_cache directly
248 if (xN_independent) {
249 double summer = 0;
250 for (int j = N - 1; j >= 0; --j) {
251 summer += x[j] * aij_term(tau, i, j, itau);
252 }
253 return 2 * summer;
254 } else {
255 double summer = 0;
256 for (int k = N - 2; k >= 0; --k) {
257 summer += x[k] * (aij_term(tau, i, k, itau) - aij_term(tau, k, N - 1, itau));
258 }
259 return 2 * (summer + x[N - 1] * (aij_term(tau, N - 1, i, itau) - aij_term(tau, N - 1, N - 1, itau)));
260 }
261}
262double AbstractCubic::d2_am_term_dxidxj(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, std::size_t j,
263 bool xN_independent) {
264 _ensure_aii_cache(tau); // see am_term: guarantees u_term can read m_aii_cache directly
265 if (xN_independent) {
266 return 2 * aij_term(tau, i, j, itau);
267 } else {
268 return 2 * (aij_term(tau, i, j, itau) - aij_term(tau, j, N - 1, itau) - aij_term(tau, N - 1, i, itau) + aij_term(tau, N - 1, N - 1, itau));
269 }
270}
271
272double AbstractCubic::d3_am_term_dxidxjdxk(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, std::size_t j, std::size_t k,
273 bool xN_independent) {
274 return 0;
275}
276
277double AbstractCubic::bm_term(const std::vector<double>& x) {
278 // b0_ii(i) is a constant (R_u*Tc[i]/pc[i] scaled); cache it once so this sum -- called many
279 // times per all() via psi_minus/psi_plus/PI_12 -- avoids the repeated division. Bit-exact:
280 // m_b0_ii_cache[i] holds exactly the value b0_ii(i) returns.
281 if (m_b0_ii_cache.size() != static_cast<std::size_t>(N)) {
282 m_b0_ii_cache.resize(N);
283 for (int i = 0; i < N; ++i) {
284 m_b0_ii_cache[i] = b0_ii(i);
285 }
286 }
287 double summer = 0;
288 for (int i = N - 1; i >= 0; --i) {
289 summer += x[i] * m_b0_ii_cache[i];
290 }
291 return summer;
292}
293double AbstractCubic::d_bm_term_dxi(const std::vector<double>& x, std::size_t i, bool xN_independent) {
294 if (xN_independent) {
295 return b0_ii(i);
296 } else {
297 return b0_ii(i) - b0_ii(N - 1);
298 }
299}
300double AbstractCubic::d2_bm_term_dxidxj(const std::vector<double>& x, std::size_t i, std::size_t j, bool xN_independent) {
301 return 0;
302}
303double AbstractCubic::d3_bm_term_dxidxjdxk(const std::vector<double>& x, std::size_t i, std::size_t j, std::size_t k, bool xN_independent) {
304 return 0;
305}
306
307double AbstractCubic::cm_term() {
308 return cm;
309}
310
311double AbstractCubic::aii_term(double tau, std::size_t i, std::size_t itau) {
312 if (itau > 4) {
313 // m_aii_cache[i] only stores derivatives for itau = 0..4 (see calc_all_terms()); higher
314 // orders are unsupported, matching the contract of AbstractCubicAlphaFunction::term().
315 throw -1;
316 }
317 _ensure_aii_cache(tau);
318 return m_aii_cache[i][itau];
319}
320double AbstractCubic::u_term(double tau, std::size_t i, std::size_t j, std::size_t itau) {
321 // Read the aii derivatives straight from the cache. The only callers (the am_term family)
322 // populate it for this tau first; the assert pins that invariant in debug builds so a future
323 // caller that forgets to do so fails loudly instead of silently returning stale values.
324 assert(std::memcmp(&tau, &m_tau_cache, sizeof(double)) == 0 && "u_term: aii cache not populated for this tau");
325 (void)tau;
326 const std::array<double, 5>& ai = m_aii_cache[i];
327 const std::array<double, 5>& aj = m_aii_cache[j];
328 const double aii = ai[0], ajj = aj[0];
329 switch (itau) {
330 case 0:
331 return aii * ajj;
332 case 1:
333 return aii * aj[1] + ajj * ai[1];
334 case 2:
335 return (aii * aj[2] + 2 * ai[1] * aj[1] + ajj * ai[2]);
336 case 3:
337 return (aii * aj[3] + 3 * ai[1] * aj[2] + 3 * ai[2] * aj[1] + ajj * ai[3]);
338 case 4:
339 return (aii * aj[4] + 4 * ai[1] * aj[3] + 6 * ai[2] * aj[2] + 4 * ai[3] * aj[1] + ajj * ai[4]);
340 default:
341 throw -1;
342 }
343}
344double AbstractCubic::aij_term(double tau, std::size_t i, std::size_t j, std::size_t itau) {
345 // Half-integer powers of u are written as u^k * sqrt(u), and integer powers as repeated
346 // multiplication, to avoid std::pow (the dominant cost of the residual hot path). u_term
347 // values are pulled into locals so each is computed once. Results differ from the previous
348 // std::pow formulation only at the last ULP.
349 const double u = u_term(tau, i, j, 0);
350 const double su = sqrt(u);
351 const double kf = 1 - k[i][j];
352
353 switch (itau) {
354 case 0:
355 return kf * su;
356 case 1:
357 return kf / (2.0 * su) * u_term(tau, i, j, 1);
358 case 2: {
359 const double u1 = u_term(tau, i, j, 1);
360 return kf / (4.0 * u * su) * (2 * u * u_term(tau, i, j, 2) - u1 * u1);
361 }
362 case 3: {
363 const double u1 = u_term(tau, i, j, 1), u2 = u_term(tau, i, j, 2);
364 return kf / (8.0 * u * u * su) * (4 * u * u * u_term(tau, i, j, 3) - 6 * u * u1 * u2 + 3 * u1 * u1 * u1);
365 }
366 case 4: {
367 const double u1 = u_term(tau, i, j, 1), u2 = u_term(tau, i, j, 2);
368 return kf / (16.0 * u * u * u * su)
369 * (-4 * u * u * (4 * u1 * u_term(tau, i, j, 3) + 3 * u2 * u2) + 8 * u * u * u * u_term(tau, i, j, 4) + 36 * u * u1 * u1 * u2
370 - 15 * u1 * u1 * u1 * u1);
371 }
372 default:
373 throw -1;
374 }
375}
376double AbstractCubic::psi_minus(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta) {
377 if (itau > 0) return 0.0;
378 double bmc = bm_term(x) - cm_term(); // appears only in the form (b-c) in the equations
379 double bracket = 1 - bmc * delta * rho_r;
380
381 switch (idelta) {
382 case 0:
383 return -log(bracket);
384 case 1:
385 return bmc * rho_r / bracket;
386 case 2: {
387 const double r = bmc * rho_r / bracket;
388 return r * r;
389 }
390 case 3: {
391 const double r = bmc * rho_r / bracket;
392 return 2 * r * r * r;
393 }
394 case 4: {
395 const double r = bmc * rho_r / bracket;
396 return 6 * r * r * r * r;
397 }
398 default:
399 throw -1;
400 }
401}
402double AbstractCubic::d_psi_minus_dxi(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
403 bool xN_independent) {
404 if (itau > 0) return 0.0;
405 double bmc = bm_term(x) - cm_term(); // appears only in the form (b-c) in the equations
406 double db_dxi = d_bm_term_dxi(x, i, xN_independent);
407 double bracket = 1 - bmc * delta * rho_r;
408
409 switch (idelta) {
410 case 0:
411 return delta * rho_r * db_dxi / bracket;
412 case 1:
413 return rho_r * db_dxi / pow(bracket, 2);
414 case 2:
415 return 2 * pow(rho_r, 2) * bmc * db_dxi / pow(bracket, 3);
416 case 3:
417 return 6 * pow(rho_r, 3) * pow(bmc, 2) * db_dxi / pow(bracket, 4);
418 case 4:
419 return 24 * pow(rho_r, 4) * pow(bmc, 3) * db_dxi / pow(bracket, 5);
420 default:
421 throw -1;
422 }
423}
424double AbstractCubic::d2_psi_minus_dxidxj(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
425 std::size_t j, bool xN_independent) {
426 if (itau > 0) return 0.0;
427 double bmc = bm_term(x) - cm_term(); // appears only in the form (b-c) in the equations
428 double db_dxi = d_bm_term_dxi(x, i, xN_independent), db_dxj = d_bm_term_dxi(x, j, xN_independent),
429 d2b_dxidxj = d2_bm_term_dxidxj(x, i, j, xN_independent);
430 double bracket = 1 - bmc * delta * rho_r;
431
432 switch (idelta) {
433 case 0:
434 return pow(delta * rho_r, 2) * db_dxi * db_dxj / pow(bracket, 2) + delta * rho_r * d2b_dxidxj / bracket;
435 case 1:
436 return 2 * delta * pow(rho_r, 2) * db_dxi * db_dxj / pow(bracket, 3) + rho_r * d2b_dxidxj / pow(bracket, 2);
437 case 2:
438 return 2 * pow(rho_r, 2) * db_dxi * db_dxj / pow(bracket, 4) * (2 * delta * rho_r * bmc + 1)
439 + 2 * pow(rho_r, 2) * bmc * d2b_dxidxj / pow(bracket, 3);
440 case 3:
441 return 12 * pow(rho_r, 3) * bmc * db_dxi * db_dxj / pow(bracket, 5) * (delta * rho_r * bmc + 1)
442 + 6 * pow(rho_r, 3) * pow(bmc, 2) * d2b_dxidxj / pow(bracket, 4);
443 case 4:
444 return 24 * pow(rho_r, 4) * pow(bmc, 2) * db_dxi * db_dxj / pow(bracket, 6) * (2 * delta * rho_r * bmc + 3)
445 + 24 * pow(rho_r, 4) * pow(bmc, 3) * d2b_dxidxj / pow(bracket, 5);
446 default:
447 throw -1;
448 }
449}
450double AbstractCubic::d3_psi_minus_dxidxjdxk(double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
451 std::size_t j, std::size_t k, bool xN_independent) {
452 if (itau > 0) return 0.0;
453 double bmc = bm_term(x) - cm_term(); // appears only in the form (b-c) in the equations
454 double db_dxi = d_bm_term_dxi(x, i, xN_independent), db_dxj = d_bm_term_dxi(x, j, xN_independent), db_dxk = d_bm_term_dxi(x, k, xN_independent),
455 d2b_dxidxj = d2_bm_term_dxidxj(x, i, j, xN_independent), d2b_dxidxk = d2_bm_term_dxidxj(x, i, k, xN_independent),
456 d2b_dxjdxk = d2_bm_term_dxidxj(x, j, k, xN_independent), d3b_dxidxjdxk = d3_bm_term_dxidxjdxk(x, i, j, k, xN_independent);
457 double bracket = 1 - bmc * delta * rho_r;
458
459 switch (idelta) {
460 case 0:
461 return delta * rho_r * d3b_dxidxjdxk / bracket + 2 * pow(delta * rho_r, 3) * db_dxi * db_dxj * db_dxk / pow(bracket, 3)
462 + pow(delta * rho_r, 2) / pow(bracket, 2) * (db_dxi * d2b_dxjdxk + db_dxj * d2b_dxidxk + db_dxk * d2b_dxidxj);
463 case 1:
464 return rho_r * d3b_dxidxjdxk / pow(bracket, 2) + 6 * pow(delta, 2) * pow(rho_r, 3) * db_dxi * db_dxj * db_dxk / pow(bracket, 4)
465 + 2 * delta * pow(rho_r, 2) / pow(bracket, 3) * (db_dxi * d2b_dxjdxk + db_dxj * d2b_dxidxk + db_dxk * d2b_dxidxj);
466 default:
467 throw -1;
468 }
469}
470double AbstractCubic::PI_12(double delta, const std::vector<double>& x, std::size_t idelta) {
471 double bm = bm_term(x);
472 double cm = cm_term();
473 switch (idelta) {
474 case 0:
475 return (1 + (Delta_1 * bm + cm) * rho_r * delta) * (1 + (Delta_2 * bm + cm) * rho_r * delta);
476 case 1:
477 return rho_r * (2 * (bm * Delta_1 + cm) * (bm * Delta_2 + cm) * delta * rho_r + (Delta_1 + Delta_2) * bm + 2 * cm);
478 case 2:
479 return 2 * (Delta_1 * bm + cm) * (Delta_2 * bm + cm) * rho_r * rho_r;
480 case 3:
481 return 0;
482 case 4:
483 return 0;
484 default:
485 throw -1;
486 }
487}
488double AbstractCubic::d_PI_12_dxi(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, bool xN_independent) {
489 double bm = bm_term(x);
490 double cm = cm_term();
491 double db_dxi = d_bm_term_dxi(x, i, xN_independent);
492 switch (idelta) {
493 case 0:
494 return delta * rho_r * db_dxi * (2 * Delta_1 * Delta_2 * bm * delta * rho_r + (Delta_1 + Delta_2) * (1 + cm * delta * rho_r));
495 case 1:
496 return rho_r * db_dxi * (4 * Delta_1 * Delta_2 * bm * delta * rho_r + (Delta_1 + Delta_2) * (1 + 2 * cm * delta * rho_r));
497 case 2:
498 return 2 * pow(rho_r, 2) * (2 * Delta_1 * Delta_2 * bm + (Delta_1 + Delta_2) * cm) * db_dxi;
499 case 3:
500 return 0;
501 case 4:
502 return 0;
503 default:
504 throw -1;
505 }
506}
507double AbstractCubic::d2_PI_12_dxidxj(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, std::size_t j,
508 bool xN_independent) {
509 double bm = bm_term(x);
510 double cm = cm_term();
511 double db_dxi = d_bm_term_dxi(x, i, xN_independent), db_dxj = d_bm_term_dxi(x, j, xN_independent),
512 d2b_dxidxj = d2_bm_term_dxidxj(x, i, j, xN_independent);
513 switch (idelta) {
514 case 0:
515 return delta * rho_r
516 * (2 * Delta_1 * Delta_2 * delta * rho_r * db_dxi * db_dxj
517 + (2 * Delta_1 * Delta_2 * bm * delta * rho_r + (Delta_1 + Delta_2) * (1 + cm * delta * rho_r)) * d2b_dxidxj);
518 case 1:
519 return rho_r
520 * (4 * Delta_1 * Delta_2 * delta * rho_r * db_dxi * db_dxj
521 + (4 * Delta_1 * Delta_2 * bm * delta * rho_r + (Delta_1 + Delta_2) * (1 + 2 * cm * delta * rho_r)) * d2b_dxidxj);
522 case 2:
523 return 2 * pow(rho_r, 2)
524 * (2 * Delta_1 * Delta_2 * db_dxi * db_dxj + (2 * Delta_1 * Delta_2 * bm + (Delta_1 + Delta_2) * cm) * d2b_dxidxj);
525 case 3:
526 return 0;
527 case 4:
528 return 0;
529 default:
530 throw -1;
531 }
532}
533double AbstractCubic::d3_PI_12_dxidxjdxk(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, std::size_t j, std::size_t k,
534 bool xN_independent) {
535 double bm = bm_term(x);
536 double cm = cm_term();
537 double db_dxi = d_bm_term_dxi(x, i, xN_independent), db_dxj = d_bm_term_dxi(x, j, xN_independent), db_dxk = d_bm_term_dxi(x, k, xN_independent),
538 d2b_dxidxj = d2_bm_term_dxidxj(x, i, j, xN_independent), d2b_dxidxk = d2_bm_term_dxidxj(x, i, k, xN_independent),
539 d2b_dxjdxk = d2_bm_term_dxidxj(x, j, k, xN_independent), d3b_dxidxjdxk = d3_bm_term_dxidxjdxk(x, i, j, k, xN_independent);
540 switch (idelta) {
541 case 0:
542 return delta * rho_r
543 * ((2 * Delta_1 * Delta_2 * bm * delta * rho_r + (Delta_1 + Delta_2) * (1 + cm * delta * rho_r)) * d3b_dxidxjdxk
544 + 2 * Delta_1 * Delta_2 * delta * rho_r * (db_dxi * d2b_dxjdxk + db_dxj * d2b_dxidxk + db_dxk * d2b_dxidxj));
545 case 1:
546 return rho_r
547 * ((4. * Delta_1 * Delta_2 * bm * delta * rho_r + (Delta_1 + Delta_2) * (1 + 2 * cm * delta * rho_r)) * d3b_dxidxjdxk
548 + 4 * Delta_1 * Delta_2 * delta * rho_r * (db_dxi * d2b_dxjdxk + db_dxj * d2b_dxidxk + db_dxk * d2b_dxidxj));
549 default:
550 throw -1;
551 }
552}
553double AbstractCubic::psi_plus(double delta, const std::vector<double>& x, std::size_t idelta) {
554 switch (idelta) {
555 case 0:
556 return A_term(delta, x) * c_term(x) / (Delta_1 - Delta_2);
557 case 1:
558 return rho_r / PI_12(delta, x, 0);
559 case 2: {
560 const double p0 = PI_12(delta, x, 0), p1 = PI_12(delta, x, 1);
561 return -rho_r / (p0 * p0) * p1;
562 }
563 case 3: {
564 const double p0 = PI_12(delta, x, 0), p1 = PI_12(delta, x, 1), p2 = PI_12(delta, x, 2);
565 return rho_r * (-p0 * p2 + 2 * p1 * p1) / (p0 * p0 * p0);
566 }
567 case 4: {
568 // Term -PI_12(delta,x,0)*PI_12(delta,x,3) in the numerator is zero (and removed) since PI_12(delta,x,3) = 0
569 const double p0 = PI_12(delta, x, 0), p1 = PI_12(delta, x, 1), p2 = PI_12(delta, x, 2);
570 return rho_r * (6 * p0 * p1 * p2 - 6 * p1 * p1 * p1) / (p0 * p0 * p0 * p0);
571 }
572 default:
573 throw -1;
574 }
575}
576double AbstractCubic::d_psi_plus_dxi(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, bool xN_independent) {
577 double bracket = 0;
578 if (idelta == 0) {
579 return (A_term(delta, x) * d_c_term_dxi(x, i, xN_independent) + c_term(x) * d_A_term_dxi(delta, x, i, xN_independent)) / (Delta_1 - Delta_2);
580 }
581 // All the terms with at least one delta derivative are multiplied by a common term of -rhor/PI12^2
582 // So we just evaluate the bracketed term and then multiply by the common factor in the front
583 switch (idelta) {
584 case 1:
585 bracket = d_PI_12_dxi(delta, x, 0, i, xN_independent);
586 break;
587 case 2:
588 bracket = (d_PI_12_dxi(delta, x, 1, i, xN_independent)
589 + 2 / rho_r * PI_12(delta, x, 0) * PI_12(delta, x, 1) * d_psi_plus_dxi(delta, x, 1, i, xN_independent));
590 break;
591 case 3: {
592 bracket =
593 (d_PI_12_dxi(delta, x, 2, i, xN_independent)
594 + 2 / rho_r * (pow(PI_12(delta, x, 1), 2) + PI_12(delta, x, 0) * PI_12(delta, x, 2)) * d_psi_plus_dxi(delta, x, 1, i, xN_independent)
595 + 4 / rho_r * PI_12(delta, x, 0) * PI_12(delta, x, 1) * d_psi_plus_dxi(delta, x, 2, i, xN_independent));
596 break;
597 }
598 case 4:
599 // d_PI_12_dxi(delta, x, 3, i, xN_independent) = 0, and PI_12(delta,x,0)*PI_12(delta,x,3) = 0, so removed from sum
600 bracket =
601 (6 / rho_r * PI_12(delta, x, 1) * PI_12(delta, x, 2) * d_psi_plus_dxi(delta, x, 1, i, xN_independent)
602 + 6 / rho_r * (pow(PI_12(delta, x, 1), 2) + PI_12(delta, x, 0) * PI_12(delta, x, 2)) * d_psi_plus_dxi(delta, x, 2, i, xN_independent)
603 + 6 / rho_r * PI_12(delta, x, 0) * PI_12(delta, x, 1) * d_psi_plus_dxi(delta, x, 3, i, xN_independent));
604 break;
605 default:
606 throw -1;
607 }
608 return -rho_r / pow(PI_12(delta, x, 0), 2) * bracket;
609}
610double AbstractCubic::d2_psi_plus_dxidxj(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, std::size_t j,
611 bool xN_independent) {
612 double bracket = 0;
613 double PI12 = PI_12(delta, x, 0);
614 if (idelta == 0) {
615 return (A_term(delta, x) * d2_c_term_dxidxj(x, i, j, xN_independent) + c_term(x) * d2_A_term_dxidxj(delta, x, i, j, xN_independent)
616 + d_A_term_dxi(delta, x, i, xN_independent) * d_c_term_dxi(x, j, xN_independent)
617 + d_A_term_dxi(delta, x, j, xN_independent) * d_c_term_dxi(x, i, xN_independent))
618 / (Delta_1 - Delta_2);
619 }
620 // All the terms with at least one delta derivative have a common factor of -1/PI_12^2 out front
621 // so we just calculate the bracketed term and then multiply later on
622 switch (idelta) {
623 case 1:
624 bracket = (rho_r * d2_PI_12_dxidxj(delta, x, 0, i, j, xN_independent)
625 + 2 * PI12 * d_PI_12_dxi(delta, x, 0, j, xN_independent) * d_psi_plus_dxi(delta, x, 1, i, xN_independent));
626 break;
627 case 2:
628 bracket = (rho_r * d2_PI_12_dxidxj(delta, x, 1, i, j, xN_independent)
629 + 2 * (PI12 * d_PI_12_dxi(delta, x, 1, j, xN_independent) + PI_12(delta, x, 1) * d_PI_12_dxi(delta, x, 0, j, xN_independent))
630 * d_psi_plus_dxi(delta, x, 1, i, xN_independent)
631 + 2 * PI12 * PI_12(delta, x, 1) * d2_psi_plus_dxidxj(delta, x, 1, i, j, xN_independent)
632 + 2 * PI12 * d_PI_12_dxi(delta, x, 0, j, xN_independent) * d_psi_plus_dxi(delta, x, 2, i, xN_independent));
633 break;
634 case 3: {
635 bracket =
636 (rho_r * d2_PI_12_dxidxj(delta, x, 2, i, j, xN_independent)
637 + 2 * (PI12 * PI_12(delta, x, 2) + pow(PI_12(delta, x, 1), 2)) * d2_psi_plus_dxidxj(delta, x, 1, i, j, xN_independent)
638 + 4 * (PI12 * d_PI_12_dxi(delta, x, 1, j, xN_independent) + PI_12(delta, x, 1) * d_PI_12_dxi(delta, x, 0, j, xN_independent))
639 * d_psi_plus_dxi(delta, x, 2, i, xN_independent)
640 + 2
641 * (PI12 * d_PI_12_dxi(delta, x, 2, j, xN_independent) + 2 * PI_12(delta, x, 1) * d_PI_12_dxi(delta, x, 1, j, xN_independent)
642 + d_PI_12_dxi(delta, x, 0, j, xN_independent) * PI_12(delta, x, 2))
643 * d_psi_plus_dxi(delta, x, 1, i, xN_independent)
644 + 4 * PI12 * PI_12(delta, x, 1) * d2_psi_plus_dxidxj(delta, x, 2, i, j, xN_independent)
645 + 2 * PI12 * d_PI_12_dxi(delta, x, 0, j, xN_independent) * d_psi_plus_dxi(delta, x, 3, i, xN_independent));
646 break;
647 }
648 case 4:
649 // rho_r*d2_PI_12_dxidxj(delta, x, 3, i, j, xN_independent) = 0
650 // PI_12(delta, x, 3) = 0
651 // PI12*d_PI_12_dxi(delta, x, 3, j, xN_independent) = 0
652 // d_PI_12_dxi(delta, x, 0, j, xN_independent)*PI_12(delta, x, 3) = 0
653 bracket =
654 (+6 * (PI12 * PI_12(delta, x, 2) + pow(PI_12(delta, x, 1), 2)) * d2_psi_plus_dxidxj(delta, x, 2, i, j, xN_independent)
655 + 6 * PI_12(delta, x, 1) * PI_12(delta, x, 2) * d2_psi_plus_dxidxj(delta, x, 1, i, j, xN_independent)
656 + 6 * (PI12 * d_PI_12_dxi(delta, x, 1, j, xN_independent) + PI_12(delta, x, 1) * d_PI_12_dxi(delta, x, 0, j, xN_independent))
657 * d_psi_plus_dxi(delta, x, 3, i, xN_independent)
658 + 6
659 * (PI12 * d_PI_12_dxi(delta, x, 2, j, xN_independent) + 2 * PI_12(delta, x, 1) * d_PI_12_dxi(delta, x, 1, j, xN_independent)
660 + d_PI_12_dxi(delta, x, 0, j, xN_independent) * PI_12(delta, x, 2))
661 * d_psi_plus_dxi(delta, x, 2, i, xN_independent)
662 + 6
663 * (PI_12(delta, x, 1) * d_PI_12_dxi(delta, x, 2, j, xN_independent)
664 + PI_12(delta, x, 2) * d_PI_12_dxi(delta, x, 1, j, xN_independent))
665 * d_psi_plus_dxi(delta, x, 1, i, xN_independent)
666 + 6 * PI12 * PI_12(delta, x, 1) * d2_psi_plus_dxidxj(delta, x, 3, i, j, xN_independent)
667 + 2 * PI12 * d_PI_12_dxi(delta, x, 0, j, xN_independent) * d_psi_plus_dxi(delta, x, 4, i, xN_independent));
668 break;
669 default:
670 throw -1;
671 }
672 return -1 / pow(PI12, 2) * bracket;
673}
674double AbstractCubic::d3_psi_plus_dxidxjdxk(double delta, const std::vector<double>& x, std::size_t idelta, std::size_t i, std::size_t j,
675 std::size_t k, bool xN_independent) {
676 double PI12 = PI_12(delta, x, 0);
677 switch (idelta) {
678 case 0:
679 return (A_term(delta, x) * d3_c_term_dxidxjdxk(x, i, j, k, xN_independent)
680 + c_term(x) * d3_A_term_dxidxjdxk(delta, x, i, j, k, xN_independent)
681 + d_A_term_dxi(delta, x, i, xN_independent) * d2_c_term_dxidxj(x, j, k, xN_independent)
682 + d_A_term_dxi(delta, x, j, xN_independent) * d2_c_term_dxidxj(x, i, k, xN_independent)
683 + d_A_term_dxi(delta, x, k, xN_independent) * d2_c_term_dxidxj(x, i, j, xN_independent)
684 + d_c_term_dxi(x, i, xN_independent) * d2_A_term_dxidxj(delta, x, j, k, xN_independent)
685 + d_c_term_dxi(x, j, xN_independent) * d2_A_term_dxidxj(delta, x, i, k, xN_independent)
686 + d_c_term_dxi(x, k, xN_independent) * d2_A_term_dxidxj(delta, x, i, j, xN_independent))
687 / (Delta_1 - Delta_2);
688 case 1:
689 return -1 / pow(PI12, 2)
690 * (rho_r * d3_PI_12_dxidxjdxk(delta, x, 0, i, j, k, xN_independent)
691 + 2
692 * (PI12 * d2_PI_12_dxidxj(delta, x, 0, j, k, xN_independent)
693 + d_PI_12_dxi(delta, x, 0, j, xN_independent) * d_PI_12_dxi(delta, x, 0, k, xN_independent))
694 * d_psi_plus_dxi(delta, x, 1, i, xN_independent)
695 + 2 * PI12 * d_PI_12_dxi(delta, x, 0, j, xN_independent) * d2_psi_plus_dxidxj(delta, x, 1, i, k, xN_independent)
696 + 2 * PI12 * d_PI_12_dxi(delta, x, 0, k, xN_independent) * d2_psi_plus_dxidxj(delta, x, 1, i, j, xN_independent));
697 default:
698 throw -1;
699 }
700}
701
702double AbstractCubic::tau_times_a(double tau, const std::vector<double>& x, std::size_t itau) {
703 if (itau == 0) {
704 return tau * am_term(tau, x, 0);
705 } else {
706 return tau * am_term(tau, x, itau) + itau * am_term(tau, x, itau - 1);
707 }
708}
709double AbstractCubic::d_tau_times_a_dxi(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, bool xN_independent) {
710 if (itau == 0) {
711 return tau * d_am_term_dxi(tau, x, 0, i, xN_independent);
712 } else {
713 return tau * d_am_term_dxi(tau, x, itau, i, xN_independent) + itau * d_am_term_dxi(tau, x, itau - 1, i, xN_independent);
714 }
715}
716double AbstractCubic::d2_tau_times_a_dxidxj(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, std::size_t j,
717 bool xN_independent) {
718 if (itau == 0) {
719 return tau * d2_am_term_dxidxj(tau, x, 0, i, j, xN_independent);
720 } else {
721 return tau * d2_am_term_dxidxj(tau, x, itau, i, j, xN_independent) + itau * d2_am_term_dxidxj(tau, x, itau - 1, i, j, xN_independent);
722 }
723}
724double AbstractCubic::d3_tau_times_a_dxidxjdxk(double tau, const std::vector<double>& x, std::size_t itau, std::size_t i, std::size_t j,
725 std::size_t k, bool xN_independent) {
726 if (itau == 0) {
727 return tau * d3_am_term_dxidxjdxk(tau, x, 0, i, j, k, xN_independent);
728 } else {
729 return tau * d3_am_term_dxidxjdxk(tau, x, itau, i, j, k, xN_independent)
730 + itau * d3_am_term_dxidxjdxk(tau, x, itau - 1, i, j, k, xN_independent);
731 }
732}
733double AbstractCubic::alphar(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta) {
734 return psi_minus(delta, x, itau, idelta) - tau_times_a(tau, x, itau) / (R_u * T_r) * psi_plus(delta, x, idelta);
735}
736double AbstractCubic::d_alphar_dxi(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
737 bool xN_independent) {
738 return (d_psi_minus_dxi(delta, x, itau, idelta, i, xN_independent)
739 - 1 / (R_u * T_r)
740 * (d_tau_times_a_dxi(tau, x, itau, i, xN_independent) * psi_plus(delta, x, idelta)
741 + tau_times_a(tau, x, itau) * d_psi_plus_dxi(delta, x, idelta, i, xN_independent)));
742}
743double AbstractCubic::d2_alphar_dxidxj(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
744 std::size_t j, bool xN_independent) {
745 return (d2_psi_minus_dxidxj(delta, x, itau, idelta, i, j, xN_independent)
746 - 1 / (R_u * T_r)
747 * (d2_tau_times_a_dxidxj(tau, x, itau, i, j, xN_independent) * psi_plus(delta, x, idelta)
748 + d_tau_times_a_dxi(tau, x, itau, i, xN_independent) * d_psi_plus_dxi(delta, x, idelta, j, xN_independent)
749 + d_tau_times_a_dxi(tau, x, itau, j, xN_independent) * d_psi_plus_dxi(delta, x, idelta, i, xN_independent)
750 + tau_times_a(tau, x, itau) * d2_psi_plus_dxidxj(delta, x, idelta, i, j, xN_independent)));
751}
752double AbstractCubic::d3_alphar_dxidxjdxk(double tau, double delta, const std::vector<double>& x, std::size_t itau, std::size_t idelta, std::size_t i,
753 std::size_t j, std::size_t k, bool xN_independent) {
754 return (d3_psi_minus_dxidxjdxk(delta, x, itau, idelta, i, j, k, xN_independent)
755 - 1 / (R_u * T_r)
756 * (d2_tau_times_a_dxidxj(tau, x, itau, i, j, xN_independent) * d_psi_plus_dxi(delta, x, idelta, k, xN_independent)
757 + d3_tau_times_a_dxidxjdxk(tau, x, itau, i, j, k, xN_independent) * psi_plus(delta, x, idelta)
758
759 + d_tau_times_a_dxi(tau, x, itau, i, xN_independent) * d2_psi_plus_dxidxj(delta, x, idelta, j, k, xN_independent)
760 + d2_tau_times_a_dxidxj(tau, x, itau, i, k, xN_independent) * d_psi_plus_dxi(delta, x, idelta, j, xN_independent)
761
762 + d_tau_times_a_dxi(tau, x, itau, j, xN_independent) * d2_psi_plus_dxidxj(delta, x, idelta, i, k, xN_independent)
763 + d2_tau_times_a_dxidxj(tau, x, itau, j, k, xN_independent) * d_psi_plus_dxi(delta, x, idelta, i, xN_independent)
764
765 + tau_times_a(tau, x, itau) * d3_psi_plus_dxidxjdxk(delta, x, idelta, i, j, k, xN_independent)
766 + d_tau_times_a_dxi(tau, x, itau, k, xN_independent) * d2_psi_plus_dxidxj(delta, x, idelta, i, j, xN_independent)));
767}
768
769double SRK::a0_ii(std::size_t i) {
770 // Exact value: 1/(9*(2^(1/3)-1)); see Bell and Deiters, IECR, 2021
771 double a = 0.42748023335403414043900347952220 * R_u * R_u * Tc[i] * Tc[i] / pc[i];
772 return a;
773}
774double SRK::b0_ii(std::size_t i) {
775 // Exact value: (2^(1/3)-1)/3; see Bell and Deiters, IECR, 2021
776 double b = 0.08664034999649577215890158147700 * R_u * Tc[i] / pc[i];
777 return b;
778}
779double SRK::m_ii(std::size_t i) {
780 // Values from Soave, 1972 (Equilibrium constants from a ..)
781 double omega = acentric[i];
782 double m = 0.480 + 1.574 * omega - 0.176 * omega * omega;
783 return m;
784}
785
786double PengRobinson::a0_ii(std::size_t i) {
787 // Exact value; see Bell and Deiters, IECR, 2021
788 double a = 0.45723552892138218938000849856422 * R_u * R_u * Tc[i] * Tc[i] / pc[i];
789 return a;
790}
791double PengRobinson::b0_ii(std::size_t i) {
792 // Exact value; see Bell and Deiters, IECR, 2021
793 double b = 0.07779607390388455972148597969400 * R_u * Tc[i] / pc[i];
794 return b;
795}
796double PengRobinson::m_ii(std::size_t i) {
797 double omega = acentric[i];
798 double m = 0.37464 + 1.54226 * omega - 0.26992 * omega * omega;
799 return m;
800}
Generated by doxygen 1.9.4