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