Loading [MathJax]/extensions/TeX/AMSsymbols.js
CoolProp 6.8.1dev
An open-source fluid property and humid air property database
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Helmholtz.cpp
Go to the documentation of this file.
1#include <numeric>
2#include "Helmholtz.h"
3
4#ifdef __ANDROID__
5# undef _A
6# undef _B
7# undef _C
8# undef _D
9#endif
10
11namespace CoolProp {
12
13CoolPropDbl kahanSum(const std::vector<CoolPropDbl>& x) {
14 CoolPropDbl sum = x[0], y, t;
15 CoolPropDbl c = 0.0; //A running compensation for lost low-order bits.
16 for (std::size_t i = 1; i < x.size(); ++i) {
17 y = x[i] - c; //So far, so good: c is zero.
18 t = sum + y; //Alas, sum is big, y small, so low-order digits of y are lost.
19 c = (t - sum) - y; //(t - sum) recovers the high-order part of y; subtracting y recovers -(low part of y)
20 sum = t; //Algebraically, c should always be zero. Beware eagerly optimising compilers!
21 }
22 return sum;
23}
24bool wayToSort(CoolPropDbl i, CoolPropDbl j) {
25 return std::abs(i) > std::abs(j);
26}
27
28// define function to be applied coefficient-wise
29double ramp(double x) {
30 if (x > 0)
31 return x;
32 else
33 return 0;
34}
35
36/*
37void ResidualHelmholtzGeneralizedExponential::allEigen(const CoolPropDbl &tau, const CoolPropDbl &delta, HelmholtzDerivatives &derivs) throw()
38{
39 double log_tau = log(tau), log_delta = log(delta),
40 one_over_delta = 1/delta, one_over_tau = 1/tau; // division is much slower than multiplication, so do one division here
41
42 Eigen::Map<Eigen::ArrayXd> nE(&(n[0]), elements.size());
43 Eigen::Map<Eigen::ArrayXd> dE(&(d[0]), elements.size());
44 Eigen::Map<Eigen::ArrayXd> tE(&(t[0]), elements.size());
45 Eigen::Map<Eigen::ArrayXd> cE(&(c[0]), elements.size());
46 Eigen::Map<Eigen::ArrayXi> l_intE(&(l_int[0]), elements.size());
47 Eigen::Map<Eigen::ArrayXd> l_doubleE(&(l_double[0]), elements.size());
48 Eigen::Map<Eigen::ArrayXd> eta1E(&(eta1[0]), elements.size());
49 Eigen::Map<Eigen::ArrayXd> eta2E(&(eta2[0]), elements.size());
50 Eigen::Map<Eigen::ArrayXd> epsilon1E(&(epsilon1[0]), elements.size());
51 Eigen::Map<Eigen::ArrayXd> epsilon2E(&(epsilon2[0]), elements.size());
52 Eigen::Map<Eigen::ArrayXd> beta1E(&(beta1[0]), elements.size());
53 Eigen::Map<Eigen::ArrayXd> beta2E(&(beta2[0]), elements.size());
54 Eigen::Map<Eigen::ArrayXd> gamma1E(&(gamma1[0]), elements.size());
55 Eigen::Map<Eigen::ArrayXd> gamma2E(&(gamma2[0]), elements.size());
56
57 // ****************************************
58 // The u part in exp(u) and its derivatives
59 // ****************************************
60
61 #if defined(EIGEN_VECTORIZE_SSE2)
62 //std::cout << "EIGEN_VECTORIZE_SSE2" << std::endl;
63 #endif
64
65 // Set the u part of exp(u) to zero
66 uE.fill(0);
67 du_ddeltaE.fill(0);
68 du_dtauE.fill(0);
69 d2u_ddelta2E.fill(0);
70 d2u_dtau2E.fill(0);
71 d3u_ddelta3E.fill(0);
72 d3u_dtau3E.fill(0);
73
74 if (delta_li_in_u){
75 Eigen::ArrayXd u_increment = -cE*(log_delta*l_doubleE).exp(); //pow(delta,L) -> exp(L*log(delta))
76 uE += u_increment;
77 du_ddeltaE += l_doubleE*u_increment*one_over_delta;
78 d2u_ddelta2E += (l_doubleE-1)*l_doubleE*u_increment*one_over_delta*one_over_delta;
79 d3u_ddelta3E += (l_doubleE-2)*(l_doubleE-1)*l_doubleE*u_increment*one_over_delta*one_over_delta*one_over_delta;
80 }
81
82// if (tau_mi_in_u){
83// CoolPropDbl omegai = el.omega, m_double = el.m_double;
84// if (std::abs(m_double) > 0){
85// CoolPropDbl u_increment = -omegai*pow(tau, m_double);
86// CoolPropDbl du_dtau_increment = m_double*u_increment*one_over_tau;
87// CoolPropDbl d2u_dtau2_increment = (m_double-1)*du_dtau_increment*one_over_tau;
88// CoolPropDbl d3u_dtau3_increment = (m_double-2)*d2u_dtau2_increment*one_over_tau;
89// u += u_increment;
90// du_dtau += du_dtau_increment;
91// d2u_dtau2 += d2u_dtau2_increment;
92// d3u_dtau3 += d3u_dtau3_increment;
93// }
94// }
95 if (eta1_in_u){
96 uE += -eta1E*(delta-epsilon1E);
97 du_ddeltaE += -eta1E;
98 }
99 if (eta2_in_u){
100 uE += -eta2E*POW2(delta-epsilon2E);
101 du_ddeltaE += -2*eta2E*(delta-epsilon2E);
102 d2u_ddelta2E += -2*eta2E;
103 }
104 if (beta1_in_u){
105 uE += -beta1E*(tau-gamma1E);
106 du_dtauE += -beta1E;
107 }
108 if (beta2_in_u){
109 uE += -beta2E*POW2(tau-gamma2E);
110 du_dtauE += -2*beta2E*(tau-gamma2E);
111 d2u_dtau2E += -2*beta2E;
112 }
113
114 Eigen::ArrayXd ndteuE = nE*exp(tE*log_tau + dE*log_delta + uE);
115 Eigen::ArrayXd B_deltaE = delta*du_ddeltaE + dE;
116 Eigen::ArrayXd B_tauE = tau*du_dtauE + tE;
117 Eigen::ArrayXd B_delta2E = POW2(delta)*(d2u_ddelta2E + du_ddeltaE.square()) + 2*dE*delta*du_ddeltaE + dE*(dE-1);
118 Eigen::ArrayXd B_tau2E = POW2(tau)*(d2u_dtau2E + du_dtauE.square()) + 2*tE*tau*du_dtauE + tE*(tE-1);
119 Eigen::ArrayXd B_delta3E = POW3(delta)*d3u_ddelta3E + 3*dE*POW2(delta)*d2u_ddelta2E+3*POW3(delta)*d2u_ddelta2E*du_ddeltaE+3*dE*POW2(delta*du_ddeltaE)+3*dE*(dE-1)*delta*du_ddeltaE+dE*(dE-1)*(dE-2)+POW3(delta*du_ddeltaE);
120 Eigen::ArrayXd B_tau3E = POW3(tau)*d3u_dtau3E + 3*tE*POW2(tau)*d2u_dtau2E+3*POW3(tau)*d2u_dtau2E*du_dtauE+3*tE*POW2(tau*du_dtauE)+3*tE*(tE-1)*tau*du_dtauE+tE*(tE-1)*(tE-2)+POW3(tau*du_dtauE);
121
122 derivs.alphar += ndteuE.sum();
123 derivs.dalphar_ddelta += (ndteuE*B_deltaE).sum()*one_over_delta;
124 derivs.dalphar_dtau += (ndteuE*B_tauE).sum()*one_over_tau;
125 derivs.d2alphar_ddelta2 += (ndteuE*B_delta2E).sum()*POW2(one_over_delta);
126 derivs.d2alphar_dtau2 += (ndteuE*B_tau2E).sum()*POW2(one_over_tau);
127 derivs.d2alphar_ddelta_dtau += (ndteuE*B_deltaE*B_tauE).sum()*one_over_delta*one_over_tau;
128
129 derivs.d3alphar_ddelta3 += (ndteuE*B_delta3E).sum()*POW3(one_over_delta);
130 derivs.d3alphar_dtau3 += (ndteuE*B_tau3E).sum()*POW3(one_over_tau);
131 derivs.d3alphar_ddelta2_dtau += (ndteuE*B_delta2E*B_tauE).sum()*POW2(one_over_delta)*one_over_tau;
132 derivs.d3alphar_ddelta_dtau2 += (ndteuE*B_deltaE*B_tau2E).sum()*one_over_delta*POW2(one_over_tau);
133
134 return;
135};
136*/
137void ResidualHelmholtzGeneralizedExponential::all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) {
138 CoolPropDbl log_tau = log(tau), log_delta = log(delta), ndteu, one_over_delta = 1 / delta,
139 one_over_tau = 1 / tau; // division is much slower than multiplication, so do one division here
140
141 // Maybe split the construction of u and other parts into two separate loops?
142 // If both loops can get vectorized, could be worth it.
143 const std::size_t N = elements.size();
144 for (std::size_t i = 0; i < N; ++i) {
145 ResidualHelmholtzGeneralizedExponentialElement& el = elements[i];
146 CoolPropDbl ni = el.n, di = el.d, ti = el.t;
147
148 // Set the u part of exp(u) to zero
149 CoolPropDbl u = 0;
150 CoolPropDbl du_ddelta = 0;
151 CoolPropDbl du_dtau = 0;
152 CoolPropDbl d2u_ddelta2 = 0;
153 CoolPropDbl d2u_dtau2 = 0;
154 CoolPropDbl d3u_ddelta3 = 0;
155 CoolPropDbl d3u_dtau3 = 0;
156 CoolPropDbl d4u_ddelta4 = 0;
157 CoolPropDbl d4u_dtau4 = 0;
158
159 if (delta_li_in_u) {
160 CoolPropDbl ci = el.c, l_double = el.l_double;
161 if (ValidNumber(l_double) && l_double > 0 && std::abs(ci) > DBL_EPSILON) {
162 const CoolPropDbl u_increment = (el.l_is_int) ? -ci * powInt(delta, el.l_int) : -ci * pow(delta, l_double);
163 const CoolPropDbl du_ddelta_increment = l_double * u_increment * one_over_delta;
164 const CoolPropDbl d2u_ddelta2_increment = (l_double - 1) * du_ddelta_increment * one_over_delta;
165 const CoolPropDbl d3u_ddelta3_increment = (l_double - 2) * d2u_ddelta2_increment * one_over_delta;
166 const CoolPropDbl d4u_ddelta4_increment = (l_double - 3) * d3u_ddelta3_increment * one_over_delta;
167 u += u_increment;
168 du_ddelta += du_ddelta_increment;
169 d2u_ddelta2 += d2u_ddelta2_increment;
170 d3u_ddelta3 += d3u_ddelta3_increment;
171 d4u_ddelta4 += d4u_ddelta4_increment;
172 }
173 }
174 if (tau_mi_in_u) {
175 CoolPropDbl omegai = el.omega, m_double = el.m_double;
176 if (std::abs(m_double) > 0) {
177 const CoolPropDbl u_increment = -omegai * pow(tau, m_double);
178 const CoolPropDbl du_dtau_increment = m_double * u_increment * one_over_tau;
179 const CoolPropDbl d2u_dtau2_increment = (m_double - 1) * du_dtau_increment * one_over_tau;
180 const CoolPropDbl d3u_dtau3_increment = (m_double - 2) * d2u_dtau2_increment * one_over_tau;
181 const CoolPropDbl d4u_dtau4_increment = (m_double - 3) * d3u_dtau3_increment * one_over_tau;
182 u += u_increment;
183 du_dtau += du_dtau_increment;
184 d2u_dtau2 += d2u_dtau2_increment;
185 d3u_dtau3 += d3u_dtau3_increment;
186 d4u_dtau4 += d4u_dtau4_increment;
187 }
188 }
189 if (eta1_in_u) {
190 CoolPropDbl eta1 = el.eta1, epsilon1 = el.epsilon1;
191 if (ValidNumber(eta1)) {
192 u += -eta1 * (delta - epsilon1);
193 du_ddelta += -eta1;
194 }
195 }
196 if (eta2_in_u) {
197 CoolPropDbl eta2 = el.eta2, epsilon2 = el.epsilon2;
198 if (ValidNumber(eta2)) {
199 u += -eta2 * POW2(delta - epsilon2);
200 du_ddelta += -2 * eta2 * (delta - epsilon2);
201 d2u_ddelta2 += -2 * eta2;
202 }
203 }
204 if (beta1_in_u) {
205 CoolPropDbl beta1 = el.beta1, gamma1 = el.gamma1;
206 if (ValidNumber(beta1)) {
207 u += -beta1 * (tau - gamma1);
208 du_dtau += -beta1;
209 }
210 }
211 if (beta2_in_u) {
212 CoolPropDbl beta2 = el.beta2, gamma2 = el.gamma2;
213 if (ValidNumber(beta2)) {
214 u += -beta2 * POW2(tau - gamma2);
215 du_dtau += -2 * beta2 * (tau - gamma2);
216 d2u_dtau2 += -2 * beta2;
217 }
218 }
219
220 ndteu = ni * exp(ti * log_tau + di * log_delta + u);
221
222 const CoolPropDbl dB_delta_ddelta = delta * d2u_ddelta2 + du_ddelta;
223 const CoolPropDbl d2B_delta_ddelta2 = delta * d3u_ddelta3 + 2 * d2u_ddelta2;
224 const CoolPropDbl d3B_delta_ddelta3 = delta * d4u_ddelta4 + 3 * d3u_ddelta3;
225
226 const CoolPropDbl B_delta = (delta * du_ddelta + di);
227 const CoolPropDbl B_delta2 = delta * dB_delta_ddelta + (B_delta - 1) * B_delta;
228 const CoolPropDbl dB_delta2_ddelta = delta * d2B_delta_ddelta2 + 2 * B_delta * dB_delta_ddelta;
229 const CoolPropDbl B_delta3 = delta * dB_delta2_ddelta + (B_delta - 2) * B_delta2;
230 const CoolPropDbl dB_delta3_ddelta = delta * delta * d3B_delta_ddelta3 + 3 * delta * B_delta * d2B_delta_ddelta2
231 + 3 * delta * POW2(dB_delta_ddelta) + 3 * B_delta * (B_delta - 1) * dB_delta_ddelta;
232 const CoolPropDbl B_delta4 = delta * dB_delta3_ddelta + (B_delta - 3) * B_delta3;
233
234 const CoolPropDbl dB_tau_dtau = tau * d2u_dtau2 + du_dtau;
235 const CoolPropDbl d2B_tau_dtau2 = tau * d3u_dtau3 + 2 * d2u_dtau2;
236 const CoolPropDbl d3B_tau_dtau3 = tau * d4u_dtau4 + 3 * d3u_dtau3;
237
238 const CoolPropDbl B_tau = (tau * du_dtau + ti);
239 const CoolPropDbl B_tau2 = tau * dB_tau_dtau + (B_tau - 1) * B_tau;
240 const CoolPropDbl dB_tau2_dtau = tau * d2B_tau_dtau2 + 2 * B_tau * dB_tau_dtau;
241 const CoolPropDbl B_tau3 = tau * dB_tau2_dtau + (B_tau - 2) * B_tau2;
242 const CoolPropDbl dB_tau3_dtau =
243 tau * tau * d3B_tau_dtau3 + 3 * tau * B_tau * d2B_tau_dtau2 + 3 * tau * POW2(dB_tau_dtau) + 3 * B_tau * (B_tau - 1) * dB_tau_dtau;
244 const CoolPropDbl B_tau4 = tau * dB_tau3_dtau + (B_tau - 3) * B_tau3;
245
246 derivs.alphar += ndteu;
247
248 derivs.dalphar_ddelta += ndteu * B_delta;
249 derivs.dalphar_dtau += ndteu * B_tau;
250
251 derivs.d2alphar_ddelta2 += ndteu * B_delta2;
252 derivs.d2alphar_ddelta_dtau += ndteu * B_delta * B_tau;
253 derivs.d2alphar_dtau2 += ndteu * B_tau2;
254
255 derivs.d3alphar_ddelta3 += ndteu * B_delta3;
256 derivs.d3alphar_ddelta2_dtau += ndteu * B_delta2 * B_tau;
257 derivs.d3alphar_ddelta_dtau2 += ndteu * B_delta * B_tau2;
258 derivs.d3alphar_dtau3 += ndteu * B_tau3;
259
260 derivs.d4alphar_ddelta4 += ndteu * B_delta4;
261 derivs.d4alphar_ddelta3_dtau += ndteu * B_delta3 * B_tau;
262 derivs.d4alphar_ddelta2_dtau2 += ndteu * B_delta2 * B_tau2;
263 derivs.d4alphar_ddelta_dtau3 += ndteu * B_delta * B_tau3;
264 derivs.d4alphar_dtau4 += ndteu * B_tau4;
265 }
266 derivs.dalphar_ddelta *= one_over_delta;
267 derivs.dalphar_dtau *= one_over_tau;
268 derivs.d2alphar_ddelta2 *= POW2(one_over_delta);
269 derivs.d2alphar_dtau2 *= POW2(one_over_tau);
270 derivs.d2alphar_ddelta_dtau *= one_over_delta * one_over_tau;
271
272 derivs.d3alphar_ddelta3 *= POW3(one_over_delta);
273 derivs.d3alphar_dtau3 *= POW3(one_over_tau);
274 derivs.d3alphar_ddelta2_dtau *= POW2(one_over_delta) * one_over_tau;
275 derivs.d3alphar_ddelta_dtau2 *= one_over_delta * POW2(one_over_tau);
276
277 derivs.d4alphar_ddelta4 *= POW4(one_over_delta);
278 derivs.d4alphar_dtau4 *= POW4(one_over_tau);
279 derivs.d4alphar_ddelta3_dtau *= POW3(one_over_delta) * one_over_tau;
280 derivs.d4alphar_ddelta2_dtau2 *= POW2(one_over_delta) * POW2(one_over_tau);
281 derivs.d4alphar_ddelta_dtau3 *= one_over_delta * POW3(one_over_tau);
282
283 return;
284};
285
286#if ENABLE_CATCH
287mcx::MultiComplex<double> ResidualHelmholtzGeneralizedExponential::one_mcx(const mcx::MultiComplex<double>& tau,
288 const mcx::MultiComplex<double>& delta) const {
289 //throw CoolProp::NotImplementedError("Nope");
290 mcx::MultiComplex<double> sum00 = 0.0*tau*delta;
291 auto ln_tau = log(tau);
292 auto ln_delta = log(delta);
293 const std::size_t N = elements.size();
294 for (std::size_t i = 0; i < N; ++i) {
295 const ResidualHelmholtzGeneralizedExponentialElement& el = elements[i];
296
297 mcx::MultiComplex<double> u = 0.0 * tau * delta;
298 if (delta_li_in_u) {
299 CoolPropDbl ci = el.c, l_double = el.l_double;
300 if (ValidNumber(l_double) && l_double > 0 && std::abs(ci) > DBL_EPSILON) {
301 const auto u_increment = -ci * pow(delta, l_double);
302 u += u_increment;
303 }
304 }
305 if (tau_mi_in_u) {
306 CoolPropDbl omegai = el.omega, m_double = el.m_double;
307 if (std::abs(m_double) > 0) {
308 const auto u_increment = -omegai * pow(tau, m_double);
309 u += u_increment;
310 }
311 }
312 if (eta1_in_u) {
313 CoolPropDbl eta1 = el.eta1, epsilon1 = el.epsilon1;
314 if (ValidNumber(eta1)) {
315 u += -eta1 * (delta - epsilon1);
316 }
317 }
318 if (eta2_in_u) {
319 CoolPropDbl eta2 = el.eta2, epsilon2 = el.epsilon2;
320 if (ValidNumber(eta2)) {
321 u += -eta2 * POW2(delta - epsilon2);
322 }
323 }
324 if (beta1_in_u) {
325 CoolPropDbl beta1 = el.beta1, gamma1 = el.gamma1;
326 if (ValidNumber(beta1)) {
327 u += -beta1 * (tau - gamma1);
328 }
329 }
330 if (beta2_in_u) {
331 CoolPropDbl beta2 = el.beta2, gamma2 = el.gamma2;
332 if (ValidNumber(beta2)) {
333 u += -beta2 * POW2(tau - gamma2);
334 }
335 }
336 sum00 += el.n * exp(el.t * ln_tau + el.d * ln_delta + u);
337 }
338 return sum00;
339}
340#endif
341
342void ResidualHelmholtzGeneralizedExponential::to_json(rapidjson::Value& el, rapidjson::Document& doc) {
343 el.AddMember("type", "GeneralizedExponential", doc.GetAllocator());
344 cpjson::set_double_array("n", n, el, doc);
345 cpjson::set_double_array("t", t, el, doc);
346 cpjson::set_double_array("d", d, el, doc);
347 cpjson::set_double_array("eta1", eta1, el, doc);
348 cpjson::set_double_array("eta2", eta2, el, doc);
349 cpjson::set_double_array("beta1", beta1, el, doc);
350 cpjson::set_double_array("beta2", beta2, el, doc);
351 cpjson::set_double_array("gamma1", gamma1, el, doc);
352 cpjson::set_double_array("gamma2", gamma2, el, doc);
353 cpjson::set_double_array("epsilon1", epsilon1, el, doc);
354 cpjson::set_double_array("epsilon2", epsilon2, el, doc);
355 cpjson::set_double_array("l_double", l_double, el, doc);
356 cpjson::set_int_array("l_int", l_int, el, doc);
357}
358
359void ResidualHelmholtzNonAnalytic::to_json(rapidjson::Value& el, rapidjson::Document& doc) {
360 el.AddMember("type", "ResidualHelmholtzNonAnalytic", doc.GetAllocator());
361
362 rapidjson::Value _n(rapidjson::kArrayType), _a(rapidjson::kArrayType), _b(rapidjson::kArrayType), _beta(rapidjson::kArrayType),
363 _A(rapidjson::kArrayType), _B(rapidjson::kArrayType), _C(rapidjson::kArrayType), _D(rapidjson::kArrayType);
364 for (unsigned int i = 0; i <= N; ++i) {
365 ResidualHelmholtzNonAnalyticElement& elem = elements[i];
366 _n.PushBack((double)elem.n, doc.GetAllocator());
367 _a.PushBack((double)elem.a, doc.GetAllocator());
368 _b.PushBack((double)elem.b, doc.GetAllocator());
369 _beta.PushBack((double)elem.beta, doc.GetAllocator());
370 _A.PushBack((double)elem.A, doc.GetAllocator());
371 _B.PushBack((double)elem.B, doc.GetAllocator());
372 _C.PushBack((double)elem.C, doc.GetAllocator());
373 _D.PushBack((double)elem.D, doc.GetAllocator());
374 }
375 el.AddMember("n", _n, doc.GetAllocator());
376 el.AddMember("a", _a, doc.GetAllocator());
377 el.AddMember("b", _b, doc.GetAllocator());
378 el.AddMember("beta", _beta, doc.GetAllocator());
379 el.AddMember("A", _A, doc.GetAllocator());
380 el.AddMember("B", _B, doc.GetAllocator());
381 el.AddMember("C", _C, doc.GetAllocator());
382 el.AddMember("D", _D, doc.GetAllocator());
383}
384
385void ResidualHelmholtzNonAnalytic::all(const CoolPropDbl& tau_in, const CoolPropDbl& delta_in, HelmholtzDerivatives& derivs) {
386 if (N == 0) {
387 return;
388 }
389
390 // Here we want to hack this function just a tiny bit to avoid evaluation AT the critical point
391 // If we are VERY close to the critical point, just offset us a tiny bit away
392 CoolPropDbl tau = tau_in, delta = delta_in;
393 if (std::abs(tau_in - 1) < 10 * DBL_EPSILON) {
394 tau = 1.0 + 10 * DBL_EPSILON;
395 }
396 if (std::abs(delta_in - 1) < 10 * DBL_EPSILON) {
397 delta = 1.0 + 10 * DBL_EPSILON;
398 }
399
400 for (unsigned int i = 0; i < N; ++i) {
401 const ResidualHelmholtzNonAnalyticElement& el = elements[i];
402 const CoolPropDbl ni = el.n, ai = el.a, bi = el.b, betai = el.beta;
403 const CoolPropDbl Ai = el.A, Bi = el.B, Ci = el.C, Di = el.D;
404
405 // Derivatives of theta (all others are zero) (OK - checked)
406 // Do not factor because then when delta = 1 you are dividing by 0
407 const CoolPropDbl theta = (1.0 - tau) + Ai * pow(POW2(delta - 1.0), 1.0 / (2.0 * betai));
408 const CoolPropDbl dtheta_dTau = -1;
409 const CoolPropDbl dtheta_dDelta = Ai / (betai)*pow(POW2(delta - 1), 1 / (2 * betai) - 1) * (delta - 1);
410
411 const CoolPropDbl d2theta_dDelta2 = Ai / betai * (1 / betai - 1) * pow(POW2(delta - 1), 1 / (2 * betai) - 1);
412 const CoolPropDbl d3theta_dDelta3 = Ai / betai * (2 - 3 / betai + 1 / POW2(betai)) * pow(POW2(delta - 1), 1 / (2 * betai)) / POW3(delta - 1);
413 const CoolPropDbl d4theta_dDelta4 =
414 Ai / betai * (-6 + 11 / betai - 6 / POW2(betai) + 1 / POW3(betai)) * pow(POW2(delta - 1), 1 / (2 * betai) - 2);
415
416 // Derivatives of PSI (OK - checked)
417 const CoolPropDbl PSI = exp(-Ci * POW2(delta - 1.0) - Di * POW2(tau - 1.0));
418 const CoolPropDbl dPSI_dDelta_over_PSI = -2.0 * Ci * (delta - 1.0);
419 const CoolPropDbl dPSI_dDelta = dPSI_dDelta_over_PSI * PSI;
420 const CoolPropDbl dPSI_dTau_over_PSI = -2.0 * Di * (tau - 1.0);
421 const CoolPropDbl dPSI_dTau = dPSI_dTau_over_PSI * PSI;
422 const CoolPropDbl d2PSI_dDelta2_over_PSI = (2.0 * Ci * POW2(delta - 1.0) - 1.0) * 2.0 * Ci;
423 const CoolPropDbl d2PSI_dDelta2 = d2PSI_dDelta2_over_PSI * PSI;
424 const CoolPropDbl d3PSI_dDelta3 = 2 * Ci * PSI * (-4 * Ci * Ci * POW3(delta - 1) + 6 * Ci * (delta - 1));
425 const CoolPropDbl d4PSI_dDelta4 = 4 * Ci * Ci * PSI * (4 * Ci * Ci * POW4(delta - 1) - 12 * Ci * POW2(delta - 1) + 3);
426 const CoolPropDbl d2PSI_dTau2 = (2.0 * Di * POW2(tau - 1.0) - 1.0) * 2.0 * Di * PSI;
427 const CoolPropDbl d3PSI_dTau3 = 2.0 * Di * PSI * (-4 * Di * Di * POW3(tau - 1) + 6 * Di * (tau - 1));
428 const CoolPropDbl d4PSI_dTau4 = 4 * Di * Di * PSI * (4 * Di * Di * POW4(tau - 1) - 12 * Di * POW2(tau - 1) + 3);
429 const CoolPropDbl d2PSI_dDelta_dTau = dPSI_dDelta * dPSI_dTau_over_PSI;
430 const CoolPropDbl d3PSI_dDelta2_dTau = d2PSI_dDelta2 * dPSI_dTau_over_PSI;
431 const CoolPropDbl d3PSI_dDelta_dTau2 = d2PSI_dTau2 * dPSI_dDelta_over_PSI;
432 const CoolPropDbl d4PSI_dDelta_dTau3 = d3PSI_dTau3 * dPSI_dDelta_over_PSI;
433 const CoolPropDbl d4PSI_dDelta2_dTau2 = d2PSI_dTau2 * d2PSI_dDelta2_over_PSI;
434 const CoolPropDbl d4PSI_dDelta3_dTau = d3PSI_dDelta3 * dPSI_dTau_over_PSI;
435
436 // Derivatives of DELTA (OK - Checked)
437 const CoolPropDbl DELTA = POW2(theta) + Bi * pow(POW2(delta - 1.0), ai);
438 const CoolPropDbl dDELTA_dTau = 2 * theta * dtheta_dTau;
439 const CoolPropDbl dDELTA_dDelta = 2 * theta * dtheta_dDelta + 2 * Bi * ai * pow(POW2(delta - 1.0), ai - 1.0) * (delta - 1);
440 const CoolPropDbl d2DELTA_dTau2 = 2; // d2theta_dTau2 is zero and (dtheta_dtau)^2 = 1
441 const CoolPropDbl d2DELTA_dDelta_dTau = 2 * dtheta_dTau * dtheta_dDelta; // d2theta_dDelta2 is zero
442 const CoolPropDbl d2DELTA_dDelta2 =
443 2 * (theta * d2theta_dDelta2 + POW2(dtheta_dDelta) + Bi * (2 * ai * ai - ai) * pow(POW2(delta - 1.0), ai - 1.0));
444 const CoolPropDbl d3DELTA_dTau3 = 0;
445 const CoolPropDbl d3DELTA_dDelta_dTau2 = 0;
446 const CoolPropDbl d3DELTA_dDelta2_dTau = 2 * dtheta_dTau * d2theta_dDelta2;
447 const CoolPropDbl d3DELTA_dDelta3 = 2
448 * (theta * d3theta_dDelta3 + 3 * dtheta_dDelta * d2theta_dDelta2
449 + 2 * Bi * ai * (2 * ai * ai - 3 * ai + 1) * pow(POW2(delta - 1.0), ai - 1.0) / (delta - 1));
450
451 const CoolPropDbl d4DELTA_dTau4 = 0;
452 const CoolPropDbl d4DELTA_dDelta_dTau3 = 0;
453 const CoolPropDbl d4DELTA_dDelta2_dTau2 = 0;
454 const CoolPropDbl d4DELTA_dDelta3_dTau = 2 * dtheta_dTau * d3theta_dDelta3;
455 const CoolPropDbl d4DELTA_dDelta4 = 2
456 * (theta * d4theta_dDelta4 + 4 * dtheta_dDelta * d3theta_dDelta3 + 3 * POW2(d2theta_dDelta2)
457 + 2 * Bi * ai * (4 * ai * ai * ai - 12 * ai * ai + 11 * ai - 3) * pow(POW2(delta - 1.0), ai - 2.0));
458
459 const CoolPropDbl dDELTAbi_dDelta = bi * pow(DELTA, bi - 1.0) * dDELTA_dDelta;
460 const CoolPropDbl dDELTAbi_dTau = -2.0 * theta * bi * pow(DELTA, bi - 1.0);
461 const CoolPropDbl d2DELTAbi_dDelta2 = bi * (pow(DELTA, bi - 1) * d2DELTA_dDelta2 + (bi - 1.0) * pow(DELTA, bi - 2.0) * pow(dDELTA_dDelta, 2));
462 const CoolPropDbl d3DELTAbi_dDelta3 =
463 bi
464 * (pow(DELTA, bi - 1) * d3DELTA_dDelta3 + d2DELTA_dDelta2 * (bi - 1) * pow(DELTA, bi - 2) * dDELTA_dDelta
465 + (bi - 1)
466 * (pow(DELTA, bi - 2) * 2 * dDELTA_dDelta * d2DELTA_dDelta2
467 + pow(dDELTA_dDelta, 2) * (bi - 2) * pow(DELTA, bi - 3) * dDELTA_dDelta));
468 const CoolPropDbl d2DELTAbi_dDelta_dTau =
469 -Ai * bi * 2.0 / betai * pow(DELTA, bi - 1.0) * (delta - 1.0) * pow(pow(delta - 1.0, 2), 1.0 / (2.0 * betai) - 1.0)
470 - 2.0 * theta * bi * (bi - 1.0) * pow(DELTA, bi - 2.0) * dDELTA_dDelta;
471 const CoolPropDbl d2DELTAbi_dTau2 = 2.0 * bi * pow(DELTA, bi - 1.0) + 4.0 * pow(theta, 2) * bi * (bi - 1.0) * pow(DELTA, bi - 2.0);
472 const CoolPropDbl d3DELTAbi_dTau3 =
473 -12.0 * theta * bi * (bi - 1.0) * pow(DELTA, bi - 2) - 8 * pow(theta, 3) * bi * (bi - 1) * (bi - 2) * pow(DELTA, bi - 3);
474 const CoolPropDbl d3DELTAbi_dDelta_dTau2 = 2 * bi * (bi - 1) * pow(DELTA, bi - 2) * dDELTA_dDelta
475 + 4 * pow(theta, 2) * bi * (bi - 1) * (bi - 2) * pow(DELTA, bi - 3) * dDELTA_dDelta
476 + 8 * theta * bi * (bi - 1) * pow(DELTA, bi - 2) * dtheta_dDelta;
477 const CoolPropDbl d3DELTAbi_dDelta2_dTau =
478 bi
479 * ((bi - 1) * pow(DELTA, bi - 2) * dDELTA_dTau * d2DELTA_dDelta2 + pow(DELTA, bi - 1) * d3DELTA_dDelta2_dTau
480 + (bi - 1)
481 * ((bi - 2) * pow(DELTA, bi - 3) * dDELTA_dTau * pow(dDELTA_dDelta, 2)
482 + pow(DELTA, bi - 2) * 2 * dDELTA_dDelta * d2DELTA_dDelta_dTau));
483
484 // Fourth partials
485 const CoolPropDbl DELTA_bi = pow(DELTA, bi);
486 const CoolPropDbl d4DELTAbi_dTau4 =
487 bi * DELTA_bi / DELTA
488 * ((POW3(bi) - 6 * POW2(bi) + 11 * bi - 6) * POW4(dDELTA_dTau) / POW3(DELTA)
489 + 6 * (bi * bi - 3 * bi + 2) * POW2(dDELTA_dTau / DELTA) * d2DELTA_dTau2 + 4 * (bi - 1) * dDELTA_dTau / DELTA * d3DELTA_dTau3
490 + 3 * (bi - 1) * POW2(d2DELTA_dTau2) / DELTA + d4DELTA_dTau4);
491 const CoolPropDbl d4DELTAbi_dDelta4 =
492 bi * DELTA_bi / DELTA
493 * ((POW3(bi) - 6 * POW2(bi) + 11 * bi - 6) * POW4(dDELTA_dDelta) / POW3(DELTA)
494 + 6 * (bi * bi - 3 * bi + 2) * POW2(dDELTA_dDelta / DELTA) * d2DELTA_dDelta2 + 4 * (bi - 1) * dDELTA_dDelta / DELTA * d3DELTA_dDelta3
495 + 3 * (bi - 1) * POW2(d2DELTA_dDelta2) / DELTA + d4DELTA_dDelta4);
496 const CoolPropDbl d4DELTAbi_dDelta_dTau3 =
497 bi * (bi - 1) * DELTA_bi / POW2(DELTA) * dDELTA_dDelta
498 * ((bi - 1) * (bi - 2) * POW3(dDELTA_dTau) / POW2(DELTA) + 3 * (bi - 1) * dDELTA_dTau / DELTA * d2DELTA_dTau2 + d3DELTA_dTau3)
499 + bi * DELTA_bi / DELTA
500 * ((bi - 1) * (bi - 2)
501 * (3 * POW2(dDELTA_dTau) / POW2(DELTA) * d2DELTA_dDelta_dTau - 2 * POW3(dDELTA_dTau) / POW3(DELTA) * dDELTA_dDelta)
502 + 3 * (bi - 1)
503 * (dDELTA_dTau / DELTA * d3DELTA_dDelta_dTau2 + d2DELTA_dTau2 / DELTA * d2DELTA_dDelta_dTau
504 - dDELTA_dDelta / POW2(DELTA) * dDELTA_dTau * d2DELTA_dTau2)
505 + d4DELTA_dDelta_dTau3);
506 const CoolPropDbl d4DELTAbi_dDelta3_dTau =
507 bi * (bi - 1) * DELTA_bi / POW2(DELTA) * dDELTA_dTau
508 * ((bi - 1) * (bi - 2) * POW3(dDELTA_dDelta) / POW2(DELTA) + 3 * (bi - 1) * dDELTA_dDelta / DELTA * d2DELTA_dDelta2 + d3DELTA_dDelta3)
509 + bi * DELTA_bi / DELTA
510 * ((bi - 1) * (bi - 2)
511 * (3 * POW2(dDELTA_dDelta) / POW2(DELTA) * d2DELTA_dDelta_dTau - 2 * POW3(dDELTA_dDelta) / POW3(DELTA) * dDELTA_dTau)
512 + 3 * (bi - 1)
513 * (dDELTA_dDelta / DELTA * d3DELTA_dDelta2_dTau + d2DELTA_dDelta2 / DELTA * d2DELTA_dDelta_dTau
514 - dDELTA_dTau / POW2(DELTA) * dDELTA_dDelta * d2DELTA_dDelta2)
515 + d4DELTA_dDelta3_dTau);
516 const CoolPropDbl d4DELTAbi_dDelta2_dTau2 = bi * DELTA_bi / POW4(DELTA)
517 * ((POW3(bi) - 6 * bi * bi + 11 * bi - 6) * POW2(dDELTA_dDelta) * POW2(dDELTA_dTau) // Yellow
518 + (bi - 1) * (bi - 2) * DELTA * POW2(dDELTA_dDelta) * d2DELTA_dTau2 // Orange
519 + (bi - 1) * (bi - 2) * DELTA * POW2(dDELTA_dTau) * d2DELTA_dDelta2 // Pink
520 + 4 * (bi - 1) * (bi - 2) * DELTA * dDELTA_dDelta * dDELTA_dTau * d2DELTA_dDelta_dTau // Green
521 + 2 * (bi - 1) * POW2(DELTA * d2DELTA_dDelta_dTau) // Blue hi
522 + 2 * (bi - 1) * POW2(DELTA) * dDELTA_dTau * d3DELTA_dDelta2_dTau // Blue sharp
523 + 2 * (bi - 1) * POW2(DELTA) * dDELTA_dDelta * d3DELTA_dDelta_dTau2 // Red sharp
524 + (bi - 1) * POW2(DELTA) * d2DELTA_dDelta2 * d2DELTA_dTau2 // black sharp
525 + POW3(DELTA) * d4DELTA_dDelta2_dTau2);
526
527 derivs.alphar += delta * ni * DELTA_bi * PSI;
528
529 // First partials
530 derivs.dalphar_dtau += ni * delta * (DELTA_bi * dPSI_dTau + dDELTAbi_dTau * PSI);
531 derivs.dalphar_ddelta += ni * (DELTA_bi * (PSI + delta * dPSI_dDelta) + dDELTAbi_dDelta * delta * PSI);
532
533 // Second partials
534 derivs.d2alphar_dtau2 += ni * delta * (d2DELTAbi_dTau2 * PSI + 2 * dDELTAbi_dTau * dPSI_dTau + DELTA_bi * d2PSI_dTau2);
535 derivs.d2alphar_ddelta_dtau += ni
536 * (DELTA_bi * (dPSI_dTau + delta * d2PSI_dDelta_dTau) + delta * dDELTAbi_dDelta * dPSI_dTau
537 + dDELTAbi_dTau * (PSI + delta * dPSI_dDelta) + d2DELTAbi_dDelta_dTau * delta * PSI);
538 derivs.d2alphar_ddelta2 += ni
539 * (DELTA_bi * (2.0 * dPSI_dDelta + delta * d2PSI_dDelta2) + 2.0 * dDELTAbi_dDelta * (PSI + delta * dPSI_dDelta)
540 + d2DELTAbi_dDelta2 * delta * PSI);
541
542 // Third partials
543 derivs.d3alphar_dtau3 +=
544 ni * delta * (d3DELTAbi_dTau3 * PSI + 3 * d2DELTAbi_dTau2 * dPSI_dTau + 3 * dDELTAbi_dTau * d2PSI_dTau2 + DELTA_bi * d3PSI_dTau3);
545 derivs.d3alphar_ddelta_dtau2 +=
546 ni * delta
547 * (d2DELTAbi_dTau2 * dPSI_dDelta + d3DELTAbi_dDelta_dTau2 * PSI + 2 * dDELTAbi_dTau * d2PSI_dDelta_dTau
548 + 2.0 * d2DELTAbi_dDelta_dTau * dPSI_dTau + DELTA_bi * d3PSI_dDelta_dTau2 + dDELTAbi_dDelta * d2PSI_dTau2)
549 + ni * (d2DELTAbi_dTau2 * PSI + 2.0 * dDELTAbi_dTau * dPSI_dTau + DELTA_bi * d2PSI_dTau2);
550 derivs.d3alphar_ddelta3 +=
551 ni
552 * (DELTA_bi * (3 * d2PSI_dDelta2 + delta * d3PSI_dDelta3) + 3 * dDELTAbi_dDelta * (2 * dPSI_dDelta + delta * d2PSI_dDelta2)
553 + 3 * d2DELTAbi_dDelta2 * (PSI + delta * dPSI_dDelta) + d3DELTAbi_dDelta3 * PSI * delta);
554 CoolPropDbl Line1 =
555 DELTA_bi * (2 * d2PSI_dDelta_dTau + delta * d3PSI_dDelta2_dTau) + dDELTAbi_dTau * (2 * dPSI_dDelta + delta * d2PSI_dDelta2);
556 CoolPropDbl Line2 = 2 * dDELTAbi_dDelta * (dPSI_dTau + delta * d2PSI_dDelta_dTau) + 2 * d2DELTAbi_dDelta_dTau * (PSI + delta * dPSI_dDelta);
557 CoolPropDbl Line3 = d2DELTAbi_dDelta2 * delta * dPSI_dTau + d3DELTAbi_dDelta2_dTau * delta * PSI;
558 derivs.d3alphar_ddelta2_dtau += ni * (Line1 + Line2 + Line3);
559
560 // Fourth partials
561 derivs.d4alphar_dtau4 += ni * delta
562 * (DELTA_bi * d4PSI_dTau4 + 4 * dDELTAbi_dTau * d3PSI_dTau3 + 6 * d2DELTAbi_dTau2 * d2PSI_dTau2
563 + 4 * d3DELTAbi_dTau3 * dPSI_dTau + PSI * d4DELTAbi_dTau4);
564 derivs.d4alphar_ddelta4 +=
565 ni
566 * (delta * DELTA_bi * d4PSI_dDelta4 + delta * PSI * d4DELTAbi_dDelta4 + 4 * delta * dDELTAbi_dDelta * d3PSI_dDelta3
567 + 4 * delta * dPSI_dDelta * d3DELTAbi_dDelta3 + 6 * delta * d2DELTAbi_dDelta2 * d2PSI_dDelta2 + 4 * DELTA_bi * d3PSI_dDelta3
568 + 4 * PSI * d3DELTAbi_dDelta3 + 12 * dDELTAbi_dDelta * d2PSI_dDelta2 + 12 * dPSI_dDelta * d2DELTAbi_dDelta2);
569
570 derivs.d4alphar_ddelta_dtau3 +=
571 ni
572 * (delta * DELTA_bi * d4PSI_dDelta_dTau3 + delta * PSI * d4DELTAbi_dDelta_dTau3 + delta * dDELTAbi_dDelta * d3PSI_dTau3
573 + 3 * delta * dDELTAbi_dTau * d3PSI_dDelta_dTau2 + delta * dPSI_dDelta * d3DELTAbi_dTau3 + 3 * delta * dPSI_dTau * d3DELTAbi_dDelta_dTau2
574 + 3 * delta * d2DELTAbi_dDelta_dTau * d2PSI_dTau2 + 3 * delta * d2DELTAbi_dTau2 * d2PSI_dDelta_dTau + DELTA_bi * d3PSI_dTau3
575 + PSI * d3DELTAbi_dTau3 + 3 * dDELTAbi_dTau * d2PSI_dTau2 + 3 * dPSI_dTau * d2DELTAbi_dTau2);
576 derivs.d4alphar_ddelta3_dtau +=
577 ni
578 * (delta * DELTA_bi * d4PSI_dDelta3_dTau + delta * PSI * d4DELTAbi_dDelta3_dTau + 3 * delta * dDELTAbi_dDelta * d3PSI_dDelta2_dTau
579 + delta * dDELTAbi_dTau * d3PSI_dDelta3 + 3 * delta * dPSI_dDelta * d3DELTAbi_dDelta2_dTau + delta * dPSI_dTau * d3DELTAbi_dDelta3
580 + 3 * delta * d2DELTAbi_dDelta2 * d2PSI_dDelta_dTau + 3 * delta * d2DELTAbi_dDelta_dTau * d2PSI_dDelta2
581 + 3 * DELTA_bi * d3PSI_dDelta2_dTau + 3 * PSI * d3DELTAbi_dDelta2_dTau + 6 * dDELTAbi_dDelta * d2PSI_dDelta_dTau
582 + 3 * dDELTAbi_dTau * d2PSI_dDelta2 + 6 * dPSI_dDelta * d2DELTAbi_dDelta_dTau + 3 * dPSI_dTau * d2DELTAbi_dDelta2);
583 derivs.d4alphar_ddelta2_dtau2 +=
584 ni
585 * (delta * DELTA_bi * d4PSI_dDelta2_dTau2 + delta * PSI * d4DELTAbi_dDelta2_dTau2 + 2 * delta * dDELTAbi_dDelta * d3PSI_dDelta_dTau2
586 + 2 * delta * dDELTAbi_dTau * d3PSI_dDelta2_dTau + 2 * delta * dPSI_dDelta * d3DELTAbi_dDelta_dTau2
587 + 2 * delta * dPSI_dTau * d3DELTAbi_dDelta2_dTau + delta * d2DELTAbi_dDelta2 * d2PSI_dTau2
588 + 4 * delta * d2DELTAbi_dDelta_dTau * d2PSI_dDelta_dTau + delta * d2DELTAbi_dTau2 * d2PSI_dDelta2 + 2 * DELTA_bi * d3PSI_dDelta_dTau2
589 + 2 * PSI * d3DELTAbi_dDelta_dTau2 + 2 * dDELTAbi_dDelta * d2PSI_dTau2 + 4 * dDELTAbi_dTau * d2PSI_dDelta_dTau
590 + 2 * dPSI_dDelta * d2DELTAbi_dTau2 + 4 * dPSI_dTau * d2DELTAbi_dDelta_dTau);
591 }
592}
593
594#if ENABLE_CATCH
595mcx::MultiComplex<double> ResidualHelmholtzNonAnalytic::one_mcx(const mcx::MultiComplex<double>& tau,
596 const mcx::MultiComplex<double>& delta) const {
597
598 mcx::MultiComplex<double> sum00 = 0.0 * tau * delta;
599 for (unsigned int i = 0; i < N; ++i) {
600 const ResidualHelmholtzNonAnalyticElement& el = elements[i];
601 const CoolPropDbl ni = el.n, ai = el.a, bi = el.b, betai = el.beta;
602 const CoolPropDbl Ai = el.A, Bi = el.B, Ci = el.C, Di = el.D;
603
604 const auto theta = (1.0 - tau) + Ai * pow(POW2(delta - 1.0), 1.0 / (2.0 * betai));
605 const auto PSI = exp(-Ci * POW2(delta - 1.0) - Di * POW2(tau - 1.0));
606 const auto DELTA = POW2(theta) + Bi * pow(POW2(delta - 1.0), ai);
607 const auto DELTA_bi = pow(DELTA, bi);
608
609 sum00 += delta * ni * DELTA_bi * PSI;
610 }
611 return sum00;
612}
613#endif
614
615void ResidualHelmholtzGeneralizedCubic::all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw() {
616 if (!enabled) {
617 return;
618 }
619
620 derivs.alphar += m_abstractcubic->alphar(tau, delta, z, 0, 0);
621
622 derivs.dalphar_ddelta += m_abstractcubic->alphar(tau, delta, z, 0, 1);
623 derivs.dalphar_dtau += m_abstractcubic->alphar(tau, delta, z, 1, 0);
624
625 derivs.d2alphar_ddelta2 += m_abstractcubic->alphar(tau, delta, z, 0, 2);
626 derivs.d2alphar_ddelta_dtau += m_abstractcubic->alphar(tau, delta, z, 1, 1);
627 derivs.d2alphar_dtau2 += m_abstractcubic->alphar(tau, delta, z, 2, 0);
628
629 derivs.d3alphar_ddelta3 += m_abstractcubic->alphar(tau, delta, z, 0, 3);
630 derivs.d3alphar_ddelta2_dtau += m_abstractcubic->alphar(tau, delta, z, 1, 2);
631 derivs.d3alphar_ddelta_dtau2 += m_abstractcubic->alphar(tau, delta, z, 2, 1);
632 derivs.d3alphar_dtau3 += m_abstractcubic->alphar(tau, delta, z, 3, 0);
633
634 derivs.d4alphar_ddelta4 += m_abstractcubic->alphar(tau, delta, z, 0, 4);
635 derivs.d4alphar_ddelta3_dtau += m_abstractcubic->alphar(tau, delta, z, 1, 3);
636 derivs.d4alphar_ddelta2_dtau2 += m_abstractcubic->alphar(tau, delta, z, 2, 2);
637 derivs.d4alphar_ddelta_dtau3 += m_abstractcubic->alphar(tau, delta, z, 3, 1);
638 derivs.d4alphar_dtau4 += m_abstractcubic->alphar(tau, delta, z, 4, 0);
639}
640
665void ResidualHelmholtzGaoB::all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) {
666 if (!enabled) {
667 return;
668 }
669
670 CoolPropDbl Ftau = 0, Fdelta = 0, taudFtaudtau = 0, tau2d2Ftaudtau2 = 0, tau3d3Ftaudtau3 = 0, tau4d4Ftaudtau4 = 0, deltadFdeltaddelta = 0,
671 delta2d2Fdeltaddelta2 = 0, delta3d3Fdeltaddelta3 = 0, delta4d4Fdeltaddelta4 = 0;
672
673 for (int i = 0; i < static_cast<int>(n.size()); ++i) {
674
675 const CoolPropDbl n = this->n[i], t = this->t[i], d = this->d[i], eta = this->eta[i], beta = this->beta[i], gamma = this->gamma[i],
676 epsilon = this->epsilon[i], b = this->b[i];
677
678 Ftau = pow(tau, t) * exp(1.0 / (b + beta * pow(-gamma + tau, 2)));
679 Fdelta = pow(delta, d) * exp(eta * pow(delta - epsilon, 2));
680 taudFtaudtau = (2 * beta * pow(tau, t + 1) * (gamma - tau) + t * pow(tau, t) * pow(b + beta * pow(gamma - tau, 2), 2))
681 * exp(1.0 / (b + beta * pow(gamma - tau, 2))) / pow(b + beta * pow(gamma - tau, 2), 2);
682 tau2d2Ftaudtau2 = pow(tau, t)
683 * (4 * beta * t * tau * pow(b + beta * pow(gamma - tau, 2), 2) * (gamma - tau)
684 + 2 * beta * pow(tau, 2)
685 * (4 * beta * (b + beta * pow(gamma - tau, 2)) * pow(gamma - tau, 2) + 2 * beta * pow(gamma - tau, 2)
686 - pow(b + beta * pow(gamma - tau, 2), 2))
687 + t * pow(b + beta * pow(gamma - tau, 2), 4) * (t - 1))
688 * exp(1.0 / (b + beta * pow(gamma - tau, 2))) / pow(b + beta * pow(gamma - tau, 2), 4);
689 tau3d3Ftaudtau3 =
690 pow(tau, t)
691 * (4 * pow(beta, 2) * pow(tau, 3) * (gamma - tau)
692 * (12 * beta * (b + beta * pow(gamma - tau, 2)) * pow(gamma - tau, 2) + 2 * beta * pow(gamma - tau, 2)
693 - 6 * pow(b + beta * pow(gamma - tau, 2), 3) + pow(b + beta * pow(gamma - tau, 2), 2) * (12 * beta * pow(gamma - tau, 2) - 3))
694 + 6 * beta * t * pow(tau, 2) * pow(b + beta * pow(gamma - tau, 2), 2)
695 * (4 * beta * (b + beta * pow(gamma - tau, 2)) * pow(gamma - tau, 2) + 2 * beta * pow(gamma - tau, 2)
696 - pow(b + beta * pow(gamma - tau, 2), 2))
697 + 6 * beta * t * tau * pow(b + beta * pow(gamma - tau, 2), 4) * (gamma - tau) * (t - 1)
698 + t * pow(b + beta * pow(gamma - tau, 2), 6) * (pow(t, 2) - 3 * t + 2))
699 * exp(1.0 / (b + beta * pow(gamma - tau, 2))) / pow(b + beta * pow(gamma - tau, 2), 6);
700 tau4d4Ftaudtau4 =
701 pow(tau, t)
702 * (16 * pow(beta, 2) * t * pow(tau, 3) * pow(b + beta * pow(gamma - tau, 2), 2) * (gamma - tau)
703 * (12 * beta * (b + beta * pow(gamma - tau, 2)) * pow(gamma - tau, 2) + 2 * beta * pow(gamma - tau, 2)
704 - 6 * pow(b + beta * pow(gamma - tau, 2), 3) + pow(b + beta * pow(gamma - tau, 2), 2) * (12 * beta * pow(gamma - tau, 2) - 3))
705 + pow(beta, 2) * pow(tau, 4)
706 * (pow(beta, 2) * (192 * b + 192 * beta * pow(gamma - tau, 2)) * pow(gamma - tau, 4) + 16 * pow(beta, 2) * pow(gamma - tau, 4)
707 + 96 * beta * pow(b + beta * pow(gamma - tau, 2), 3) * pow(gamma - tau, 2) * (4 * beta * pow(gamma - tau, 2) - 3)
708 + 48 * beta * pow(b + beta * pow(gamma - tau, 2), 2) * pow(gamma - tau, 2) * (12 * beta * pow(gamma - tau, 2) - 1)
709 + 24 * pow(b + beta * pow(gamma - tau, 2), 5) + pow(b + beta * pow(gamma - tau, 2), 4) * (-288 * beta * pow(gamma - tau, 2) + 12))
710 + 12 * beta * t * pow(tau, 2) * pow(b + beta * pow(gamma - tau, 2), 4) * (t - 1)
711 * (4 * beta * (b + beta * pow(gamma - tau, 2)) * pow(gamma - tau, 2) + 2 * beta * pow(gamma - tau, 2)
712 - pow(b + beta * pow(gamma - tau, 2), 2))
713 + 8 * beta * t * tau * pow(b + beta * pow(gamma - tau, 2), 6) * (gamma - tau) * (pow(t, 2) - 3 * t + 2)
714 + t * pow(b + beta * pow(gamma - tau, 2), 8) * (pow(t, 3) - 6 * pow(t, 2) + 11 * t - 6))
715 * exp(1.0 / (b + beta * pow(gamma - tau, 2))) / pow(b + beta * pow(gamma - tau, 2), 8);
716 deltadFdeltaddelta = (d * pow(delta, d) + 2 * pow(delta, d + 1) * eta * (delta - epsilon)) * exp(eta * pow(delta - epsilon, 2));
717 delta2d2Fdeltaddelta2 =
718 pow(delta, d) * (4 * d * delta * eta * (delta - epsilon) + d * (d - 1) + 2 * pow(delta, 2) * eta * (2 * eta * pow(delta - epsilon, 2) + 1))
719 * exp(eta * pow(delta - epsilon, 2));
720 delta3d3Fdeltaddelta3 =
721 pow(delta, d)
722 * (6 * d * pow(delta, 2) * eta * (2 * eta * pow(delta - epsilon, 2) + 1) + 6 * d * delta * eta * (d - 1) * (delta - epsilon)
723 + d * (pow(d, 2) - 3 * d + 2) + 4 * pow(delta, 3) * pow(eta, 2) * (delta - epsilon) * (2 * eta * pow(delta - epsilon, 2) + 3))
724 * exp(eta * pow(delta - epsilon, 2));
725 delta4d4Fdeltaddelta4 =
726 pow(delta, d)
727 * (16 * d * pow(delta, 3) * pow(eta, 2) * (delta - epsilon) * (2 * eta * pow(delta - epsilon, 2) + 3)
728 + 12 * d * pow(delta, 2) * eta * (d - 1) * (2 * eta * pow(delta - epsilon, 2) + 1)
729 + 8 * d * delta * eta * (delta - epsilon) * (pow(d, 2) - 3 * d + 2) + d * (pow(d, 3) - 6 * pow(d, 2) + 11 * d - 6)
730 + pow(delta, 4) * pow(eta, 2) * (16 * pow(eta, 2) * pow(delta - epsilon, 4) + 48 * eta * pow(delta - epsilon, 2) + 12))
731 * exp(eta * pow(delta - epsilon, 2));
732
733 derivs.alphar += n * Ftau * Fdelta;
734
735 derivs.dalphar_ddelta += n * Ftau * deltadFdeltaddelta / delta;
736 derivs.dalphar_dtau += n * Fdelta * taudFtaudtau / tau;
737
738 derivs.d2alphar_ddelta2 += n * Ftau * delta2d2Fdeltaddelta2 / POW2(delta);
739 derivs.d2alphar_ddelta_dtau += n * taudFtaudtau * deltadFdeltaddelta / tau / delta;
740 derivs.d2alphar_dtau2 += n * Fdelta * tau2d2Ftaudtau2 / POW2(tau);
741
742 derivs.d3alphar_ddelta3 += n * Ftau * delta3d3Fdeltaddelta3 / POW3(delta);
743 derivs.d3alphar_ddelta2_dtau += n * taudFtaudtau * delta2d2Fdeltaddelta2 / POW2(delta) / tau;
744 derivs.d3alphar_ddelta_dtau2 += n * tau2d2Ftaudtau2 * deltadFdeltaddelta / POW2(tau) / delta;
745 derivs.d3alphar_dtau3 += n * Fdelta * tau3d3Ftaudtau3 / POW3(tau);
746
747 derivs.d4alphar_ddelta4 += n * Ftau * delta4d4Fdeltaddelta4 / POW4(delta);
748 derivs.d4alphar_ddelta3_dtau += n * taudFtaudtau * delta3d3Fdeltaddelta3 / POW3(delta) / tau;
749 derivs.d4alphar_ddelta2_dtau2 += n * tau2d2Ftaudtau2 * delta2d2Fdeltaddelta2 / POW2(delta) / POW2(tau);
750 derivs.d4alphar_ddelta_dtau3 += n * tau3d3Ftaudtau3 * deltadFdeltaddelta / POW3(tau) / delta;
751 derivs.d4alphar_dtau4 += n * Fdelta * tau4d4Ftaudtau4 / POW4(tau);
752 }
753}
754
755
756#if ENABLE_CATCH
757mcx::MultiComplex<double> ResidualHelmholtzGaoB::one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const {
758 mcx::MultiComplex<double> sum00 = 0.0;
759 for (std::size_t i = 0; i < static_cast<int>(n.size()); ++i) {
760 auto u = b[i] + beta[i] * (tau-gamma[i])*(tau-gamma[i]);
761 auto Ftau = pow(tau, t[i]) * exp(1.0 / u);
762 auto Fdelta = pow(delta, d[i]) * exp(eta[i] * pow(delta - epsilon[i], 2));
763 sum00 += n[i] * Ftau * Fdelta;
764 }
765 return sum00;
766}
767#endif
768
769ResidualHelmholtzXiangDeiters::ResidualHelmholtzXiangDeiters(const CoolPropDbl Tc, const CoolPropDbl pc, const CoolPropDbl rhomolarc,
770 const CoolPropDbl acentric, const CoolPropDbl R)
771 : Tc(Tc), pc(pc), rhomolarc(rhomolarc), acentric(acentric), R(R) {
772 double Zc = pc / (R * Tc * rhomolarc);
773 theta = POW2(Zc - 0.29);
774
775 // From Xiang & Deiters, doi:10.1016/j.ces.2007.11.029
776 double _d[] = {1, 1, 1, 2, 3, 7, 1, 1, 2, 5, 1, 1, 4, 2};
777 std::vector<CoolPropDbl> d(_d, _d + sizeof(_d) / sizeof(double));
778 double _t[] = {0.25, 1.25, 1.5, 1.375, 0.25, 0.875, 0, 2.375, 2, 2.125, 3.5, 6.5, 4.75, 12.5};
779 std::vector<CoolPropDbl> t(_t, _t + sizeof(_t) / sizeof(double));
780 double _l[] = {0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3};
781 std::vector<CoolPropDbl> l(_l, _l + sizeof(_l) / sizeof(double));
782 double _g[] = {0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1};
783 std::vector<CoolPropDbl> g(_g, _g + sizeof(_g) / sizeof(double));
784 double _a0[] = {8.5740489E-01, -3.2863233E+00, 1.6480939E+00, -5.4524817E-02, 6.1623592E-02, 2.7389266E-04, -6.0655087E-02,
785 -3.1811852E-02, -1.1550422E-01, -1.8610466E-02, -1.8348671E-01, 5.5071325E-03, -1.2268039E-02, -5.0433436E-03};
786 std::vector<CoolPropDbl> a0(_a0, _a0 + sizeof(_a0) / sizeof(double));
787 double _a1[] = {5.6200117E-01, 3.2439544E+00, -4.9628768E+00, -2.2132851E-01, 9.3388356E-02, 2.4940171E-05, -1.7519204E-01,
788 8.9325660E-01, 2.9886613E+00, 1.0881387E-01, -6.7166746E-01, 1.4477326E-01, -2.8716809E-01, -1.1478402E-01};
789 std::vector<CoolPropDbl> a1(_a1, _a1 + sizeof(_a1) / sizeof(double));
790 double _a2[] = {-8.1680511E+01, 4.6384732E+02, -2.7970850E+02, 2.9317364E+01, -2.2324825E+01, -5.0932691E-02, -7.2836590E+00,
791 -2.2063100E+02, -3.0435126E+02, 5.8514719E+00, 1.7995451E+02, -1.0178400E+02, 4.0848053E+01, 1.2411984E+01};
792 std::vector<CoolPropDbl> a2(_a2, _a2 + sizeof(_a2) / sizeof(double));
793
794 phi0.add_Exponential(a0, d, t, g, l);
795 phi1.add_Exponential(a1, d, t, g, l);
796 phi2.add_Exponential(a2, d, t, g, l);
797
798 enabled = true;
799};
800
801void ResidualHelmholtzXiangDeiters::all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) {
802 if (!enabled) {
803 return;
804 }
805
806 HelmholtzDerivatives derivs0, derivs1, derivs2;
807
808 // Calculate each of the derivative terms
809 phi0.all(tau, delta, derivs0);
810 phi1.all(tau, delta, derivs1);
811 phi2.all(tau, delta, derivs2);
812
813 // Add up the contributions
814 derivs = derivs + derivs0 + derivs1 * acentric + derivs2 * theta;
815}
816#if ENABLE_CATCH
817mcx::MultiComplex<double> ResidualHelmholtzXiangDeiters::one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const {
818
819 return phi0.one_mcx(tau, delta) + phi1.one_mcx(tau, delta) * acentric + phi2.one_mcx(tau, delta) * theta;
820}
821#endif
822
823void ResidualHelmholtzSAFTAssociating::to_json(rapidjson::Value& el, rapidjson::Document& doc) {
824 el.AddMember("type", "ResidualHelmholtzSAFTAssociating", doc.GetAllocator());
825 el.AddMember("a", a, doc.GetAllocator());
826 el.AddMember("m", m, doc.GetAllocator());
827 el.AddMember("epsilonbar", epsilonbar, doc.GetAllocator());
828 el.AddMember("vbarn", vbarn, doc.GetAllocator());
829 el.AddMember("kappabar", kappabar, doc.GetAllocator());
830}
831CoolPropDbl ResidualHelmholtzSAFTAssociating::Deltabar(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
832 return this->g(this->eta(delta)) * (exp(this->epsilonbar * tau) - 1) * this->kappabar;
833}
834CoolPropDbl ResidualHelmholtzSAFTAssociating::dDeltabar_ddelta__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
835 return this->dg_deta(this->eta(delta)) * (exp(this->epsilonbar * tau) - 1) * this->kappabar * this->vbarn;
836}
837CoolPropDbl ResidualHelmholtzSAFTAssociating::d2Deltabar_ddelta2__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
838 return this->d2g_deta2(this->eta(delta)) * (exp(this->epsilonbar * tau) - 1) * this->kappabar * pow(this->vbarn, (int)2);
839}
840CoolPropDbl ResidualHelmholtzSAFTAssociating::dDeltabar_dtau__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
841 return this->g(this->eta(delta)) * this->kappabar * exp(this->epsilonbar * tau) * this->epsilonbar;
842}
843CoolPropDbl ResidualHelmholtzSAFTAssociating::d2Deltabar_dtau2__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
844 return this->g(this->eta(delta)) * this->kappabar * exp(this->epsilonbar * tau) * pow(this->epsilonbar, (int)2);
845}
846CoolPropDbl ResidualHelmholtzSAFTAssociating::d2Deltabar_ddelta_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
847 return this->dg_deta(this->eta(delta)) * exp(this->epsilonbar * tau) * this->epsilonbar * this->kappabar * this->vbarn;
848}
849CoolPropDbl ResidualHelmholtzSAFTAssociating::d3Deltabar_dtau3__constdelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
850 return this->g(this->eta(delta)) * this->kappabar * exp(this->epsilonbar * tau) * pow(this->epsilonbar, (int)3);
851}
852CoolPropDbl ResidualHelmholtzSAFTAssociating::d3Deltabar_ddelta_dtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
853 return this->dg_deta(this->eta(delta)) * this->kappabar * exp(this->epsilonbar * tau) * pow(this->epsilonbar, (int)2) * this->vbarn;
854}
855CoolPropDbl ResidualHelmholtzSAFTAssociating::d3Deltabar_ddelta2_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
856 return this->d2g_deta2(this->eta(delta)) * exp(this->epsilonbar * tau) * this->epsilonbar * this->kappabar * pow(this->vbarn, (int)2);
857}
858CoolPropDbl ResidualHelmholtzSAFTAssociating::d3Deltabar_ddelta3__consttau(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
859 return this->d3g_deta3(this->eta(delta)) * (exp(this->epsilonbar * tau) - 1) * this->kappabar * pow(this->vbarn, (int)3);
860}
861
862CoolPropDbl ResidualHelmholtzSAFTAssociating::X(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const {
863 return 2 / (sqrt(1 + 4 * Deltabar * delta) + 1);
864}
865CoolPropDbl ResidualHelmholtzSAFTAssociating::dX_dDeltabar__constdelta(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const {
866 CoolPropDbl X = this->X(delta, Deltabar);
867 return -delta * X * X / (2 * Deltabar * delta * X + 1);
868}
869CoolPropDbl ResidualHelmholtzSAFTAssociating::dX_ddelta__constDeltabar(const CoolPropDbl& delta, const CoolPropDbl& Deltabar) const {
870 CoolPropDbl X = this->X(delta, Deltabar);
871 return -Deltabar * X * X / (2 * Deltabar * delta * X + 1);
872}
873CoolPropDbl ResidualHelmholtzSAFTAssociating::dX_dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
874 CoolPropDbl Deltabar = this->Deltabar(tau, delta);
875 return this->dX_dDeltabar__constdelta(delta, Deltabar) * this->dDeltabar_dtau__constdelta(tau, delta);
876}
877CoolPropDbl ResidualHelmholtzSAFTAssociating::dX_ddelta(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
878 CoolPropDbl Deltabar = this->Deltabar(tau, delta);
879 return (this->dX_ddelta__constDeltabar(delta, Deltabar)
880 + this->dX_dDeltabar__constdelta(delta, Deltabar) * this->dDeltabar_ddelta__consttau(tau, delta));
881}
882CoolPropDbl ResidualHelmholtzSAFTAssociating::d2X_dtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
883 CoolPropDbl Deltabar = this->Deltabar(tau, delta);
884 CoolPropDbl X = this->X(delta, Deltabar);
885 CoolPropDbl beta = this->dDeltabar_dtau__constdelta(tau, delta);
886 CoolPropDbl d_dXdtau_dbeta = -delta * X * X / (2 * Deltabar * delta * X + 1);
887 CoolPropDbl d_dXdtau_dDeltabar = 2 * delta * delta * X * X * X / pow(2 * Deltabar * delta * X + 1, 2) * beta;
888 CoolPropDbl d_dXdtau_dX = -2 * beta * delta * X * (Deltabar * delta * X + 1) / pow(2 * Deltabar * delta * X + 1, 2);
889 CoolPropDbl dbeta_dtau = this->d2Deltabar_dtau2__constdelta(tau, delta);
890 CoolPropDbl dX_dDeltabar = this->dX_dDeltabar__constdelta(delta, Deltabar);
891 return d_dXdtau_dX * dX_dDeltabar * beta + d_dXdtau_dDeltabar * beta + d_dXdtau_dbeta * dbeta_dtau;
892}
893CoolPropDbl ResidualHelmholtzSAFTAssociating::d2X_ddeltadtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
894 CoolPropDbl Deltabar = this->Deltabar(tau, delta);
895 CoolPropDbl X = this->X(delta, Deltabar);
896 CoolPropDbl alpha = this->dDeltabar_ddelta__consttau(tau, delta);
897 CoolPropDbl beta = this->dDeltabar_dtau__constdelta(tau, delta);
898 CoolPropDbl dalpha_dtau = this->d2Deltabar_ddelta_dtau(tau, delta);
899 CoolPropDbl d_dXddelta_dDeltabar = X * X * (2 * delta * delta * X * alpha - 1) / pow(2 * Deltabar * delta * X + 1, 2);
900 CoolPropDbl d_dXddelta_dalpha = -delta * X * X / (2 * Deltabar * delta * X + 1);
901 CoolPropDbl d_dXddelta_dX = -(Deltabar + delta * alpha) * 2 * (Deltabar * delta * X * X + X) / pow(2 * Deltabar * delta * X + 1, 2);
902 CoolPropDbl dX_dDeltabar = this->dX_dDeltabar__constdelta(delta, Deltabar);
903 return d_dXddelta_dX * dX_dDeltabar * beta + d_dXddelta_dDeltabar * beta + d_dXddelta_dalpha * dalpha_dtau;
904}
905CoolPropDbl ResidualHelmholtzSAFTAssociating::d2X_ddelta2(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
906 CoolPropDbl Deltabar = this->Deltabar(tau, delta);
907 CoolPropDbl X = this->X(delta, Deltabar);
908 CoolPropDbl alpha = this->dDeltabar_ddelta__consttau(tau, delta);
909 CoolPropDbl dalpha_ddelta = this->d2Deltabar_ddelta2__consttau(tau, delta);
910
911 CoolPropDbl dX_ddelta_constall = X * X * (2 * Deltabar * Deltabar * X - alpha) / pow(2 * Deltabar * delta * X + 1, 2);
912 CoolPropDbl d_dXddelta_dX = -(Deltabar + delta * alpha) * 2 * (Deltabar * delta * X * X + X) / pow(2 * Deltabar * delta * X + 1, 2);
913 CoolPropDbl d_dXddelta_dDeltabar = X * X * (2 * delta * delta * X * alpha - 1) / pow(2 * Deltabar * delta * X + 1, 2);
914 CoolPropDbl d_dXddelta_dalpha = -delta * X * X / (2 * Deltabar * delta * X + 1);
915
916 CoolPropDbl dX_dDeltabar = this->dX_dDeltabar__constdelta(delta, Deltabar);
917 CoolPropDbl dX_ddelta = this->dX_ddelta__constDeltabar(delta, Deltabar);
918
919 CoolPropDbl val = (dX_ddelta_constall + d_dXddelta_dX * dX_ddelta + d_dXddelta_dX * dX_dDeltabar * alpha + d_dXddelta_dDeltabar * alpha
920 + d_dXddelta_dalpha * dalpha_ddelta);
921 return val;
922}
923CoolPropDbl ResidualHelmholtzSAFTAssociating::d3X_dtau3(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
924 CoolPropDbl Delta = this->Deltabar(tau, delta);
925 CoolPropDbl X = this->X(delta, Delta);
926 CoolPropDbl dX_dDelta = this->dX_dDeltabar__constdelta(delta, Delta);
927 CoolPropDbl Delta_t = this->dDeltabar_dtau__constdelta(tau, delta);
928 CoolPropDbl Delta_tt = this->d2Deltabar_dtau2__constdelta(tau, delta);
929 CoolPropDbl Delta_ttt = this->d3Deltabar_dtau3__constdelta(tau, delta);
930 CoolPropDbl dXtt_dX = 2 * X * delta
931 * (-6 * Delta * pow(Delta_t, 2) * pow(X, 2) * pow(delta, 2) * (Delta * X * delta + 1)
932 + 3 * pow(Delta_t, 2) * X * delta * (2 * Delta * X * delta + 1) - Delta_tt * pow(2 * Delta * X * delta + 1, 3)
933 + X * delta * (Delta * Delta_tt + 3 * pow(Delta_t, 2)) * pow(2 * Delta * X * delta + 1, 2))
934 / pow(2 * Delta * X * delta + 1, 4);
935 CoolPropDbl dXtt_dDelta = 2 * pow(X, 3) * pow(delta, 2)
936 * (-6 * pow(Delta_t, 2) * X * delta * (Delta * X * delta + 1)
937 - 3 * pow(Delta_t, 2) * X * delta * (2 * Delta * X * delta + 1) + Delta_tt * pow(2 * Delta * X * delta + 1, 2))
938 / pow(2 * Delta * X * delta + 1, 4);
939 CoolPropDbl dXtt_dDelta_t = 4 * Delta_t * pow(X, 3) * pow(delta, 2) * (3 * Delta * X * delta + 2) / pow(2 * Delta * X * delta + 1, 3);
940 CoolPropDbl dXtt_dDelta_tt = -pow(X, 2) * delta / (2 * Delta * X * delta + 1);
941 return dXtt_dX * dX_dDelta * Delta_t + dXtt_dDelta * Delta_t + dXtt_dDelta_t * Delta_tt + dXtt_dDelta_tt * Delta_ttt;
942}
943CoolPropDbl ResidualHelmholtzSAFTAssociating::d3X_ddeltadtau2(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
944 CoolPropDbl Delta = this->Deltabar(tau, delta);
945 CoolPropDbl X = this->X(delta, Delta);
946 CoolPropDbl dX_ddelta = this->dX_ddelta__constDeltabar(delta, Delta);
947 CoolPropDbl dX_dDelta = this->dX_dDeltabar__constdelta(delta, Delta);
948 CoolPropDbl Delta_t = this->dDeltabar_dtau__constdelta(tau, delta);
949 CoolPropDbl Delta_d = this->dDeltabar_ddelta__consttau(tau, delta);
950 CoolPropDbl Delta_dt = this->d2Deltabar_ddelta_dtau(tau, delta);
951 CoolPropDbl Delta_tt = this->d2Deltabar_dtau2__constdelta(tau, delta);
952 CoolPropDbl Delta_dtt = this->d3Deltabar_ddelta_dtau2(tau, delta);
953 CoolPropDbl dXtt_ddelta =
954 pow(X, 2)
955 * (-12 * Delta * pow(Delta_t, 2) * pow(X, 2) * pow(delta, 2) * (Delta * X * delta + 1)
956 + 2 * pow(Delta_t, 2) * X * delta * (-Delta * X * delta + 2) * (2 * Delta * X * delta + 1) - Delta_tt * pow(2 * Delta * X * delta + 1, 3)
957 + 2 * X * delta * (Delta * Delta_tt + 2 * pow(Delta_t, 2)) * pow(2 * Delta * X * delta + 1, 2))
958 / pow(2 * Delta * X * delta + 1, 4);
959 CoolPropDbl dXtt_dX = 2 * X * delta
960 * (-6 * Delta * pow(Delta_t, 2) * pow(X, 2) * pow(delta, 2) * (Delta * X * delta + 1)
961 + 3 * pow(Delta_t, 2) * X * delta * (2 * Delta * X * delta + 1) - Delta_tt * pow(2 * Delta * X * delta + 1, 3)
962 + X * delta * (Delta * Delta_tt + 3 * pow(Delta_t, 2)) * pow(2 * Delta * X * delta + 1, 2))
963 / pow(2 * Delta * X * delta + 1, 4);
964 CoolPropDbl dXtt_dDelta = 2 * pow(X, 3) * pow(delta, 2)
965 * (-6 * pow(Delta_t, 2) * X * delta * (Delta * X * delta + 1)
966 - 3 * pow(Delta_t, 2) * X * delta * (2 * Delta * X * delta + 1) + Delta_tt * pow(2 * Delta * X * delta + 1, 2))
967 / pow(2 * Delta * X * delta + 1, 4);
968 CoolPropDbl dXtt_dDelta_t = 4 * Delta_t * pow(X, 3) * pow(delta, 2) * (3 * Delta * X * delta + 2) / pow(2 * Delta * X * delta + 1, 3);
969 CoolPropDbl dXtt_dDelta_tt = -pow(X, 2) * delta / (2 * Delta * X * delta + 1);
970 return dXtt_ddelta + dXtt_dX * dX_ddelta + dXtt_dX * dX_dDelta * Delta_d + dXtt_dDelta * Delta_d + dXtt_dDelta_t * Delta_dt
971 + dXtt_dDelta_tt * Delta_dtt;
972}
973
974CoolPropDbl ResidualHelmholtzSAFTAssociating::d3X_ddelta2dtau(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
975 CoolPropDbl Delta = this->Deltabar(tau, delta);
976 CoolPropDbl X = this->X(delta, Delta);
977 CoolPropDbl dX_dDelta = this->dX_dDeltabar__constdelta(delta, Delta);
978 CoolPropDbl Delta_t = this->dDeltabar_dtau__constdelta(tau, delta);
979 CoolPropDbl Delta_d = this->dDeltabar_ddelta__consttau(tau, delta);
980 CoolPropDbl Delta_dd = this->d2Deltabar_ddelta2__consttau(tau, delta);
981 CoolPropDbl Delta_dt = this->d2Deltabar_ddelta_dtau(tau, delta);
982 CoolPropDbl Delta_ddt = this->d3Deltabar_ddelta2_dtau(tau, delta);
983 CoolPropDbl dXdd_dX =
984 2 * X
985 * (-6 * Delta * pow(X, 2) * delta * pow(Delta + Delta_d * delta, 2) * (Delta * X * delta + 1)
986 - Delta_dd * delta * pow(2 * Delta * X * delta + 1, 3)
987 + 2 * X * (2 * Delta * X * delta + 1)
988 * (-Delta * Delta_d * delta * (2 * Delta_d * X * pow(delta, 2) - 1) - Delta * delta * (2 * pow(Delta, 2) * X - Delta_d)
989 + Delta * (Delta + Delta_d * delta) * (Delta * X * delta + 1) + Delta_d * delta * (Delta + Delta_d * delta) * (Delta * X * delta + 1))
990 + pow(2 * Delta * X * delta + 1, 2)
991 * (3 * pow(Delta, 2) * X + Delta * Delta_dd * X * pow(delta, 2) + Delta * X * (Delta + Delta_d * delta)
992 + pow(Delta_d, 2) * X * pow(delta, 2) + Delta_d * X * delta * (Delta + Delta_d * delta)
993 + Delta_d * (2 * Delta_d * X * pow(delta, 2) - 1) - Delta_d))
994 / pow(2 * Delta * X * delta + 1, 4);
995 CoolPropDbl dXdd_dDelta = pow(X, 3)
996 * (-8 * pow(Delta, 2) * Delta_d * pow(X, 2) * pow(delta, 3) + 8 * pow(Delta, 2) * Delta_dd * pow(X, 2) * pow(delta, 4)
997 + 10 * pow(Delta, 2) * X * delta - 24 * Delta * pow(Delta_d, 2) * pow(X, 2) * pow(delta, 4)
998 + 8 * Delta * Delta_d * X * pow(delta, 2) + 8 * Delta * Delta_dd * X * pow(delta, 3) + 8 * Delta
999 - 18 * pow(Delta_d, 2) * X * pow(delta, 3) + 12 * Delta_d * delta + 2 * Delta_dd * pow(delta, 2))
1000 / (16 * pow(Delta, 4) * pow(X, 4) * pow(delta, 4) + 32 * pow(Delta, 3) * pow(X, 3) * pow(delta, 3)
1001 + 24 * pow(Delta, 2) * pow(X, 2) * pow(delta, 2) + 8 * Delta * X * delta + 1);
1002 CoolPropDbl dXdd_dDelta_d =
1003 2 * pow(X, 2)
1004 * (2 * X * delta * (Delta + Delta_d * delta) * (Delta * X * delta + 1) + (2 * Delta * X * delta + 1) * (2 * Delta_d * X * pow(delta, 2) - 1))
1005 / pow(2 * Delta * X * delta + 1, 3);
1006 CoolPropDbl dXdd_dDelta_dd = -pow(X, 2) * delta / (2 * Delta * X * delta + 1);
1007
1008 return dXdd_dX * dX_dDelta * Delta_t + dXdd_dDelta * Delta_t + dXdd_dDelta_d * Delta_dt + dXdd_dDelta_dd * Delta_ddt;
1009}
1010
1011double Xdd(double X, double delta, double Delta, double Delta_d, double Delta_dd) {
1012 return Delta * pow(X, 2) * (2 * Delta + 2 * Delta_d * delta) * (Delta * pow(X, 2) * delta + X) / pow(2 * Delta * X * delta + 1, 3)
1013 + Delta_d * pow(X, 2) * delta * (2 * Delta + 2 * Delta_d * delta) * (Delta * pow(X, 2) * delta + X) / pow(2 * Delta * X * delta + 1, 3)
1014 + Delta_d * pow(X, 2) * (2 * Delta_d * X * pow(delta, 2) - 1) / pow(2 * Delta * X * delta + 1, 2)
1015 - Delta_dd * pow(X, 2) * delta / (2 * Delta * X * delta + 1)
1016 + pow(X, 2) * (2 * pow(Delta, 2) * X - Delta_d) / pow(2 * Delta * X * delta + 1, 2);
1017}
1018
1019CoolPropDbl ResidualHelmholtzSAFTAssociating::d3X_ddelta3(const CoolPropDbl& tau, const CoolPropDbl& delta) const {
1020 CoolPropDbl Delta = this->Deltabar(tau, delta);
1021 CoolPropDbl X = this->X(delta, Delta);
1022 CoolPropDbl dX_ddelta = this->dX_ddelta__constDeltabar(delta, Delta);
1023 CoolPropDbl dX_dDelta = this->dX_dDeltabar__constdelta(delta, Delta);
1024 CoolPropDbl Delta_d = this->dDeltabar_ddelta__consttau(tau, delta);
1025 CoolPropDbl Delta_dd = this->d2Deltabar_ddelta2__consttau(tau, delta);
1026 CoolPropDbl Delta_ddd = this->d3Deltabar_ddelta3__consttau(tau, delta);
1027
1028 CoolPropDbl dXdd_dX =
1029 2 * X
1030 * (-6 * Delta * pow(X, 2) * delta * pow(Delta + Delta_d * delta, 2) * (Delta * X * delta + 1)
1031 - Delta_dd * delta * pow(2 * Delta * X * delta + 1, 3)
1032 + 2 * X * (2 * Delta * X * delta + 1)
1033 * (-Delta * Delta_d * delta * (2 * Delta_d * X * pow(delta, 2) - 1) - Delta * delta * (2 * pow(Delta, 2) * X - Delta_d)
1034 + Delta * (Delta + Delta_d * delta) * (Delta * X * delta + 1) + Delta_d * delta * (Delta + Delta_d * delta) * (Delta * X * delta + 1))
1035 + pow(2 * Delta * X * delta + 1, 2)
1036 * (3 * pow(Delta, 2) * X + Delta * Delta_dd * X * pow(delta, 2) + Delta * X * (Delta + Delta_d * delta)
1037 + pow(Delta_d, 2) * X * pow(delta, 2) + Delta_d * X * delta * (Delta + Delta_d * delta)
1038 + Delta_d * (2 * Delta_d * X * pow(delta, 2) - 1) - Delta_d))
1039 / pow(2 * Delta * X * delta + 1, 4);
1040 CoolPropDbl dXdd_ddelta = pow(X, 2)
1041 * (-24 * pow(Delta, 4) * pow(X, 3) * delta - 8 * pow(Delta, 3) * Delta_d * pow(X, 3) * pow(delta, 2)
1042 - 18 * pow(Delta, 3) * pow(X, 2) + 8 * pow(Delta, 2) * Delta_d * pow(X, 2) * delta
1043 - 4 * pow(Delta, 2) * Delta_dd * pow(X, 2) * pow(delta, 2) + 10 * Delta * pow(Delta_d, 2) * pow(X, 2) * pow(delta, 2)
1044 + 12 * Delta * Delta_d * X - 4 * Delta * Delta_dd * X * delta + 8 * pow(Delta_d, 2) * X * delta - Delta_dd)
1045 / (16 * pow(Delta, 4) * pow(X, 4) * pow(delta, 4) + 32 * pow(Delta, 3) * pow(X, 3) * pow(delta, 3)
1046 + 24 * pow(Delta, 2) * pow(X, 2) * pow(delta, 2) + 8 * Delta * X * delta + 1);
1047 CoolPropDbl dXdd_dDelta = pow(X, 3)
1048 * (-8 * pow(Delta, 2) * Delta_d * pow(X, 2) * pow(delta, 3) + 8 * pow(Delta, 2) * Delta_dd * pow(X, 2) * pow(delta, 4)
1049 + 10 * pow(Delta, 2) * X * delta - 24 * Delta * pow(Delta_d, 2) * pow(X, 2) * pow(delta, 4)
1050 + 8 * Delta * Delta_d * X * pow(delta, 2) + 8 * Delta * Delta_dd * X * pow(delta, 3) + 8 * Delta
1051 - 18 * pow(Delta_d, 2) * X * pow(delta, 3) + 12 * Delta_d * delta + 2 * Delta_dd * pow(delta, 2))
1052 / (16 * pow(Delta, 4) * pow(X, 4) * pow(delta, 4) + 32 * pow(Delta, 3) * pow(X, 3) * pow(delta, 3)
1053 + 24 * pow(Delta, 2) * pow(X, 2) * pow(delta, 2) + 8 * Delta * X * delta + 1);
1054 CoolPropDbl dXdd_dDelta_d =
1055 2 * pow(X, 2)
1056 * (2 * X * delta * (Delta + Delta_d * delta) * (Delta * X * delta + 1) + (2 * Delta * X * delta + 1) * (2 * Delta_d * X * pow(delta, 2) - 1))
1057 / pow(2 * Delta * X * delta + 1, 3);
1058 CoolPropDbl dXdd_dDelta_dd = -pow(X, 2) * delta / (2 * Delta * X * delta + 1);
1059
1060 return dXdd_ddelta + dXdd_dX * (dX_ddelta + dX_dDelta * Delta_d) + dXdd_dDelta * Delta_d + dXdd_dDelta_d * Delta_dd + dXdd_dDelta_dd * Delta_ddd;
1061}
1062CoolPropDbl ResidualHelmholtzSAFTAssociating::g(const CoolPropDbl& eta) const {
1063 return 0.5 * (2 - eta) / pow(1 - eta, (int)3);
1064}
1065CoolPropDbl ResidualHelmholtzSAFTAssociating::dg_deta(const CoolPropDbl& eta) const {
1066 return 0.5 * (5 - 2 * eta) / pow(1 - eta, (int)4);
1067}
1068CoolPropDbl ResidualHelmholtzSAFTAssociating::d2g_deta2(const CoolPropDbl& eta) const {
1069 return 3 * (3 - eta) / pow(1 - eta, (int)5);
1070}
1071CoolPropDbl ResidualHelmholtzSAFTAssociating::d3g_deta3(const CoolPropDbl& eta) const {
1072 return 6 * (7 - 2 * eta) / pow(1 - eta, (int)6);
1073}
1074CoolPropDbl ResidualHelmholtzSAFTAssociating::eta(const CoolPropDbl& delta) const {
1075 return this->vbarn * delta;
1076}
1077
1078void ResidualHelmholtzSAFTAssociating::all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& deriv) throw() {
1079 if (disabled) {
1080 return;
1081 }
1082 CoolPropDbl X = this->X(delta, this->Deltabar(tau, delta));
1083 CoolPropDbl X_t = this->dX_dtau(tau, delta);
1084 CoolPropDbl X_d = this->dX_ddelta(tau, delta);
1085 CoolPropDbl X_tt = this->d2X_dtau2(tau, delta);
1086 CoolPropDbl X_dd = this->d2X_ddelta2(tau, delta);
1087 CoolPropDbl X_dt = this->d2X_ddeltadtau(tau, delta);
1088 CoolPropDbl X_ttt = this->d3X_dtau3(tau, delta);
1089 CoolPropDbl X_dtt = this->d3X_ddeltadtau2(tau, delta);
1090 CoolPropDbl X_ddt = this->d3X_ddelta2dtau(tau, delta);
1091 CoolPropDbl X_ddd = this->d3X_ddelta3(tau, delta);
1092
1093 deriv.alphar += this->m * this->a * ((log(X) - X / 2.0 + 0.5));
1094 deriv.dalphar_ddelta += this->m * this->a * (1 / X - 0.5) * this->dX_ddelta(tau, delta);
1095 deriv.dalphar_dtau += this->m * this->a * (1 / X - 0.5) * this->dX_dtau(tau, delta);
1096 deriv.d2alphar_dtau2 += this->m * this->a * ((1 / X - 0.5) * X_tt - pow(X_t / X, 2));
1097 deriv.d2alphar_ddelta2 += this->m * this->a * ((1 / X - 0.5) * X_dd - pow(X_d / X, 2));
1098 deriv.d2alphar_ddelta_dtau += this->m * this->a * ((-X_t / X / X) * X_d + X_dt * (1 / X - 0.5));
1099 deriv.d3alphar_dtau3 += this->m * this->a
1100 * ((1 / X - 1.0 / 2.0) * X_ttt + (-X_t / pow(X, (int)2)) * X_tt
1101 - 2 * (pow(X, (int)2) * (X_t * X_tt) - pow(X_t, (int)2) * (X * X_t)) / pow(X, (int)4));
1102 deriv.d3alphar_ddelta_dtau2 += this->m * this->a
1103 * ((1 / X - 1.0 / 2.0) * X_dtt - X_d / pow(X, (int)2) * X_tt
1104 - 2 * (pow(X, (int)2) * (X_t * X_dt) - pow(X_t, (int)2) * (X * X_d)) / pow(X, (int)4));
1105 deriv.d3alphar_ddelta2_dtau += this->m * this->a
1106 * ((1 / X - 1.0 / 2.0) * X_ddt - X_t / pow(X, (int)2) * X_dd
1107 - 2 * (pow(X, (int)2) * (X_d * X_dt) - pow(X_d, (int)2) * (X * X_t)) / pow(X, (int)4));
1108 deriv.d3alphar_ddelta3 += this->m * this->a
1109 * ((1 / X - 1.0 / 2.0) * X_ddd - X_d / pow(X, (int)2) * X_dd
1110 - 2 * (pow(X, (int)2) * (X_d * X_dd) - pow(X_d, (int)2) * (X * X_d)) / pow(X, (int)4));
1111}
1112
1113#if ENABLE_CATCH
1114mcx::MultiComplex<double> IdealHelmholtzCP0PolyT::one_mcx(const mcx::MultiComplex<double>& tau,
1115 const mcx::MultiComplex<double>& delta) const {
1116 mcx::MultiComplex<double> sum = 0.0;
1117 for (std::size_t i = 0; i < N; ++i) {
1118 if (std::abs(t[i]) < 10 * DBL_EPSILON) {
1119 sum += c[i] - c[i] * tau / tau0 + c[i] * log(tau / tau0);
1120 } else if (std::abs(t[i] + 1) < 10 * DBL_EPSILON) {
1121 sum += c[i] * tau / Tc * log(tau0 / tau) + c[i] / Tc * (tau - tau0);
1122 } else {
1123 sum += -c[i] * pow(Tc, t[i]) * pow(tau, -t[i]) / (t[i] * (t[i] + 1)) - c[i] * pow(T0, t[i] + 1) * tau / (Tc * (t[i] + 1))
1124 + c[i] * pow(T0, t[i]) / t[i];
1125 }
1126 }
1127 return sum;
1128}
1129#endif
1130
1131void IdealHelmholtzCP0PolyT::to_json(rapidjson::Value& el, rapidjson::Document& doc) {
1132 el.AddMember("type", "IdealGasCP0Poly", doc.GetAllocator());
1133
1134 rapidjson::Value _c(rapidjson::kArrayType), _t(rapidjson::kArrayType);
1135 for (std::size_t i = 0; i < N; ++i) {
1136 _c.PushBack(static_cast<double>(c[i]), doc.GetAllocator());
1137 _t.PushBack(static_cast<double>(t[i]), doc.GetAllocator());
1138 }
1139 el.AddMember("c", _c, doc.GetAllocator());
1140 el.AddMember("t", _t, doc.GetAllocator());
1141 el.AddMember("Tc", static_cast<double>(Tc), doc.GetAllocator());
1142 el.AddMember("T0", static_cast<double>(T0), doc.GetAllocator());
1143}
1144
1145void IdealHelmholtzLead::all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw() {
1146 if (!enabled) {
1147 return;
1148 }
1149 derivs.alphar += log(delta) + a1 + a2 * tau;
1150 derivs.dalphar_ddelta += 1.0 / delta;
1151 derivs.dalphar_dtau += a2;
1152 derivs.d2alphar_ddelta2 += -1.0 / delta / delta;
1153 derivs.d3alphar_ddelta3 += 2 / delta / delta / delta;
1154 derivs.d4alphar_ddelta4 += -6 / POW4(delta);
1155}
1156void IdealHelmholtzEnthalpyEntropyOffset::all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw() {
1157 if (!enabled) {
1158 return;
1159 }
1160 derivs.alphar += a1 + a2 * tau;
1161 derivs.dalphar_dtau += a2;
1162}
1163void IdealHelmholtzLogTau::all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw() {
1164
1165 if (!enabled) {
1166 return;
1167 }
1168 derivs.alphar += a1 * log(tau);
1169 derivs.dalphar_dtau += a1 / tau;
1170 derivs.d2alphar_dtau2 += -a1 / tau / tau;
1171 derivs.d3alphar_dtau3 += 2 * a1 / tau / tau / tau;
1172 derivs.d4alphar_dtau4 += -6 * a1 / POW4(tau);
1173}
1174void IdealHelmholtzPower::all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw() {
1175 if (!enabled) {
1176 return;
1177 }
1178 {
1179 CoolPropDbl s = 0;
1180 for (std::size_t i = 0; i < N; ++i) {
1181 s += n[i] * pow(tau, t[i]);
1182 }
1183 derivs.alphar += s;
1184 }
1185 {
1186 CoolPropDbl s = 0;
1187 for (std::size_t i = 0; i < N; ++i) {
1188 s += n[i] * t[i] * pow(tau, t[i] - 1);
1189 }
1190 derivs.dalphar_dtau += s;
1191 }
1192 {
1193 CoolPropDbl s = 0;
1194 for (std::size_t i = 0; i < N; ++i) {
1195 s += n[i] * t[i] * (t[i] - 1) * pow(tau, t[i] - 2);
1196 }
1197 derivs.d2alphar_dtau2 += s;
1198 }
1199 {
1200 CoolPropDbl s = 0;
1201 for (std::size_t i = 0; i < N; ++i) {
1202 s += n[i] * t[i] * (t[i] - 1) * (t[i] - 2) * pow(tau, t[i] - 3);
1203 }
1204 derivs.d3alphar_dtau3 += s;
1205 }
1206 {
1207 CoolPropDbl s = 0;
1208 for (std::size_t i = 0; i < N; ++i) {
1209 s += n[i] * t[i] * (t[i] - 1) * (t[i] - 2) * (t[i] - 3) * pow(tau, t[i] - 4);
1210 }
1211 derivs.d4alphar_dtau4 += s;
1212 }
1213}
1214void IdealHelmholtzPlanckEinsteinGeneralized::all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw() {
1215 // First pre-calculate exp(theta[i]*tau) for each contribution; used in each term
1216 std::vector<double> expthetatau(N);
1217 for (std::size_t i = 0; i < N; ++i) {
1218 expthetatau[i] = exp(theta[i] * tau);
1219 }
1220
1221 if (!enabled) {
1222 return;
1223 }
1224 {
1225 CoolPropDbl s = 0;
1226 for (std::size_t i = 0; i < N; ++i) {
1227 s += n[i] * log(c[i] + d[i] * expthetatau[i]);
1228 }
1229 derivs.alphar += s;
1230 }
1231 {
1232 CoolPropDbl s = 0;
1233 for (std::size_t i = 0; i < N; ++i) {
1234 s += n[i] * theta[i] * d[i] * expthetatau[i] / (c[i] + d[i] * expthetatau[i]);
1235 }
1236 derivs.dalphar_dtau += s;
1237 }
1238 {
1239 CoolPropDbl s = 0;
1240 for (std::size_t i = 0; i < N; ++i) {
1241 s += n[i] * POW2(theta[i]) * c[i] * d[i] * expthetatau[i] / pow(c[i] + d[i] * expthetatau[i], 2);
1242 }
1243 derivs.d2alphar_dtau2 += s;
1244 }
1245 {
1246 CoolPropDbl s = 0;
1247 for (std::size_t i = 0; i < N; ++i) {
1248 s += n[i] * POW3(theta[i]) * c[i] * d[i] * (c[i] - d[i] * expthetatau[i]) * expthetatau[i] / pow(c[i] + d[i] * expthetatau[i], 3);
1249 }
1250 derivs.d3alphar_dtau3 += s;
1251 }
1252 {
1253 CoolPropDbl s = 0;
1254 for (std::size_t i = 0; i < N; ++i) {
1255 const CoolPropDbl para = c[i] + d[i] * expthetatau[i];
1256 const CoolPropDbl bracket = 6 * POW3(d[i]) * POW3(expthetatau[i]) - 12 * d[i] * d[i] * para * POW2(expthetatau[i])
1257 + 7 * d[i] * POW2(para) * expthetatau[i] - POW3(para);
1258 s += -n[i] * d[i] * pow(theta[i], 4) * bracket * expthetatau[i] / pow(c[i] + d[i] * expthetatau[i], 4);
1259 }
1260 derivs.d4alphar_dtau4 += s;
1261 }
1262}
1263void IdealHelmholtzCP0Constant::all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) throw() {
1264 if (!enabled) {
1265 return;
1266 }
1267 derivs.alphar += cp_over_R - cp_over_R * tau / tau0 + cp_over_R * log(tau / tau0);
1268 derivs.dalphar_dtau += cp_over_R / tau - cp_over_R / tau0;
1269 derivs.d2alphar_dtau2 += -cp_over_R / (tau * tau);
1270 derivs.d3alphar_dtau3 += 2 * cp_over_R / (tau * tau * tau);
1271 derivs.d4alphar_dtau4 += -6 * cp_over_R / POW4(tau);
1272}
1273void IdealHelmholtzCP0PolyT::all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) {
1274 if (!enabled) {
1275 return;
1276 }
1277 {
1278 double sum = 0;
1279 for (std::size_t i = 0; i < N; ++i) {
1280 if (std::abs(t[i]) < 10 * DBL_EPSILON) {
1281 sum += c[i] - c[i] * tau / tau0 + c[i] * log(tau / tau0);
1282 } else if (std::abs(t[i] + 1) < 10 * DBL_EPSILON) {
1283 sum += c[i] * tau / Tc * log(tau0 / tau) + c[i] / Tc * (tau - tau0);
1284 } else {
1285 sum += -c[i] * pow(Tc, t[i]) * pow(tau, -t[i]) / (t[i] * (t[i] + 1)) - c[i] * pow(T0, t[i] + 1) * tau / (Tc * (t[i] + 1))
1286 + c[i] * pow(T0, t[i]) / t[i];
1287 }
1288 }
1289 derivs.alphar += sum;
1290 }
1291 {
1292 double sum = 0;
1293 for (std::size_t i = 0; i < N; ++i) {
1294 if (std::abs(t[i]) < 10 * DBL_EPSILON) {
1295 sum += c[i] / tau - c[i] / tau0;
1296 } else if (std::abs(t[i] + 1) < 10 * DBL_EPSILON) {
1297 sum += c[i] / Tc * log(tau0 / tau);
1298 } else {
1299 sum += c[i] * pow(Tc, t[i]) * pow(tau, -t[i] - 1) / (t[i] + 1) - c[i] * pow(Tc, t[i]) / (pow(tau0, t[i] + 1) * (t[i] + 1));
1300 }
1301 }
1302 derivs.dalphar_dtau += sum;
1303 }
1304 {
1305 double sum = 0;
1306 for (std::size_t i = 0; i < N; ++i) {
1307 if (std::abs(t[i]) < 10 * DBL_EPSILON) {
1308 sum += -c[i] / (tau * tau);
1309 } else if (std::abs(t[i] + 1) < 10 * DBL_EPSILON) {
1310 sum += -c[i] / (tau * Tc);
1311 } else {
1312 sum += -c[i] * pow(Tc / tau, t[i]) / (tau * tau);
1313 }
1314 }
1315 derivs.d2alphar_dtau2 += sum;
1316 }
1317 {
1318 double sum = 0;
1319 for (std::size_t i = 0; i < N; ++i) {
1320 if (std::abs(t[i]) < 10 * DBL_EPSILON) {
1321 sum += 2 * c[i] / (tau * tau * tau);
1322 } else if (std::abs(t[i] + 1) < 10 * DBL_EPSILON) {
1323 sum += c[i] / (tau * tau * Tc);
1324 } else {
1325 sum += c[i] * pow(Tc / tau, t[i]) * (t[i] + 2) / (tau * tau * tau);
1326 }
1327 }
1328 derivs.d3alphar_dtau3 += sum;
1329 }
1330 {
1331 double sum = 0;
1332 for (std::size_t i = 0; i < N; ++i) {
1333 if (std::abs(t[i]) < 10 * DBL_EPSILON) {
1334 sum += -6 * c[i] / POW4(tau);
1335 } else if (std::abs(t[i] + 1) < 10 * DBL_EPSILON) {
1336 sum += -3 * c[i] / (POW3(tau) * Tc);
1337 } else {
1338 sum += -c[i] * (t[i] + 2) * (t[i] + 3) * pow(Tc / tau, t[i]) / (tau * tau * tau * tau);
1339 }
1340 }
1341 derivs.d4alphar_dtau4 += sum;
1342 }
1343}
1344
1345void IdealHelmholtzGERG2004Sinh::all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) {
1346 if (!enabled) {
1347 return;
1348 }
1349 // Check that the reducing temperature value is provided
1350 CoolPropDbl T_red = _HUGE;
1351 if (ValidNumber(_Tr)) {
1352 T_red = _Tr; // Primarily useful for testing
1353 } else if (ValidNumber(derivs.T_red)) {
1354 T_red = derivs.T_red;
1355 } else {
1356 throw ValueError("T_red needs to be stored somewhere for GERG2004Sinh");
1357 }
1358 CoolPropDbl Tci_over_Tr = Tc / T_red;
1359
1360 double sum00 = 0, sum10 = 0, sum20 = 0, sum30 = 0, sum40 = 0;
1361 for (std::size_t i = 0; i < N; ++i) {
1362 CoolPropDbl t = theta[i] * Tci_over_Tr;
1363 sum00 += n[i] * log(std::abs(sinh(t * tau)));
1364 sum10 += n[i] * t / tanh(t * tau);
1365 sum20 += -n[i] * POW2(t) / POW2(sinh(t * tau));
1366 sum30 += -2 * n[i] * POW3(t) * (1 / tanh(t * tau) - 1 / POW3(tanh(t * tau)));
1367 sum40 += -2 * n[i] * POW4(t) * (1 - 4 / POW2(tanh(t * tau)) + 3 / POW4(tanh(t * tau)));
1368 }
1369 derivs.alphar += sum00;
1370 derivs.dalphar_dtau += sum10;
1371 derivs.d2alphar_dtau2 += sum20;
1372 derivs.d3alphar_dtau3 += sum30;
1373 derivs.d4alphar_dtau4 += sum40;
1374}
1375
1376#if ENABLE_CATCH
1377mcx::MultiComplex<double> IdealHelmholtzGERG2004Sinh::one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const {
1378 // Check that the reducing temperature value is provided in the derivs structure
1379 CoolPropDbl T_red = _HUGE;
1380 if (ValidNumber(_Tr)) {
1381 T_red = _Tr; // Primarily useful for testing
1382 } else {
1383 throw ValueError("T_red needs to be stored somewhere for GERG2004Cosh");
1384 }
1385 auto Tci_over_Tr = Tc / T_red;
1386 auto diff_abs = [](const mcx::MultiComplex<double>& x) { return sqrt(x * x); }; // differentiable absolute value
1387
1388 mcx::MultiComplex<double> sum00 = 0.0;
1389 for (std::size_t i = 0; i < N; ++i) {
1390 auto t = theta[i] * Tci_over_Tr;
1391 sum00 += n[i] * log(diff_abs(sinh(t * tau)));
1392 }
1393 return sum00;
1394}
1395#endif
1396
1397void IdealHelmholtzGERG2004Cosh::all(const CoolPropDbl& tau, const CoolPropDbl& delta, HelmholtzDerivatives& derivs) {
1398 if (!enabled) {
1399 return;
1400 }
1401 // Check that the reducing temperature value is provided in the derivs structure
1402 CoolPropDbl T_red = _HUGE;
1403 if (ValidNumber(_Tr)) {
1404 T_red = _Tr; // Primarily useful for testing
1405 } else if (ValidNumber(derivs.T_red)) {
1406 T_red = derivs.T_red;
1407 } else {
1408 throw ValueError("T_red needs to be stored somewhere for GERG2004Cosh");
1409 }
1410 CoolPropDbl Tci_over_Tr = Tc / T_red;
1411
1412 double sum00 = 0, sum10 = 0, sum20 = 0, sum30 = 0, sum40 = 0;
1413 for (std::size_t i = 0; i < N; ++i) {
1414 CoolPropDbl t = theta[i] * Tci_over_Tr;
1415 sum00 += -n[i] * log(std::abs(cosh(t * tau)));
1416 sum10 += -n[i] * t * tanh(t * tau);
1417 sum20 += -n[i] * POW2(t) / POW2(cosh(t * tau));
1418 sum30 += -2 * n[i] * POW3(t) * (POW3(tanh(t * tau)) - tanh(t * tau));
1419 sum40 += 2 * n[i] * POW4(t) * (3 * POW4(tanh(t * tau)) - 4 * POW2(tanh(t * tau)) + 1);
1420 }
1421 derivs.alphar += sum00;
1422 derivs.dalphar_dtau += sum10;
1423 derivs.d2alphar_dtau2 += sum20;
1424 derivs.d3alphar_dtau3 += sum30;
1425 derivs.d4alphar_dtau4 += sum40;
1426}
1427
1428#if ENABLE_CATCH
1429mcx::MultiComplex<double> IdealHelmholtzGERG2004Cosh::one_mcx(const mcx::MultiComplex<double>& tau, const mcx::MultiComplex<double>& delta) const {
1430 // Check that the reducing temperature value is provided in the derivs structure
1431 CoolPropDbl T_red = _HUGE;
1432 if (ValidNumber(_Tr)) {
1433 T_red = _Tr; // Primarily useful for testing
1434 } else {
1435 throw ValueError("T_red needs to be stored somewhere for GERG2004Cosh");
1436 }
1437 auto Tci_over_Tr = Tc / T_red;
1438 auto diff_abs = [](const mcx::MultiComplex<double>& x) { return sqrt(x * x); }; // differentiable absolute value
1439
1440 mcx::MultiComplex<double> sum00 = 0.0;
1441 for (std::size_t i = 0; i < N; ++i) {
1442 auto t = theta[i] * Tci_over_Tr;
1443 sum00 += -n[i] * log(diff_abs(cosh(t * tau)));
1444 }
1445 return sum00;
1446}
1447#endif
1448
1449
1450//void IdealHelmholtzCP0AlyLee::to_json(rapidjson::Value &el, rapidjson::Document &doc){
1451// el.AddMember("type","IdealGasHelmholtzCP0AlyLee",doc.GetAllocator());
1452// rapidjson::Value _n(rapidjson::kArrayType);
1453// for (std::size_t i=0; i<=4; ++i)
1454// {
1455// _n.PushBack(static_cast<double>(c[i]),doc.GetAllocator());
1456// }
1457// el.AddMember("c",_n,doc.GetAllocator());
1458// el.AddMember("Tc",static_cast<double>(Tc),doc.GetAllocator());
1459// el.AddMember("T0",static_cast<double>(T0),doc.GetAllocator());
1460//}
1461//CoolPropDbl IdealHelmholtzCP0AlyLee::base(const CoolPropDbl &tau, const CoolPropDbl &delta) throw()
1462//{
1463// if (!enabled){ return 0.0;}
1464// return -tau*(anti_deriv_cp0_tau2(tau)-anti_deriv_cp0_tau2(tau0))+(anti_deriv_cp0_tau(tau)-anti_deriv_cp0_tau(tau0));
1465//}
1466//CoolPropDbl IdealHelmholtzCP0AlyLee::dTau(const CoolPropDbl &tau, const CoolPropDbl &delta) throw()
1467//{
1468// if (!enabled){ return 0.0;}
1469// return -(anti_deriv_cp0_tau2(tau) - anti_deriv_cp0_tau2(tau0));
1470//}
1471//CoolPropDbl IdealHelmholtzCP0AlyLee::anti_deriv_cp0_tau2(const CoolPropDbl &tau)
1472//{
1473// return -c[0]/tau + 2*c[1]*c[2]/Tc/(exp(-2*c[2]*tau/Tc)-1) - 2*c[3]*c[4]/Tc*(exp(2*c[4]*tau/Tc)+1);
1474//}
1475//CoolPropDbl IdealHelmholtzCP0AlyLee::anti_deriv_cp0_tau(const CoolPropDbl &tau)
1476//{
1477// CoolPropDbl term1 = c[0]*log(tau);
1478// CoolPropDbl term2 = 2*c[1]*c[2]*tau/(-Tc + Tc*exp(-2*c[2]*tau/Tc)) + c[1]*log(-1 + exp(-2*c[2]*tau/Tc)) + 2*c[1]*c[2]*tau/Tc;
1479// CoolPropDbl term3 = -c[3]*(Tc*exp(2*c[4]*tau/Tc)*log(exp(2*c[4]*tau/Tc) + 1) + Tc*log(exp(2*c[4]*tau/Tc) + 1) - 2*c[4]*tau*exp(2*c[4]*tau/Tc))/(Tc*(exp(2*c[4]*tau/Tc) + 1));
1480// return term1 + term2 + term3;
1481//}
1482//CoolPropDbl IdealHelmholtzCP0AlyLee::dTau2(const CoolPropDbl &tau, const CoolPropDbl &delta) throw()
1483//{
1484// if (!enabled){ return 0.0;}
1485// return -c[0]/pow(tau,2) - c[1]*pow(c[2]/Tc/sinh(c[2]*tau/Tc),2) - c[3]*pow(c[4]/Tc/cosh(c[4]*tau/Tc),2);
1486//}
1487//CoolPropDbl IdealHelmholtzCP0AlyLee::dTau3(const CoolPropDbl &tau, const CoolPropDbl &delta) throw()
1488//{
1489// if (!enabled){ return 0.0;}
1490// return 2*c[0]/pow(tau,3) + 2*c[1]*pow(c[2]/Tc,3)*cosh(c[2]*tau/Tc)/pow(sinh(c[2]*tau/Tc),3) + 2*c[3]*pow(c[4]/Tc,3)*sinh(c[4]*tau/Tc)/pow(cosh(c[4]*tau/Tc),3);
1491//}
1492
1493}; /* namespace CoolProp */
1494
1495/*
1496IdealHelmholtzEnthalpyEntropyOffset EnthalpyEntropyOffset;
1497*/
1498
1499#ifdef ENABLE_CATCH
1500# include <math.h>
1501# include <catch2/catch_all.hpp>
1502# include "crossplatform_shared_ptr.h"
1503
1504class HelmholtzConsistencyFixture
1505{
1506 public:
1507 CoolPropDbl numerical, analytic;
1508
1509 shared_ptr<CoolProp::BaseHelmholtzTerm> PlanckEinstein, Lead, LogTau, IGPower, CP0Constant, CP0PolyT, SAFT, NonAnalytic, Soave, PR, XiangDeiters,
1510 GaoB, GERG2004Cosh, GERG2004Sinh;
1511 shared_ptr<CoolProp::ResidualHelmholtzGeneralizedExponential> Gaussian, Lemmon2005, Exponential, GERG2008, Power;
1512
1513 HelmholtzConsistencyFixture() {
1514 shared_ptr<AbstractCubic> _SRK(new SRK(300, 4e6, 0.3, 8.314461));
1515 _SRK->set_Tr(300);
1516 _SRK->set_rhor(4000);
1517 Soave.reset(new CoolProp::ResidualHelmholtzGeneralizedCubic(_SRK));
1518
1519 shared_ptr<AbstractCubic> _PR(new PengRobinson(300, 4e6, 0.3, 8.314461));
1520 _PR->set_Tr(300);
1521 _PR->set_rhor(4000);
1522 PR.reset(new CoolProp::ResidualHelmholtzGeneralizedCubic(_PR));
1523
1524 {
1525 // Signs of eta are flipped relative to paper from Gao et al., implemented with opposite sign in CoolProp
1526 std::vector<CoolPropDbl> beta = {0.3696, 0.2962}, epsilon = {0.4478, 0.44689}, eta = {-2.8452, -2.8342}, gamma = {1.108, 1.313},
1527 n = {-1.6909858, 0.93739074}, t = {4.3315, 4.015}, d = {1, 1}, b = {1.244, 0.6826};
1528 GaoB.reset(new CoolProp::ResidualHelmholtzGaoB(n, t, d, eta, beta, gamma, epsilon, b));
1529 }
1530
1531 XiangDeiters.reset(new CoolProp::ResidualHelmholtzXiangDeiters(300, 4e6, 4000, 0.3, 8.3144621));
1532
1533 Lead.reset(new CoolProp::IdealHelmholtzLead(1, 3));
1534 LogTau.reset(new CoolProp::IdealHelmholtzLogTau(1.5));
1535 {
1536 std::vector<CoolPropDbl> n(4, 0), t(4, 1);
1537 n[0] = -0.1;
1538 n[2] = 0.1;
1539 t[1] = -1;
1540 t[2] = -2;
1541 t[3] = 2;
1542 IGPower.reset(new CoolProp::IdealHelmholtzPower(n, t));
1543 }
1544 {
1545 std::vector<CoolPropDbl> n(4, 0), t(4, 1), c(4, 1), d(4, -1);
1546 n[0] = 0.1;
1547 n[2] = 0.5;
1548 t[0] = -1.5;
1549 t[1] = -1;
1550 t[2] = -2;
1551 t[3] = -2;
1552 PlanckEinstein.reset(new CoolProp::IdealHelmholtzPlanckEinsteinGeneralized(n, t, c, d));
1553 }
1554 {
1555 std::vector<CoolPropDbl> c(3, 1), t(3, 0);
1556 t[1] = 1;
1557 t[2] = 2;
1558 c[1] = 2;
1559 c[2] = 3;
1560 CoolPropDbl T0 = 273.15, Tc = 345.857;
1561 CP0PolyT.reset(new CoolProp::IdealHelmholtzCP0PolyT(c, t, Tc, T0));
1562 }
1563 CP0Constant.reset(new CoolProp::IdealHelmholtzCP0Constant(4 / 8.314472, 300, 250));
1564 {
1565 // Nitrogen
1566 std::vector<CoolPropDbl> n(2, 0.0);
1567 n[0] = 0.137320000;
1568 n[1] = 0.900660000;
1569 std::vector<CoolPropDbl> theta(2, 0.0);
1570 theta[0] = 5.251822620;
1571 theta[1] = 13.788988208;
1572 CoolPropDbl rhomolar_crit = 11183.900000, T_crit = 126.192000000;
1573 GERG2004Cosh.reset(new CoolProp::IdealHelmholtzGERG2004Cosh(n, theta, T_crit));
1574 static_cast<CoolProp::IdealHelmholtzGERG2004Cosh*>(GERG2004Cosh.get())->set_Tred(T_crit * 1.3);
1575 }
1576 {
1577 // Nitrogen
1578 std::vector<CoolPropDbl> n(1, 0.0);
1579 n[0] = -0.146600000;
1580 std::vector<CoolPropDbl> theta(1, 0.0);
1581 theta[0] = -5.393067706;
1582 CoolPropDbl rhomolar_crit = 11183.900000, T_crit = 126.192000000;
1583 GERG2004Sinh.reset(new CoolProp::IdealHelmholtzGERG2004Sinh(n, theta, T_crit));
1584 static_cast<CoolProp::IdealHelmholtzGERG2004Sinh*>(GERG2004Sinh.get())->set_Tred(T_crit * 1.3);
1585 }
1586
1587 {
1588 CoolPropDbl beta[] = {1.24, 0.821, 15.45, 2.21, 437, 0.743}, d[] = {1, 1, 2, 2, 3, 3},
1589 epsilon[] = {0.6734, 0.9239, 0.8636, 1.0507, 0.8482, 0.7522}, eta[] = {0.9667, 1.5154, 1.0591, 1.6642, 12.4856, 0.9662},
1590 gamma[] = {1.2827, 0.4317, 1.1217, 1.1871, 1.1243, 0.4203},
1591 n[] = {1.2198, -0.4883, -0.0033293, -0.0035387, -0.51172, -0.16882}, t[] = {1, 2.124, 0.4, 3.5, 0.5, 2.7};
1592 Gaussian.reset(new CoolProp::ResidualHelmholtzGeneralizedExponential());
1593 Gaussian->add_Gaussian(
1594 std::vector<CoolPropDbl>(n, n + sizeof(n) / sizeof(n[0])), std::vector<CoolPropDbl>(d, d + sizeof(d) / sizeof(d[0])),
1595 std::vector<CoolPropDbl>(t, t + sizeof(t) / sizeof(t[0])), std::vector<CoolPropDbl>(eta, eta + sizeof(eta) / sizeof(eta[0])),
1596 std::vector<CoolPropDbl>(epsilon, epsilon + sizeof(epsilon) / sizeof(epsilon[0])),
1597 std::vector<CoolPropDbl>(beta, beta + sizeof(beta) / sizeof(beta[0])),
1598 std::vector<CoolPropDbl>(gamma, gamma + sizeof(gamma) / sizeof(gamma[0])));
1599 }
1600 {
1601 CoolPropDbl d[] = {1, 1, 1, 2, 4, 1, 1, 2, 2, 3, 4, 5, 1, 5, 1, 2, 3, 5}, l[] = {0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 2, 3, 3},
1602 m[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.7, 7, 6},
1603 n[] = {5.28076, -8.67658, 0.7501127, 0.7590023, 0.01451899, 4.777189, -3.330988, 3.775673, -2.290919,
1604 0.8888268, -0.6234864, -0.04127263, -0.08455389, -0.1308752, 0.008344962, -1.532005, -0.05883649, 0.02296658},
1605 t[] = {0.669, 1.05, 2.75, 0.956, 1, 2, 2.75, 2.38, 3.37, 3.47, 2.63, 3.45, 0.72, 4.23, 0.2, 4.5, 29, 24};
1606 Lemmon2005.reset(new CoolProp::ResidualHelmholtzGeneralizedExponential());
1607 Lemmon2005->add_Lemmon2005(
1608 std::vector<CoolPropDbl>(n, n + sizeof(n) / sizeof(n[0])), std::vector<CoolPropDbl>(d, d + sizeof(d) / sizeof(d[0])),
1609 std::vector<CoolPropDbl>(t, t + sizeof(t) / sizeof(t[0])), std::vector<CoolPropDbl>(l, l + sizeof(l) / sizeof(l[0])),
1610 std::vector<CoolPropDbl>(m, m + sizeof(m) / sizeof(m[0])));
1611 }
1612 {
1613 CoolPropDbl d[] = {1, 1, 1, 3, 7, 1, 2, 5, 1, 1, 4, 2}, l[] = {0, 0, 0, 0, 0, 1, 1, 1, 2, 2, 2, 3},
1614 n[] = {1.0038, -2.7662, 0.42921, 0.081363, 0.00024174, 0.48246, 0.75542, -0.00743, -0.4146, -0.016558, -0.10644, -0.021704},
1615 t[] = {0.25, 1.25, 1.5, 0.25, 0.875, 2.375, 2, 2.125, 3.5, 6.5, 4.75, 12.5};
1616 Power.reset(new CoolProp::ResidualHelmholtzGeneralizedExponential());
1617 Power->add_Power(std::vector<CoolPropDbl>(n, n + sizeof(n) / sizeof(n[0])), std::vector<CoolPropDbl>(d, d + sizeof(d) / sizeof(d[0])),
1618 std::vector<CoolPropDbl>(t, t + sizeof(t) / sizeof(t[0])), std::vector<CoolPropDbl>(l, l + sizeof(l) / sizeof(l[0])));
1619 }
1620 {
1621
1622 CoolPropDbl a = 1, epsilonbar = 12.2735737, kappabar = 1.09117041e-05, m = 1.01871348, vbarn = 0.0444215309;
1623 SAFT.reset(new CoolProp::ResidualHelmholtzSAFTAssociating(a, m, epsilonbar, vbarn, kappabar));
1624 }
1625 {
1626 CoolPropDbl n[] = {-0.666422765408, 0.726086323499, 0.0550686686128}, A[] = {0.7, 0.7, 0.7}, B[] = {0.3, 0.3, 1}, C[] = {10, 10, 12.5},
1627 D[] = {275, 275, 275}, a[] = {3.5, 3.5, 3}, b[] = {0.875, 0.925, 0.875}, beta[] = {0.3, 0.3, 0.3};
1628 NonAnalytic.reset(new CoolProp::ResidualHelmholtzNonAnalytic(
1629 std::vector<CoolPropDbl>(n, n + sizeof(n) / sizeof(n[0])), std::vector<CoolPropDbl>(a, a + sizeof(a) / sizeof(a[0])),
1630 std::vector<CoolPropDbl>(b, b + sizeof(b) / sizeof(b[0])), std::vector<CoolPropDbl>(beta, beta + sizeof(beta) / sizeof(beta[0])),
1631 std::vector<CoolPropDbl>(A, A + sizeof(A) / sizeof(A[0])), std::vector<CoolPropDbl>(B, B + sizeof(B) / sizeof(B[0])),
1632 std::vector<CoolPropDbl>(C, C + sizeof(C) / sizeof(C[0])), std::vector<CoolPropDbl>(D, D + sizeof(D) / sizeof(D[0]))));
1633 }
1634 {
1635 CoolPropDbl d[] = {2, 2, 2, 0, 0, 0}, g[] = {1.65533788, 1.65533788, 1.65533788, 1.65533788, 1.65533788, 1.65533788},
1636 l[] = {2, 2, 2, 2, 2, 2},
1637 n[] = {-3.821884669859, 8.30345065618981, -4.4832307260286, -1.02590136933231, 2.20786016506394, -1.07889905203761},
1638 t[] = {3, 4, 5, 3, 4, 5};
1639 Exponential.reset(new CoolProp::ResidualHelmholtzGeneralizedExponential());
1640 Exponential->add_Exponential(
1641 std::vector<CoolPropDbl>(n, n + sizeof(n) / sizeof(n[0])), std::vector<CoolPropDbl>(d, d + sizeof(d) / sizeof(n[0])),
1642 std::vector<CoolPropDbl>(t, t + sizeof(t) / sizeof(d[0])), std::vector<CoolPropDbl>(g, g + sizeof(g) / sizeof(t[0])),
1643 std::vector<CoolPropDbl>(l, l + sizeof(l) / sizeof(l[0])));
1644 }
1645 {
1646 CoolPropDbl d[] = {1, 4, 1, 2, 2, 2, 2, 2, 3}, t[] = {0.0, 1.85, 7.85, 5.4, 0.0, 0.75, 2.8, 4.45, 4.25},
1647 n[] = {-0.0098038985517335, 0.00042487270143005, -0.034800214576142, -0.13333813013896, -0.011993694974627,
1648 0.069243379775168, -0.31022508148249, 0.24495491753226, 0.22369816716981},
1649 eta[] = {0.0, 0.0, 1.0, 1.0, 0.25, 0.0, 0.0, 0.0, 0.0}, epsilon[] = {0.0, 0.0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5},
1650 beta[] = {0.0, 0.0, 1.0, 1.0, 2.5, 3.0, 3.0, 3.0, 3.0}, gamma[] = {0.0, 0.0, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5};
1651 GERG2008.reset(new CoolProp::ResidualHelmholtzGeneralizedExponential());
1652 GERG2008->add_GERG2008Gaussian(
1653 std::vector<CoolPropDbl>(n, n + sizeof(n) / sizeof(n[0])), std::vector<CoolPropDbl>(d, d + sizeof(d) / sizeof(n[0])),
1654 std::vector<CoolPropDbl>(t, t + sizeof(t) / sizeof(d[0])), std::vector<CoolPropDbl>(eta, eta + sizeof(eta) / sizeof(eta[0])),
1655 std::vector<CoolPropDbl>(epsilon, epsilon + sizeof(epsilon) / sizeof(epsilon[0])),
1656 std::vector<CoolPropDbl>(beta, beta + sizeof(beta) / sizeof(beta[0])),
1657 std::vector<CoolPropDbl>(gamma, gamma + sizeof(gamma) / sizeof(gamma[0])));
1658 }
1659 }
1660 void call(std::string d, shared_ptr<CoolProp::BaseHelmholtzTerm> term, CoolPropDbl tau, CoolPropDbl delta, CoolPropDbl ddelta) {
1661 if (!d.compare("dTau")) {
1662 return dTau(term, tau, delta, ddelta);
1663 } else if (!d.compare("dTau2")) {
1664 return dTau2(term, tau, delta, ddelta);
1665 } else if (!d.compare("dTau3")) {
1666 return dTau3(term, tau, delta, ddelta);
1667 } else if (!d.compare("dTau4")) {
1668 return dTau4(term, tau, delta, ddelta);
1669 } else if (!d.compare("dDelta")) {
1670 return dDelta(term, tau, delta, ddelta);
1671 } else if (!d.compare("dDelta2")) {
1672 return dDelta2(term, tau, delta, ddelta);
1673 } else if (!d.compare("dDelta3")) {
1674 return dDelta3(term, tau, delta, ddelta);
1675 } else if (!d.compare("dDelta4")) {
1676 return dDelta4(term, tau, delta, ddelta);
1677 } else if (!d.compare("dDelta_dTau")) {
1678 return dDelta_dTau(term, tau, delta, ddelta);
1679 } else if (!d.compare("dDelta_dTau2")) {
1680 return dDelta_dTau2(term, tau, delta, ddelta);
1681 } else if (!d.compare("dDelta2_dTau")) {
1682 return dDelta2_dTau(term, tau, delta, ddelta);
1683 } else if (!d.compare("dDelta3_dTau")) {
1684 return dDelta3_dTau(term, tau, delta, ddelta);
1685 } else if (!d.compare("dDelta2_dTau2")) {
1686 return dDelta2_dTau2(term, tau, delta, ddelta);
1687 } else if (!d.compare("dDelta_dTau3")) {
1688 return dDelta_dTau3(term, tau, delta, ddelta);
1689 } else {
1690 throw CoolProp::ValueError("don't understand deriv type");
1691 }
1692 }
1693 shared_ptr<CoolProp::BaseHelmholtzTerm> get(std::string t) {
1694 if (!t.compare("Lead")) {
1695 return Lead;
1696 } else if (!t.compare("LogTau")) {
1697 return LogTau;
1698 } else if (!t.compare("IGPower")) {
1699 return IGPower;
1700 } else if (!t.compare("PlanckEinstein")) {
1701 return PlanckEinstein;
1702 } else if (!t.compare("CP0Constant")) {
1703 return CP0Constant;
1704 } else if (!t.compare("CP0PolyT")) {
1705 return CP0PolyT;
1706 } else if (!t.compare("GERG2004Cosh")) {
1707 return GERG2004Cosh;
1708 } else if (!t.compare("GERG2004Sinh")) {
1709 return GERG2004Sinh;
1710 } else if (!t.compare("SRK")) {
1711 return Soave;
1712 } else if (!t.compare("PengRobinson")) {
1713 return PR;
1714 } else if (!t.compare("XiangDeiters")) {
1715 return XiangDeiters;
1716 } else if (!t.compare("GaoB")) {
1717 return GaoB;
1718 }
1719
1720 else if (!t.compare("Gaussian")) {
1721 return Gaussian;
1722 } else if (!t.compare("Lemmon2005")) {
1723 return Lemmon2005;
1724 } else if (!t.compare("Power")) {
1725 return Power;
1726 } else if (!t.compare("SAFT")) {
1727 return SAFT;
1728 } else if (!t.compare("NonAnalytic")) {
1729 return NonAnalytic;
1730 } else if (!t.compare("Exponential")) {
1731 return Exponential;
1732 } else if (!t.compare("GERG2008")) {
1733 return GERG2008;
1734 } else {
1735 throw CoolProp::ValueError(format("don't understand helmholtz type: %s", t.c_str()));
1736 }
1737 }
1738 void dTau(shared_ptr<CoolProp::BaseHelmholtzTerm> term, CoolPropDbl tau, CoolPropDbl delta, CoolPropDbl dtau) {
1739 CoolPropDbl term_plus = term->base(tau + dtau, delta);
1740 CoolPropDbl term_minus = term->base(tau - dtau, delta);
1741 numerical = (term_plus - term_minus) / (2 * dtau);
1742 analytic = term->dTau(tau, delta);
1743 };
1744 void dTau2(shared_ptr<CoolProp::BaseHelmholtzTerm> term, CoolPropDbl tau, CoolPropDbl delta, CoolPropDbl dtau) {
1745 CoolPropDbl term_plus = term->dTau(tau + dtau, delta);
1746 CoolPropDbl term_minus = term->dTau(tau - dtau, delta);
1747 numerical = (term_plus - term_minus) / (2 * dtau);
1748 analytic = term->dTau2(tau, delta);
1749 };
1750 void dTau3(shared_ptr<CoolProp::BaseHelmholtzTerm> term, CoolPropDbl tau, CoolPropDbl delta, CoolPropDbl dtau) {
1751 CoolPropDbl term_plus = term->dTau2(tau + dtau, delta);
1752 CoolPropDbl term_minus = term->dTau2(tau - dtau, delta);
1753 numerical = (term_plus - term_minus) / (2 * dtau);
1754 analytic = term->dTau3(tau, delta);
1755 };
1756 void dTau4(shared_ptr<CoolProp::BaseHelmholtzTerm> term, CoolPropDbl tau, CoolPropDbl delta, CoolPropDbl dtau) {
1757 CoolPropDbl term_plus = term->dTau3(tau + dtau, delta);
1758 CoolPropDbl term_minus = term->dTau3(tau - dtau, delta);
1759 numerical = (term_plus - term_minus) / (2 * dtau);
1760 analytic = term->dTau4(tau, delta);
1761 };
1762 void dDelta(shared_ptr<CoolProp::BaseHelmholtzTerm> term, CoolPropDbl tau, CoolPropDbl delta, CoolPropDbl ddelta) {
1763 CoolPropDbl term_plus = term->base(tau, delta + ddelta);
1764 CoolPropDbl term_minus = term->base(tau, delta - ddelta);
1765 numerical = (term_plus - term_minus) / (2 * ddelta);
1766 analytic = term->dDelta(tau, delta);
1767 };
1768 void dDelta2(shared_ptr<CoolProp::BaseHelmholtzTerm> term, CoolPropDbl tau, CoolPropDbl delta, CoolPropDbl ddelta) {
1769 CoolPropDbl term_plus = term->dDelta(tau, delta + ddelta);
1770 CoolPropDbl term_minus = term->dDelta(tau, delta - ddelta);
1771 numerical = (term_plus - term_minus) / (2 * ddelta);
1772 analytic = term->dDelta2(tau, delta);
1773 };
1774 void dDelta3(shared_ptr<CoolProp::BaseHelmholtzTerm> term, CoolPropDbl tau, CoolPropDbl delta, CoolPropDbl ddelta) {
1775 CoolPropDbl term_plus = term->dDelta2(tau, delta + ddelta);
1776 CoolPropDbl term_minus = term->dDelta2(tau, delta - ddelta);
1777 numerical = (term_plus - term_minus) / (2 * ddelta);
1778 analytic = term->dDelta3(tau, delta);
1779 };
1780 void dDelta4(shared_ptr<CoolProp::BaseHelmholtzTerm> term, CoolPropDbl tau, CoolPropDbl delta, CoolPropDbl ddelta) {
1781 CoolPropDbl term_plus = term->dDelta3(tau, delta + ddelta);
1782 CoolPropDbl term_minus = term->dDelta3(tau, delta - ddelta);
1783 numerical = (term_plus - term_minus) / (2 * ddelta);
1784 analytic = term->dDelta4(tau, delta);
1785 };
1786 void dDelta_dTau(shared_ptr<CoolProp::BaseHelmholtzTerm> term, CoolPropDbl tau, CoolPropDbl delta, CoolPropDbl ddelta) {
1787 CoolPropDbl term_plus = term->dTau(tau, delta + ddelta);
1788 CoolPropDbl term_minus = term->dTau(tau, delta - ddelta);
1789 numerical = (term_plus - term_minus) / (2 * ddelta);
1790 analytic = term->dDelta_dTau(tau, delta);
1791 };
1792 void dDelta_dTau2(shared_ptr<CoolProp::BaseHelmholtzTerm> term, CoolPropDbl tau, CoolPropDbl delta, CoolPropDbl ddelta) {
1793 CoolPropDbl term_plus = term->dTau2(tau, delta + ddelta);
1794 CoolPropDbl term_minus = term->dTau2(tau, delta - ddelta);
1795 numerical = (term_plus - term_minus) / (2 * ddelta);
1796 analytic = term->dDelta_dTau2(tau, delta);
1797 };
1798 void dDelta2_dTau(shared_ptr<CoolProp::BaseHelmholtzTerm> term, CoolPropDbl tau, CoolPropDbl delta, CoolPropDbl ddelta) {
1799 CoolPropDbl term_plus = term->dDelta_dTau(tau, delta + ddelta);
1800 CoolPropDbl term_minus = term->dDelta_dTau(tau, delta - ddelta);
1801 numerical = (term_plus - term_minus) / (2 * ddelta);
1802 analytic = term->dDelta2_dTau(tau, delta);
1803 };
1804 void dDelta3_dTau(shared_ptr<CoolProp::BaseHelmholtzTerm> term, CoolPropDbl tau, CoolPropDbl delta, CoolPropDbl ddelta) {
1805 CoolPropDbl term_plus = term->dDelta2_dTau(tau, delta + ddelta);
1806 CoolPropDbl term_minus = term->dDelta2_dTau(tau, delta - ddelta);
1807 numerical = (term_plus - term_minus) / (2 * ddelta);
1808 analytic = term->dDelta3_dTau(tau, delta);
1809 };
1810 void dDelta2_dTau2(shared_ptr<CoolProp::BaseHelmholtzTerm> term, CoolPropDbl tau, CoolPropDbl delta, CoolPropDbl ddelta) {
1811 CoolPropDbl term_plus = term->dDelta_dTau2(tau, delta + ddelta);
1812 CoolPropDbl term_minus = term->dDelta_dTau2(tau, delta - ddelta);
1813 numerical = (term_plus - term_minus) / (2 * ddelta);
1814 analytic = term->dDelta2_dTau2(tau, delta);
1815 };
1816 void dDelta_dTau3(shared_ptr<CoolProp::BaseHelmholtzTerm> term, CoolPropDbl tau, CoolPropDbl delta, CoolPropDbl ddelta) {
1817 CoolPropDbl term_plus = term->dTau3(tau, delta + ddelta);
1818 CoolPropDbl term_minus = term->dTau3(tau, delta - ddelta);
1819 numerical = (term_plus - term_minus) / (2 * ddelta);
1820 analytic = term->dDelta_dTau3(tau, delta);
1821 };
1822 double get_analytic_mcx(shared_ptr<CoolProp::BaseHelmholtzTerm> term, CoolPropDbl tau, CoolPropDbl delta, int ntau, int ndelta) {
1823 try {
1824 using fcn_t = std::function<mcx::MultiComplex<double>(const std::vector<mcx::MultiComplex<double>>&)>;
1825 fcn_t f = [&term](const std::vector<mcx::MultiComplex<double>>& x) { return term->one_mcx(x[0], x[1]); };
1826 return mcx::diff_mcxN(f, std::vector<double>{tau, delta}, std::vector<int>{ntau, ndelta});
1827 } catch (CoolProp::NotImplementedError e) {
1828 return _HUGE;
1829 }
1830 }
1831 double err(double v1, double v2) {
1832 if (std::abs(v2) > 1e-15) {
1833 return std::abs((v1 - v2) / v2);
1834 } else {
1835 return std::abs(v1 - v2);
1836 }
1837 }
1838};
1839
1840std::string terms[] = {"Lead", "LogTau", "IGPower", "PlanckEinstein", "CP0Constant", "CP0PolyT", "Gaussian",
1841 "Lemmon2005", "Power", "SAFT", "NonAnalytic", "Exponential", "GERG2008", "SRK",
1842 "PengRobinson", "XiangDeiters", "GaoB", "GERG2004Cosh", "GERG2004Sinh"};
1843std::string derivs[] = {"dTau", "dTau2", "dTau3", "dDelta", "dDelta2", "dDelta3", "dDelta_dTau",
1844 "dDelta_dTau2", "dDelta2_dTau", "dTau4", "dDelta_dTau3", "dDelta2_dTau2", "dDelta3_dTau", "dDelta4"};
1845std::map<std::string, std::tuple<int, int>> counts = {
1846 {"dTau", {1, 0}},
1847 {"dTau2", {2, 0}},
1848 {"dTau3", {3, 0}},
1849 {"dTau4", {4, 0}},
1850 {"dDelta", {0, 1}},
1851 {"dDelta2", {0, 2}},
1852 {"dDelta3", {0, 3}},
1853 {"dDelta4", {0, 4}},
1854 {"dDelta_dTau", {1, 1}},
1855 {"dDelta_dTau2", {2, 1}},
1856 {"dDelta2_dTau", {1, 2}},
1857 {"dDelta_dTau3", {3, 1}},
1858 {"dDelta2_dTau2", {2, 2}},
1859 {"dDelta3_dTau", {1, 3}},
1860};
1861
1862TEST_CASE_METHOD(HelmholtzConsistencyFixture, "Helmholtz energy derivatives", "[helmholtz]") {
1863 shared_ptr<CoolProp::BaseHelmholtzTerm> term;
1864 std::size_t n = sizeof(terms) / sizeof(terms[0]);
1865 for (std::size_t i = 0; i < n; ++i) {
1866 term = get(terms[i]);
1867 for (std::size_t j = 0; j < sizeof(derivs) / sizeof(derivs[0]); ++j) {
1868 if (terms[i] == "SAFT"
1869 && (derivs[j] == "dTau4" || derivs[j] == "dDelta_dTau3" || derivs[j] == "dDelta2_dTau2" || derivs[j] == "dDelta3_dTau"
1870 || derivs[j] == "dDelta4")) {
1871 continue;
1872 }
1873 double tau = 1.3, delta = 0.9;
1874 call(derivs[j], term, tau, delta, 1e-5);
1875 double alphar = term->base(tau, delta);
1876
1877 // Do calculations with multicomplex, if the one_mcx function has been implemented
1878 auto [ntau, ndelta] = counts.at(derivs[j]);
1879 double numerical_mcx = _HUGE;
1880 double alphar_mcx = _HUGE;
1881 try {
1882 numerical_mcx = get_analytic_mcx(term, tau, delta, ntau, ndelta);
1883 alphar_mcx = term->one_mcx(tau, delta).real();
1884 } catch (std::exception& /* e */) {
1885 //std::cout << e.what() << std::endl;
1886 }
1887 CAPTURE(alphar_mcx);
1888 CAPTURE(numerical_mcx);
1889
1890 CAPTURE(derivs[j]);
1891 CAPTURE(alphar);
1892 double numerical_finitediff = numerical;
1893 CAPTURE(numerical_finitediff);
1894 CAPTURE(analytic);
1895 CAPTURE(terms[i]);
1896 double deriv_tolerance = 1e-9;
1897 if (terms[i] == "GERG2004Cosh" || terms[i] == "GERG2004Sinh" || terms[i] == "CP0PolyT") {
1898 deriv_tolerance = 1e-7; // due to, I think, a loss in precision in the log function of multicomplex
1899 }
1900 double val_tolerance = 1e-14;
1901 if (terms[i] == "CP0PolyT") {
1902 val_tolerance = 1e-10; // due to, I think, a loss in precision in the log function of multicomplex
1903 }
1904 if (terms[i] == "GaoB") {
1905 val_tolerance = 1e-12; // due to, I think, a loss in precision in the log function of multicomplex
1906 }
1907
1908 if (std::isfinite(numerical_mcx)) {
1909 CHECK(err(analytic, numerical_mcx) < deriv_tolerance);
1910 CHECK(err(alphar, alphar_mcx) < val_tolerance);
1911 } else {
1912 CHECK(err(analytic, numerical_finitediff) < deriv_tolerance);
1913 }
1914 }
1915 }
1916}
1917
1918#endif
Generated on Wed Apr 23 2025 03:23:32 for CoolProp by doxygen 1.9.4