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